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

Explore the Multi-Taper Method’s unique ability to refine spectral estimates in seismology. Leveraging multiple orthogonal tapers, this approach minimizes variance and mitigates spectral leakage, offering a more reliable analysis of non-stationary seismic signals.

## Introduction

Multi-taper (MT) spectral analysis is a method that is widely used in various scientific fields, including seismology, for the time-frequency analysis of signals. This technique offers a means to assess the power spectrum of a signal in both the time and frequency domains, with the aim of reducing bias and variance in spectral estimates (see Thomson, 1982). It serves as an important tool for the spectral analysis of non-stationary signals.

## Multitaper Spectrum Analysis

The multitaper algorithm shares similarities with other nonparametric spectral estimation techniques. Specifically, an $N$-length time sequence $y_n$ where ($n = 0, 1,…, N$) is multiplied by the kth taper $v_k$ before being Fourier transformed

$Y_k = \sum_{n=0}^{N-1} y_n v_k e^{-2\pi ifn}$

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

$\hat{S} = \frac{\sum_{k=0}^{K-1}w_k|Y_k (f)|}{\sum_{k=0}^{K-1}w_k}$

In this formulation, $\hat{S}$ is a weighted average of the $Y_k$ values, where $w_k$ serves as the weights designed to mitigate spectral leakage while maintaining a low variance in the estimate. The choice of weights $w_k$ can be guided by the eigenvalues of the tapers or through adaptive weighting, as outlined by Thomson (1982). For more advanced applications, the quadratic multitaper estimate, introduced by Prieto et al. (2007), utilizes the second derivative of the spectrum to yield a higher-resolution estimate (See Prieto et. al., 2022).

## Core Concepts for Understanding the Multi-taper Spectral Analysis Method

### Tapers

A taper is essentially a window function that multiplies a signal in the time domain, effectively localizing it both temporally and spectrally. The most straightforward example of a taper is a rectangular window, but there are various other types such as the Hanning or Hamming windows. In multi-taper analysis, multiple tapers are applied to the same data segment.

### Slepian Functions

In multi-taper analysis, the tapers commonly employed are Slepian functions, or Discrete Prolate Spheroidal Sequences (DPSS). These are a set of orthogonal tapers optimized to minimize leakage from other frequencies during the estimation of the power spectrum.

### Spectral Estimation

For each taper, the data are multiplied by the taper in the time domain, and then a Fourier transform is applied to analyze the resulting sequence in the frequency domain. These individual spectral estimates are subsequently aggregated, often by taking an average, to yield a more robust estimate of the power spectrum.

## Analyze a time series to understand the concepts of Multitaper method

• 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 respository 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. (2007). 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)

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 MT method on a 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

Attribute Multi-taper Method Short-Time Fourier Transform (STFT) Wavelet Transform
Resolution Fixed time-frequency resolution Fixed time-frequency resolution Variable time-frequency resolution
Computational Complexity Moderately high Low High
Robustness High (reduced bias and variance) Moderate Moderate to High (for transients)
Interpretability Easier (Fourier-based) Easier (Fourier-based) Complex
Leakage Reduction High (Slepian functions) Moderate (depends on window) Moderate to High
Flexibility Moderate (parameters can be tuned) Low (fixed window) High (many wavelet types)
Applicability Spectral analysis Time-frequency representation Time-scale representation

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.

## Conclusion

Multi-taper spectral analysis offers a robust and versatile method for time-frequency analysis. By using multiple tapers, it provides more reliable and less biased spectral estimates, which are valuable in various applications, including but not limited to seismology.

## References

1. Thomson, D. J. (1982). Spectrum estimation and harmonic analysis, Proc. IEEE 70, no. 9, 1055–1096.
2. Prieto, G. A., R. L. Parker, D. J. Thomson, F. L. Vernon, and R. L. Graham (2007). Reducing the bias of multitaper spectrum estimates, Geophys. J. Int. 171, no. 3, 1269–1281. doi: 10.1111/j.1365-246X.2007.03592.x.
3. Germán A. Prieto; The Multitaper Spectrum Analysis Package in Python. Seismological Research Letters 2022;; 93 (3): 1922–1929. doi: https://doi.org/10.1785/0220210332

Tags:

Categories:

Created on:

Updated on:

#### Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.