Exploring the Multi-Taper Method: A Robust Approach for Time-Frequency Analysis

Introduction

Multi-taper (MT) spectral analysis is a method widely used across scientific fields, including seismology, for the time-frequency analysis of signals. It estimates the power spectrum of a signal while reducing bias and variance in the estimate [1], which makes it an important tool for analyzing non-stationary signals.

Multitaper spectrum analysis

The multitaper algorithm resembles other nonparametric spectral estimators. An $N$-length time sequence $y_n$ (where $n = 0, 1, \ldots, N$) is multiplied by the $k$th taper $v_k$ before being Fourier transformed:

Unlike single-taper methods, the multitaper approach employs $K$ orthogonal tapers to refine the Power Spectral Density (PSD) estimate $\hat{S}$:

Here $\hat{S}$ is a weighted average of the $Y_k$ values, where $w_k$ are weights designed to mitigate spectral leakage while keeping the variance low. The weights can be guided by the eigenvalues of the tapers or by adaptive weighting [1]. For more advanced applications, the quadratic multitaper estimate introduced by Prieto et al. [2] uses the second derivative of the spectrum to yield a higher-resolution estimate (see also the multitaper package paper [3]).

Core concepts for understanding the multi-taper spectral analysis method

Tapers

A taper is a window function that multiplies a signal in the time domain, localizing it temporally and spectrally. The simplest taper is a rectangular window; others include the Hanning or Hamming windows. In multi-taper analysis, multiple tapers are applied to the same data segment.

Slepian functions

The tapers commonly used are Slepian functions, or Discrete Prolate Spheroidal Sequences (DPSS) — a set of orthogonal tapers optimized to minimize leakage from other frequencies during power-spectrum estimation.

Spectral estimation

For each taper, the data are multiplied by the taper in the time domain and Fourier transformed. These individual spectral estimates are then aggregated (often averaged) to yield a more robust estimate of the power spectrum.

Analyze a time series to understand the concepts of the multitaper method

• Download the MSEED data here
• Import Python Libraries:

import numpy as np
from scipy.signal import welch
import matplotlib.pyplot as plt
from scipy.signal.windows import dpss  

• Define DPSS parameters and generate tapers

# DPSS parameters
NW = 4.5  # Time-half bandwidth
K = 6  # Number of tapers

# Generate tapers
tapers = dpss(N, NW, K)

• Visualize the trace data

N = tr.stats.npts  # Number of data points
fs = tr.stats.sampling_rate  # Sampling frequency

# Time vector
t = tr.times()

noisy_signal = tr.data

fig, ax = plt.subplots(1, 1, figsize=(10, 6), sharex=True)
ax.plot(t, noisy_signal)
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.set_title('Waveform')
plt.savefig("multitaper_input_waveform.png", dpi=300, bbox_inches='tight')
plt.show()

• Apply the vanilla multitaper method on the noisy_signal

# Apply tapers
tapered_data = tapers * noisy_signal[np.newaxis, :]

# Initialize for accumulating the power spectrum
average_power_spectrum = np.zeros(N)

# Welch parameters
nperseg = 512  # Length of each segment
noverlap = 256  # Overlap between segments

# Loop over each taper
for k in range(K):
    # Apply Welch Periodogram on each tapered data
    f, Pxx = welch(tapered_data[k, :], fs, nperseg=nperseg, noverlap=noverlap)
    
    # Accumulate the power spectrum
    average_power_spectrum[:len(f)] += 10 * np.log10(Pxx)

# Average the power spectrum
average_power_spectrum[:len(f)] /= K

# Visualize
fig, ax = plt.subplots(1, 1, figsize=(10, 6), sharex=True)
ax.plot(f, average_power_spectrum[:len(f)])
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power/Frequency (dB/Hz)')
ax.set_title('Multitaper Power Spectrum using Welch\'s Method')
ax.set_xlim([0, 10])
plt.savefig("multitaper_output_spectrum.png", dpi=300, bbox_inches='tight')
plt.show()

In the above script, equal weights are given to each taper but this approach will not be most efficient in practice. It is possible to decide the weights of the tapers based on its eigenvalues [tapers, eigvals = dpss(N, NW, 4, return_ratios=True)].

Implementing Thomson’s multitaper method

The multitaper approach previously described incorporates Welch’s method, a classical technique in spectral estimation. In this framework, multiple tapers are computed and subsequently applied to the input signal. The power spectral density (PSD) for each tapered signal is then calculated using Welch’s method. The final PSD is obtained by averaging these individual estimates, essentially forming a straightforward mean of the spectral estimates from each taper.

The “MTSpec” class available in Prieto’s Github repository extends beyond mere spectral estimates to include advanced functionalities like adaptive spectral estimation and quadratic inverse spectrum. In the adaptive spectral estimation approach, an iterative algorithm is deployed to identify optimal weightings for each taper. This algorithm ensures that the weights fall within a range of 0 to 1, thereby normalizing them.

Furthermore, Prieto’s code provides Quadratic Inversion (QI) spectral estimation specifically tailored for seismic data analysis. This method calculates the QI spectrum using the methodology proposed by Prieto et al. [2]. It is designed to mitigate the bias introduced by the curvature characteristics of the multitaper spectrum in seismic data.

## Applying Prieto's method
import multitaper.mtspec as spec


nw = 4.5 #time-bandwidth product (twin x freq_res)/2
kspec = 6 # the standard number of tapers is calculated as K=2×NW−1. However, we could opt to use fewer tapers if computational resources are a concern.
twin = 40. #window length (time resolution is inversely proportional to freq resolution (coming from the uncertainty principle))
# Also, a shorter time window can make the spectrogram more responsive to changes in frequency content, potentially improving frequency resolution
# Decide based on how fast the data is changing
olap = 0.8 #Lower overlap may capture more dynamic changes in your signal
fmax = 25.

freq_res = (nw/twin)*2
dt = tr.stats.delta

t_vals,freq_vals,quad_vals,thomp_vals = spec.spectrogram(noisy_signal,dt,twin = twin, olap = olap, nw = nw, fmax = fmax)

max_qi_vals = np.max(quad_vals, axis=1)
log_psd_db = 10*np.log10(max_qi_vals)

fig, ax = plt.subplots(1, 1, figsize=(10, 6), sharex=True)
ax.plot(freq_vals, log_psd_db, 'b', label='Quad PSD (dB)')
ax.set_xlabel('Frequency [Hz]')
ax.set_xlim([0, 10])
plt.savefig("prieto_method_output_spectrum.png", dpi=300, bbox_inches='tight')

What are the advantages of applying the MT method on time-series data?

1. Reduced Bias and Variance: By combining spectral estimates from multiple tapers, the method reduces both the bias and variance of the spectral estimates, compared to using a single taper.
2. Improved Resolution: The use of multiple tapers allows for better frequency resolution than a single-taper method, especially useful for detecting closely-spaced frequencies.
3. Leakage Reduction: Slepian functions minimize the effect of spectral leakage, which is the contamination of the spectral estimate by power from adjacent frequencies.

Comparison with other methods

When it comes to time-frequency resolution, the Wavelet Transform has an edge, especially for non-stationary signals. Multi-taper methods offer better frequency resolution than STFT but do not have the variable resolution advantage of wavelets. However, Multi-taper methods offer robust spectral estimates, particularly in low signal-to-noise ratio scenarios, surpassing STFT in this regard. Wavelets are more suited for capturing transient phenomena rather than providing robust spectral estimates.

Recap

Without scrolling up — what makes the multi-taper method robust?

• It multiplies the signal by $K$ orthogonal Slepian (DPSS) tapers,
• Fourier-transforms each tapered copy, and averages the spectra,
• which reduces bias and variance and suppresses spectral leakage,
• with adaptive/eigenvalue weights and Prieto’s quadratic-inverse estimate for even higher fidelity.

For noisy, low-SNR spectral estimation it’s a workhorse — trading the wavelet’s variable resolution for statistical robustness.

Where to go next

• The multitaper Python package (Prieto) — adaptive and quadratic-inverse estimators.
• Related posts here: Towards multi-resolution analysis with wavelet transform and Visualizing power spectral density using Obspy.

References

1. Spectrum Estimation and Harmonic Analysis — Thomson, 1982, Proceedings of the IEEE, 70(9), 1055–1096.
2. Reducing the Bias of Multitaper Spectrum Estimates — Prieto, Parker, Thomson, Vernon & Graham, 2007, Geophysical Journal International, 171(3), 1269–1281.
3. The Multitaper Spectrum Analysis Package in Python — Prieto, 2022, Seismological Research Letters, 93(3), 1922–1929.