2.8 Spectral Transforms

Introduction¶

In geosciences, much of the data we deal with — seismic waves, ocean tides, atmospheric pressure variations — comes in the form of time series. These signals often contain information at multiple frequencies, and a simple time-domain analysis may obscure important features of the data.

Spectrograms are useful tools for analyzing and visualizing the frequency content of time-series data. A spectrogram represents the spectral density of a signal as it changes over time. It is especially useful in geoscientific applications where signals are non-stationary, i.e., their frequency content changes over time. Examples include:

• Earthquake seismology: where both low and high-frequency signals are relevant at different stages of an event.
• Atmospheric science: where periodic patterns such as tides and waves have distinct frequency components.
• Remote sensing: where different processes, such as soil moisture fluctuations, can exhibit characteristic frequencies over time.

Why spectrograms matter

• Time-frequency analysis: Spectrograms show how the frequency content of a signal evolves over time, making them suitable for studying non-stationary data.
• Feature extraction: Spectrograms highlight transient events and long-duration patterns, useful for detecting and characterizing geophysical phenomena like earthquakes, landslides, or atmospheric waves.
• Multiscale analysis: Many natural processes operate at different time scales. Spectrograms let us visualize and extract features across these scales.
• Visualization: They offer a compact and interpretable visualization of complex time-series data.

In this section, we transform the data by projecting it onto a basis of functions. The two most used transforms are the Fourier and the wavelet transforms.

Warning. Filtering any data needs to be done carefully. Filtering artifacts can lead to complete misinterpretation. Common pitfalls:

• Not all filters preserve causality: some signals may appear before events and be misinterpreted as precursors.
• Filtering over data gaps brings high-frequency artifacts.
• Edge effects when filtering time series are difficult to mitigate.

The lecture covers several levels and methods for transforming data.

1. Fourier Transforms: 1D [Level 1]
2. Fourier Transforms: 2D [Level 3]
3. Spectrograms [Level 2]
4. Wavelet Transforms [Level 3]

🖥️ Lecture slides — Session 07 (Wed Oct 14)

We first download data: seismograms recorded in the Puget Sound for the M8.2 Chignik, Alaska earthquake of July 29, 2021. We query the FDSN data center with the IRIS client. Note that IRIS data services are now operated by EarthScope, but the client name IRIS still works.

The channel code HHZ denotes a high-gain broadband seismometer (H), sampled at 100 samples per second (H), vertical component (Z). This is a broadband channel: it records ground motion over a wide frequency band, not a single frequency.

We also download the station metadata (the instrument response) and remove the response with remove_response(output="VEL"). This converts the raw digitizer counts to ground velocity in m/s, so the amplitudes have physical units.

1. Fourier Transforms [Level 1]¶

We use the scipy.fft module to transform the two time series (earthquake and noise). The older scipy.fftpack module is legacy and should no longer be used.

The Fourier transform is a decomposition of the time series onto an orthonormal basis of cosine and sine functions. The Fourier transform of a time series $f(t)$ (similarly if the variable is space $x$) is:

$\hat{F}(f) = \int_{-\infty}^\infty f(t) \exp^{-i2\pi ft} dt$

$\hat{F}(f)$ is the complex Fourier value at frequency $f$. The Fourier transform determines what frequencies dominate the time series.

1.1 Nyquist¶

The Fourier transform we use in this class takes a discrete time series of real numbers. If the time series spans $T$ seconds with $N$ regularly spaced samples, the sampling interval is $dt = T/N$. The highest frequency that can be resolved in a discrete time series, called the Nyquist frequency, is limited by $dt$:

$F_{Nyq} = \frac{1}{2 dt}$

Effectively, one cannot constrain signals that vary faster than two time samples. Here $dt = 0.01$ s, so $F_{Nyq} = 50$ Hz.

1.2 Uncertainties¶

• The discrete Fourier Transform yields an approximation of the FT. The shorter the time series, the less accurate the FT. This means that the FT on short time windows is less accurate.
• The FT assumes (and requires) periodicity of the series, meaning that the finite/trimmed time series would repeat in time. To enforce this, we taper the time series so that the first and last points are equal (to zero).

Please see the Obspy documentation to find out about the taper function. Plot the amplitude and phase spectra.

You will note above that the phase values are randomly distributed between $-\pi$ and π. We can check it by showing the distribution of the phase and amplitude spectra.

We can also analyze the spectral characteristics of the noise time series. Below:

1. compute the Fourier transform
2. plot the phase and amplitude spectra
3. plot the distribution of the phase and amplitude values

Overlay their PDFs.

You notice that their statistical differences are in the tails of the distributions. Therefore, statistical metrics such as mean or variance may not be discriminatory, but kurtosis might.

2. 2D Fourier Transforms [Level 3]¶

The 2D Fourier transform is applied to a 2D matrix. It first applies a 1D Fourier transform to every row of the matrix, then applies a 1D Fourier transform to every column of the intermediate matrix.

2D Fourier transforms give the Fourier coefficients that dominate an image. This can be used for filtering the data. Another application is to compress the data by keeping a few coefficients instead of storing the whole image.

We will practice on a synthetic topography-like field. Real topography has a “red” wavenumber spectrum: most of the power sits at long wavelengths (mountain ranges), with progressively less power at short wavelengths (small-scale roughness). The spectral amplitude decays roughly as a power law of the wavenumber $k$, $|\hat{F}(k)| \propto k^{-\beta}$. We can build such a fractal terrain directly in the Fourier domain: draw random phases, scale the amplitudes by $k^{-\beta}$, and inverse transform.

Consider elevation as a 2D data set. We can perform a 2D transform, which gives a spectrum in the spatial dimensions. The axes of the transformed image are wavenumbers (cycles per km). We use fftshift to place the zero wavenumber at the center of the image.

The energy is concentrated near the center of the plot (low wavenumbers, long wavelengths), as designed. Now we will compress the image by keeping only the largest Fourier coefficients.

Now we compare the original 2D data set with the Fourier-compressed data. Keeping 1% of the coefficients is a compression ratio of 100:1, and the reconstruction still captures the large-scale structure of the terrain. The price is the loss of the small-scale roughness, which lives in the discarded high-wavenumber coefficients.

3. Spectrograms [Level 2]¶

In time-dependent and multi-scale problems, it may be interesting to extract data features from the short time Fourier transform (STFT).

The STFT is a Fourier Transform applied to short (overlapping) windows to resolve the frequencies over different times in the series.

The spectrogram Zxx is a transform of the original data. It is common to use spectrograms as input to neural networks as 2D arrays.

4. Continuous Wavelet Transform [Level 3]¶

Introduction¶

The Continuous Wavelet Transform (CWT) is an important tool in geoscientific data analysis, particularly for time-frequency analysis of non-stationary signals. Like the spectrogram, the CWT provides insight into how the frequency content of a signal varies over time. However, the CWT offers better resolution at different frequencies, making it more suitable for analyzing signals with transient or localized frequency changes.

The wavelet transform breaks a signal down into scaled and shifted versions of a small, oscillating function known as the wavelet. This makes the CWT well-suited for geoscientific applications, where many phenomena, such as earthquakes, volcanic eruptions, and weather patterns, can manifest at different scales and frequencies.

Figure: Fourier and wavelet basis functions. Image from this article.

There exist many canonical wavelet families. The difference between families is typically their shape, compactness, and smoothness. Typically, one chooses one family for the specific time series. Wavelets have finite energy and zero mean.

Figure: Families of wavelet basis functions. Image from this article

The wavelet transform is:

$\hat{F}(a,b) = \frac{1}{\sqrt{a}} \int_{-\infty}^\infty f(t) \bar{\Psi} (\frac{t-b}{a}) dt$

where $\bar{\Psi}$ is the mother wavelet scaled by a factor of $a$ and translated/shifted by $b$. In the continuous transform, $a$ and $b$ take continuous values. The Discrete Wavelet Transform is the wavelet transform performed on a finite number of scales and shifts.

The time-scale representation of a time series is a scaleogram. Scales can be converted to pseudo-frequencies: if $f_c$ is the central frequency of the wavelet, the scale is $a$, and the sampling interval is $dt$, then the pseudo-frequency is $f_a = f_c/(a \, dt)$.

Why Continuous Wavelet Transforms are Important¶

• Multiresolution analysis: The CWT captures both low-frequency, long-duration trends and high-frequency, short-duration features. This is particularly valuable in geosciences, where processes occur on different temporal and spatial scales.
• Non-stationary data: Many geoscientific signals are non-stationary, meaning their statistical properties change over time. The CWT reveals these time-varying frequency components.
• Local feature detection: CWT is well suited for identifying localized events, such as seismic waves, landslides, or atmospheric disturbances, by analyzing how the signal changes in both time and frequency.
• Better time-frequency resolution: The wavelet transform provides more precise time and frequency localization than methods like the Fourier transform or spectrogram, especially for short-lived events.

Use in Seismology¶

In seismic data, different types of seismic waves (P-waves, S-waves, and surface waves) occur at different frequencies and durations. By applying CWT, seismologists can detect these different wave phases and their precise arrival times, which matter for earthquake characterization.

Python Example: Continuous Wavelet Transform of Seismic Data¶

We use the PyWavelets (pywt) package. The older scipy.signal.cwt and scipy.signal.morlet2 functions were removed from SciPy, so pywt is now the standard tool. We choose the complex Morlet wavelet 'cmor1.5-1.0' (bandwidth 1.5, center frequency 1.0), a common choice for seismic data because its complex form gives both amplitude and phase.

pywt.cwt works in scales. We pick the frequencies we want to resolve (0.1 to 40 Hz, log-spaced), then convert them to scales with the relation $a = f_c / (f \, dt)$, which is what pywt.frequency2scale computes. To keep the computation light, we apply the CWT to the first 20 minutes of the earthquake record, which contains the P, S, and surface waves.

The scalogram shows the P wave arriving first with high-frequency energy, followed by the S wave and long-period surface waves below 0.1-0.5 Hz. Compare this with the STFT spectrogram above: the wavelet transform resolves the low-frequency surface waves better because its time window adapts to each frequency.

Advantages of CWT in Geosciences:¶

• Time-localized event detection: CWT identifies short-lived geophysical events, such as the arrival of seismic waves during an earthquake.
• Multiscale phenomena: Natural processes like tectonic activity and oceanic tides operate over a broad range of temporal and spatial scales. CWT can analyze signals from these processes by using a wide range of scales.
• Edge detection: The wavelet transform is particularly good at detecting changes or discontinuities in signals, such as the sharp onset of seismic waves or the boundaries of different geophysical layers.

Use Cases in Geoscience:¶

• Earthquake early warning: By applying CWT, geoscientists can more accurately detect the first arrival of seismic waves, which matters for early warning systems.
• Landslide detection: The CWT helps in identifying high-frequency signals indicative of landslides, as these events are often characterized by short bursts of energy.
• Volcanic tremors: CWT can be used to analyze volcanic tremor signals, which are often complex and exhibit both short and long-duration features.

Conclusion¶

The Continuous Wavelet Transform (CWT) is a useful tool for time-frequency analysis in geosciences. Its ability to resolve signals across different scales and frequencies makes it well suited for studying non-stationary processes. Python libraries like pywt, obspy, and matplotlib let geoscientists apply the CWT to their data and extract insight into complex geophysical phenomena.

Time-frequency transforms take computational time in a workflow. Let’s compare the cost of the CWT and the STFT on the same 20-minute segment.

The STFT is much cheaper than the CWT. From these transforms, we can extract similar statistical features.