Introduction
Multitaper (MT) spectral analysis is a method that is widely used in various scientific fields, including seismology, for the timefrequency 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 nonstationary signals.
Multitaper Spectrum Analysis
The multitaper algorithm shares similarities with other nonparametric spectral estimation techniques. Specifically, an length time sequence where () is multiplied by the taper before being Fourier transformed
Unlike singletaper methods, the multitaper approach employs orthogonal tapers to refine the Power Spectral Density (PSD) estimate :
In this formulation, is a weighted average of the values, where serves as the weights designed to mitigate spectral leakage while maintaining a low variance in the estimate. The choice of weights 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 higherresolution estimate (See Prieto et. al., 2022).
Core Concepts for Understanding the Multitaper 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 multitaper analysis, multiple tapers are applied to the same data segment.
Slepian Functions
In multitaper 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
1. Download the MSEED data here
2. Import Python Libraries:
import numpy as np
from scipy.signal import welch
import matplotlib.pyplot as plt
from scipy.signal.windows import dpss
3. Define DPSS parameters and generate tapers
# DPSS parameters
NW = 4.5 # Timehalf bandwidth
K = 6 # Number of tapers
# Generate tapers
tapers = dpss(N, NW, K)
4. Visualize the trace data
from obspy import read
# Load the MSEED file
st = read("path_to_your_mseed_file.mseed")
# Extract the first trace
tr = st[0]
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()
5. 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&amp;amp;amp;#039;s method import multitaper.mtspec as spec nw = 4.5 #timebandwidth 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, &amp;amp;amp;#039;b&amp;amp;amp;#039;, label=&amp;amp;amp;#039;Quad PSD (dB)&amp;amp;amp;#039;) ax.set_xlabel(&amp;amp;amp;#039;Frequency [Hz]&amp;amp;amp;#039;) ax.set_xlim([0, 10]) plt.savefig(&amp;amp;amp;quot;prieto_method_output_spectrum.png&amp;amp;amp;quot;, dpi=300, bbox_inches=&amp;amp;amp;#039;tight&amp;amp;amp;#039;)
What are the advantages of applying MT method on a timeseries data

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.

Improved Resolution: The use of multiple tapers allows for better frequency resolution than a singletaper method, especially useful for detecting closelyspaced frequencies.

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  Multitaper Method  ShortTime Fourier Transform (STFT)  Wavelet Transform 
Resolution  Fixed timefrequency resolution  Fixed timefrequency resolution  Variable timefrequency resolution 
Computational Complexity  Moderately high  Low  High 
Robustness  High (reduced bias and variance)  Moderate  Moderate to High (for transients) 
Interpretability  Easier (Fourierbased)  Easier (Fourierbased)  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  Timefrequency representation  Timescale representation 
When it comes to timefrequency resolution, the Wavelet Transform has an edge, especially for nonstationary signals. Multitaper methods offer better frequency resolution than STFT but do not have the variable resolution advantage of wavelets. However, Multitaper methods offer robust spectral estimates, particularly in low signaltonoise ratio scenarios, surpassing STFT in this regard. Wavelets are more suited for capturing transient phenomena rather than providing robust spectral estimates.
Conclusion
Multitaper spectral analysis offers a robust and versatile method for timefrequency 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.