Wavelet Analysis Applied to Real Geophysical Datasets: A Quick and Easy Approach

Explore the application of wavelet analysis on real-world geophysical datasets, like ENSO sea surface temperature and Indian monsoon rainfall. This post provides a comparison to Fourier analysis and includes Python code for easy replication.
Kumar, U. (2021). Wavelet & Fourier Analysis on the ENSO and monsoon data in Python (Version version 0.1). Zenodo. https://doi.org/10.5281/zenodo.4667892

Wavelet analysis is a powerful tool for time-frequency analysis of signals, providing insight into both when and what frequencies occur in the data. Unlike the Fourier Transform, which decomposes a signal into a sum of sine waves, wavelets offer localized frequency analysis—meaning you can pinpoint both the frequency and the time at which it appears. This is especially valuable when working with non-stationary signals, as is often the case in geophysical studies.

Fourier Transform

The Fourier Transform is one of the most widely used methods in signal processing. It transforms a signal from the time domain to the frequency domain, allowing us to analyze the dominant frequencies in the data. However, one significant limitation of Fourier Transform is that it assumes the signal is stationary, meaning that the frequencies do not change over time. This can be a major drawback when analyzing real-world signals, which often exhibit time-varying behavior.

To address this issue, we can use a Short-Time Fourier Transform (STFT), which applies Fourier Transform within a sliding window across the signal, thus providing time-localized frequency information. However, this method suffers from a trade-off between time and frequency resolution due to the uncertainty principle (Papoulis 1977; Franks 1969), similar to the Heisenberg uncertainty principle in quantum mechanics.

Wavelet Transform

In geophysics, where we often deal with non-stationary signals, Wavelet Transform is a superior alternative to Fourier Transform. The Wavelet Transform preserves high resolution in both time and frequency domains, making it especially useful for analyzing signals that exhibit dynamic, time-varying behavior.

Wavelet analysis works by decomposing a signal into smaller wavelet components, each with its own scale and frequency. Unlike sine waves in the Fourier method, wavelets are localized in both time and frequency, enabling more detailed analysis. This process, called convolution, slides a wavelet across the signal and iterates by adjusting the wavelet’s scale. The resulting wavelet spectrogram or scaleogram provides a clear picture of how the signal’s frequency content evolves over time.

The pseudo-frequency for each wavelet scale is calculated as:

f_{\alpha} = \frac{f_{c}}{\alpha}

where f_{\alpha} is the pseudo-frequency, f_{c} is the central frequency of the mother wavelet, and \alpha is the scaling factor. This means that larger scales (longer wavelets) correspond to smaller frequencies, enabling multi-resolution analysis.

Wavelet Analysis Applied to Sea Surface Temperature & Indian Monsoon Rainfall Data

In this section, I apply the Wavelet Transform to real geophysical datasets: the El Niño-Southern Oscillation (ENSO) sea surface temperature data (1871-1997) and Indian monsoon rainfall data (1871-1995). The wavelet used for this analysis is the complex Morlet wavelet, with a bandwidth of 1.5 and a normalized center frequency of 1.0. This analysis will help visualize how the power of different frequencies in these datasets evolves over time.

import pywt
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.style.use('seaborn')

# Load and process the dataset
dataset = "monsoon.txt"
df_nino = pd.read_table(dataset, skiprows=19, header=None)

N = df_nino.shape[0]
t0 = 1871
dt = 0.25
time = np.arange(0, N) * dt + t0
signal = df_nino.values.squeeze()
signal = signal - np.mean(signal)

# Set wavelet scales
scales = np.arange(1, 128)

# Plot the original signal
def plot_signal(time, signal, figname=None):
    fig, ax = plt.subplots(figsize=(15, 3))
    ax.plot(time, signal, label='signal')
    ax.set_xlim([time[0], time[-1]])
    ax.set_ylabel('Amplitude', fontsize=18)
    ax.set_xlabel('Time', fontsize=18)
    ax.legend()
    plt.savefig(figname if figname else 'signal_plot.png', dpi=300, bbox_inches='tight')
    plt.close()

plot_signal(time, signal)

Signal

Next, we apply Fourier Transform to the signal to compare the time-invariant frequency spectrum.

def get_fft_values(y_values, T, N, f_s):
    f_values = np.linspace(0.0, 1.0/(2.0*T), N//2)
    fft_values_ = np.fft.fft(y_values)
    fft_values = 2.0/N * np.abs(fft_values_[0:N//2])
    return f_values, fft_values

def plot_fft_plus_power(time, signal, figname=None):
    dt = time[1] - time[0]
    N = len(signal)
    fs = 1/dt
    f_values, fft_values = get_fft_values(signal, dt, N, fs)

    fig, ax = plt.subplots(2, 1, figsize=(15, 3))
    ax[0].plot(f_values, fft_values, label='Fourier Transform', color='r')
    ax[1].plot(f_values, fft_values**2, label='FFT Power Spectrum', color='k')
    ax[1].set_xlabel('Frequency [Hz]', fontsize=18)
    ax[1].set_ylabel('Power', fontsize=18)
    plt.savefig(figname if figname else 'fft_plus_power.png', dpi=300, bbox_inches='tight')
    plt.close()

plot_fft_plus_power(time, signal)

Fourier Transform plot

Finally, we use Wavelet Transform to visualize time-frequency localization.

def plot_wavelet(time, signal, scales, waveletname='cmor1.5-1.0', cmap=plt.cm.seismic, figname=None):
    dt = time[1] - time[0]
    coefficients, frequencies = pywt.cwt(signal, scales, waveletname, dt)
    power = np.abs(coefficients)**2
    period = 1. / frequencies

    fig, ax = plt.subplots(figsize=(15, 10))
    im = ax.contourf(time, np.log2(period), np.log2(power), cmap=cmap)
    ax.set_title('Wavelet Transform (Power Spectrum)', fontsize=20)
    ax.set_xlabel('Time', fontsize=18)
    ax.set_ylabel('Period (years)', fontsize=18)
    plt.colorbar(im, ax=ax)
    plt.savefig(figname if figname else f'wavelet_{waveletname}.png', dpi=300, bbox_inches='tight')
    plt.close()

plot_wavelet(time, signal, scales)

Final Observations

For the ENSO dataset, the analysis reveals significant power concentration in the 2–8 year period range, particularly before 1920. The power decreases thereafter, with the wavelet spectrogram showing a transition to shorter periods as time progresses.

In the Indian monsoon rainfall dataset, the power is more evenly distributed across periods but exhibits a similar shift to shorter periods over time. The Wavelet Transform effectively highlights these time-varying features that remain obscured in the Fourier Transform.

The wavelet analysis for the ENSO (a-d) and Indian monsoon rainfall (e-h) time series
(a) and (e) shows the time series in °C and mm, respectively for ENSO and Indian monsoon. (b) and (f) shows the Fourier Transform, (c) and (g) shows the Power spectrum of ENSO and Indian monsoon respectively.

Please cite this work if you use this code:

Kumar, U., Chao, B. F., & Chang, E. T.-Y. Y. (2020). What Causes the Common‐Mode Error in Array GPS Displacement Fields: Case Study for Taiwan in Relation to Atmospheric Mass Loading. Earth and Space Science, 0–2. https://doi.org/10.1029/2020ea001159

References

Yang, Y., & Song, X. (2023). Multidecadal variation of the Earth’s inner-core rotation. Nature Geoscience, 16(2), 182–187. https://doi.org/10.1038/s41561-022-01112-z
Kumar, U., Chao, B. F., & Chang, E. T.-Y. Y. (2020). What Causes the Common‐Mode Error in Array GPS Displacement Fields: Case Study for Taiwan in Relation to Atmospheric Mass Loading. Earth and Space Science, 0–2. https://doi.org/10.1029/2020ea001159
Franks, L. E. (1969). Signal theory.
Papoulis, A. (1977). Signal analysis (Vol. 191). McGraw-Hill New York.
Torrence, C., & Compo, G. P. (1998). A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 79(1), 61–78.
Utpal Kumar
Utpal Kumar

Geophysicist | Geodesist | Seismologist | Open-source Developer
I am a geophysicist with a background in computational geophysics, currently working as a postdoctoral researcher at UC Berkeley. My research focuses on seismic data analysis, structural health monitoring, and understanding deep Earth structures. I have had the opportunity to work on diverse projects, from investigating building characteristics using smartphone data to developing 3D models of the Earth's mantle beneath the Yellowstone hotspot.

In addition to my research, I have experience in cloud computing, high-performance computing, and single-board computers, which I have applied in various projects. This includes working with platforms like AWS, Docker, and Kubernetes, as well as supercomputing environments such as STAMPEDE2, ANVIL, Savio and PERLMUTTER (and CORI). My work involves developing innovative solutions for structural health monitoring and advancing real-time seismic response analysis. I am committed to applying these skills to further research in computational seismology and structural health monitoring.

Articles: 30

Leave a Reply

Your email address will not be published. Required fields are marked *