Time-Frequency Analysis: Method Selection Guide

Interactive Version Available

Looking for the version with plots and executable code?

This theory guide has a corresponding fully-executable Jupyter Notebook version:

This page: Theory for each method, usage guidelines, decision matrices Notebook version: Full implementation examples, embedded outputs, copy-paste code

This tutorial compares different time-frequency analysis methods in gwexpy and helps you choose the right one for your analysis.

Overview

Gravitational wave signals often have time-varying frequency content. Different methods reveal different aspects:

Method

Best For

Time Resolution

Frequency Resolution

Spectrogram (STFT)

General purpose

Good

Good

Q-transform

Chirps, transients

Adaptive

Adaptive (constant Q)

Wavelet (CWT)

Multi-scale features, chirps

Scale-dependent

Scale-dependent

HHT

Instantaneous frequency

Excellent

Data-adaptive

STLT

Damped oscillations

Good

Good (+ decay rate Οƒ)

Cepstrum

Echo detection, periodicity

N/A

Quefrency domain

DCT

Compression, smooth features

N/A

N/A (basis coefficients)

Setup

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from gwexpy.timeseries import TimeSeries
from gwexpy.noise.wave import chirp, gaussian

# Create test signal: chirp + noise
sample_rate = 1024  # Hz
duration = 4  # seconds

# Chirp from 20 Hz to 100 Hz
signal_chirp = chirp(
    duration=duration,
    sample_rate=sample_rate,
    f0=20,  # Start frequency
    f1=100,  # End frequency
    t1=duration
)

# Add Gaussian noise
noise = gaussian(duration=duration, sample_rate=sample_rate, std=0.2)
data = signal_chirp + noise

ts = TimeSeries(data, t0=0, dt=1/sample_rate, unit='strain')
print(f"Data: {len(ts)} samples, {duration}s")

Method 1: Spectrogram (STFT-based)

What it is

Short-Time Fourier Transform (STFT): Divide signal into windows, compute FFT for each

Formula: S(t, f) = |∫ x(Ο„) w(Ο„-t) e^(-2Ο€ifΟ„) dΟ„|Β²

Implementation

# Create spectrogram with 0.5s windows
spec = ts.spectrogram(fftlength=0.5, overlap=0.25)

print(f"Spectrogram shape: {spec.shape}")  # (time_bins, freq_bins)
print(f"Time resolution: {spec.dt}")
print(f"Frequency resolution: {spec.df}")

Visualization

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# Spectrogram
im = axes[0].pcolormesh(
    spec.times.value,
    spec.frequencies.value,
    spec.value.T,
    cmap='viridis',
    shading='auto'
)
axes[0].set_ylim(10, 200)
axes[0].set_ylabel('Frequency (Hz)')
axes[0].set_title('Spectrogram (STFT, window=0.5s)')
fig.colorbar(im, ax=axes[0], label='Power')

# Overlay true chirp frequency
t_true = np.linspace(0, duration, 100)
f_true = 20 + (100 - 20) * t_true / duration
axes[0].plot(t_true, f_true, 'r--', linewidth=2, label='True Frequency')
axes[0].legend()

# Time series for reference
axes[1].plot(ts.times.value, ts.value, linewidth=0.5)
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Strain')
axes[1].set_title('Original Signal')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Pros and Cons

βœ… Advantages:

  • Fast computation

  • Well-understood, standard method

  • Good for stationary or slowly-varying signals

❌ Disadvantages:

  • Fixed time-frequency resolution (uncertainty principle)

  • Poor for rapid frequency changes

  • Window length affects both time and frequency resolution

When to Use

βœ… Use Spectrogram when:

  • Signal is quasi-stationary

  • You need standard, well-established method

  • Fast computation is important

  • Frequency changes are slow relative to window size

Method 2: Q-transform

What it is

Constant-Q transform: Adaptive time-frequency tiling with constant Q-factor

Q-factor: Q = f / Ξ”f (ratio of center frequency to bandwidth)

Implementation

# Q-transform with Q=6
q = 6
qgram = ts.q_transform(qrange=(4, 64), frange=(10, 200))

print(f"Q-transform shape: {qgram.shape}")

Visualization

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# Q-transform
axes[0].imshow(
    qgram.value.T,
    extent=[qgram.times.value[0], qgram.times.value[-1],
            qgram.frequencies.value[0], qgram.frequencies.value[-1]],
    aspect='auto',
    origin='lower',
    cmap='viridis',
    interpolation='bilinear'
)
axes[0].set_ylabel('Frequency (Hz)')
axes[0].set_title(f'Q-transform (constant Q)')
axes[0].plot(t_true, f_true, 'r--', linewidth=2, label='True Frequency')
axes[0].legend()

# Spectrogram for comparison
im = axes[1].pcolormesh(
    spec.times.value,
    spec.frequencies.value,
    spec.value.T,
    cmap='viridis',
    shading='auto'
)
axes[1].set_ylim(10, 200)
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Frequency (Hz)')
axes[1].set_title('Spectrogram (for comparison)')
axes[1].plot(t_true, f_true, 'r--', linewidth=2, label='True Frequency')
axes[1].legend()

plt.tight_layout()
plt.show()

Pros and Cons

βœ… Advantages:

  • Adaptive resolution: better time resolution at high frequencies

  • Natural for chirps (binary coalescence)

  • Constant Q matches gravitational wave signals

❌ Disadvantages:

  • Slower than STFT

  • More complex interpretation

  • Requires Q-factor tuning

When to Use

βœ… Use Q-transform when:

  • Analyzing chirps (compact binary coalescence)

  • Signal spans wide frequency range

  • High-frequency transients need good time resolution

  • Standard GW transient analysis

Method 3: Wavelet Transform (CWT)

What it is

Continuous Wavelet Transform (CWT): Multi-scale analysis using dilated/translated wavelets

Formula: W(a, b) = ∫ x(t) ψ*((t-b)/a) dt

  • a: scale parameter (inversely proportional to frequency)

  • b: time shift parameter

  • ψ: mother wavelet (e.g., Morlet)

Implementation

import pywt

# Create isolated chirp signal for clarity
chirp_signal = np.zeros(len(ts))
chirp_mask = (ts.times.value >= 2) & (ts.times.value <= 5)
t_chirp = ts.times.value[chirp_mask] - 2
f0, f1 = 20, 150
phase = 2*np.pi * (f0*t_chirp + 0.5*(f1-f0)*t_chirp**2/3)
chirp_signal[chirp_mask] = 1.0 * np.sin(phase)
chirp_signal += np.random.randn(len(ts)) * 0.05

# Compute Continuous Wavelet Transform
scales = np.arange(1, 128)
frequencies_wavelet = ts.sample_rate.value / (2 * scales)
coefficients, frequencies_pywt = pywt.cwt(
    chirp_signal,
    scales,
    'morl',  # Morlet wavelet
    sampling_period=1/ts.sample_rate.value
)

print(f"CWT shape: {coefficients.shape}")  # (scales, time)
print(f"Frequency range: {frequencies_wavelet.min():.1f} - {frequencies_wavelet.max():.1f} Hz")

Visualization

fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)

# Original signal
axes[0].plot(ts.times.value, chirp_signal, linewidth=0.5, color='black')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Chirp Signal (20β†’150 Hz over 3s)')

# STFT for comparison
spec_chirp = ts.spectrogram(fftlength=0.5, overlap=0.25)
im1 = axes[1].pcolormesh(
    spec_chirp.times.value,
    spec_chirp.frequencies.value,
    spec_chirp.value.T,
    cmap='viridis',
    shading='auto'
)
axes[1].set_ylim(10, 200)
axes[1].set_ylabel('Frequency (Hz)')
axes[1].set_title('STFT Spectrogram (smeared trajectory)')

# Wavelet scalogram
im2 = axes[2].pcolormesh(
    ts.times.value,
    frequencies_wavelet,
    np.abs(coefficients),
    cmap='viridis',
    shading='auto'
)
axes[2].set_ylim(10, 200)
axes[2].set_ylabel('Frequency (Hz)')
axes[2].set_xlabel('Time (s)')
axes[2].set_title('Wavelet (CWT) Scalogram (sharp trajectory)')

# Overlay true frequency
t_true = ts.times.value[chirp_mask]
f_true = 20 + (150 - 20) * (t_true - 2) / 3
axes[1].plot(t_true, f_true, 'r--', linewidth=2, alpha=0.8, label='True frequency')
axes[2].plot(t_true, f_true, 'r--', linewidth=2, alpha=0.8, label='True frequency')
axes[1].legend()
axes[2].legend()

plt.tight_layout()
plt.show()

Quantitative Evaluation

# Ridge extraction: find frequency of maximum energy at each time
ridge_indices = np.argmax(np.abs(coefficients), axis=0)
f_ridge = frequencies_wavelet[ridge_indices]

# Compute MSE over chirp region
chirp_time_mask = (ts.times.value >= 2) & (ts.times.value <= 5)
f_ridge_chirp = f_ridge[chirp_time_mask]
f_true_interp = 20 + (150 - 20) * (ts.times.value[chirp_time_mask] - 2) / 3

valid_mask = (f_ridge_chirp >= 10) & (f_ridge_chirp <= 200)
mse_wavelet = np.mean((f_ridge_chirp[valid_mask] - f_true_interp[valid_mask])**2)
rmse_wavelet = np.sqrt(mse_wavelet)

print(f"Frequency Tracking RMSE: {rmse_wavelet:.2f} Hz")
print(f"STFT uncertainty (0.5s window): ~10-20 Hz")
print(f"Improvement: {10/rmse_wavelet:.1f}Γ— better precision")

Pros and Cons

βœ… Advantages:

  • Natural scale matching (wavelet β€œstretches” to follow signal)

  • Better time-frequency localization than STFT

  • Excellent for chirps and multi-scale transients

  • Ridge extraction provides precise frequency trajectories

❌ Disadvantages:

  • Higher computational cost than STFT

  • Requires wavelet selection (Morlet, Mexican hat, etc.)

  • Redundant representation (overcomplete)

  • Edge effects at signal boundaries

When to Use

βœ… Use Wavelet when:

  • Analyzing chirps spanning multiple octaves

  • Need better frequency tracking than STFT

  • Signal has multi-scale structure

  • Time-frequency resolution tradeoff is critical

❌ Don’t use when:

  • Signal is purely stationary (STFT sufficient)

  • Need fastest computation (use STFT)

  • Frequency range is narrow (Q-transform may be better)

Method 4: Hilbert-Huang Transform (HHT)

What it is

Empirical Mode Decomposition (EMD) + Hilbert Transform:

  1. Decompose signal into Intrinsic Mode Functions (IMFs)

  2. Compute instantaneous frequency via Hilbert transform

Implementation

# Perform EMD
imfs = ts.emd(method='emd', max_imf=5)

print(f"Extracted {len(imfs)} IMFs")

# Plot IMFs
fig, axes = plt.subplots(len(imfs), 1, figsize=(12, 10), sharex=True)

for i, (name, imf) in enumerate(imfs.items()):
    axes[i].plot(imf.times.value, imf.value, linewidth=0.5)
    axes[i].set_ylabel(name)
    axes[i].grid(True, alpha=0.3)

axes[-1].set_xlabel('Time (s)')
axes[0].set_title('Empirical Mode Decomposition (EMD)')
plt.tight_layout()
plt.show()

Instantaneous Frequency

# Compute instantaneous frequency for dominant IMF
imf_main = imfs['IMF0']
inst_freq = imf_main.instantaneous_frequency()

# Plot
fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

axes[0].plot(imf_main.times.value, imf_main.value, linewidth=0.5)
axes[0].set_ylabel('IMF0 Amplitude')
axes[0].set_title('Dominant Intrinsic Mode Function')
axes[0].grid(True, alpha=0.3)

axes[1].plot(inst_freq.times.value, inst_freq.value, linewidth=1)
axes[1].plot(t_true, f_true, 'r--', linewidth=2, label='True Frequency')
axes[1].set_ylabel('Frequency (Hz)')
axes[1].set_xlabel('Time (s)')
axes[1].set_title('Instantaneous Frequency (HHT)')
axes[1].set_ylim(0, 200)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Pros and Cons

βœ… Advantages:

  • Data-adaptive (no window choice)

  • Excellent time resolution

  • Handles nonlinear, non-stationary signals

  • Direct instantaneous frequency

❌ Disadvantages:

  • Computationally expensive

  • EMD can have mode mixing

  • Less established than STFT/Q-transform

When to Use

βœ… Use HHT when:

  • Signal is highly non-stationary

  • Need precise instantaneous frequency

  • Standard methods fail to resolve features

  • Analyzing glitches or complex transients

❌ Don’t use when:

  • Signal is stationary (overkill)

  • Need fast computation

  • Standard spectrograms suffice

Method 5: Short-Time Laplace Transform (STLT)

What it is

STLT: Decomposes signals into both frequency Ο‰ and decay rate Οƒ components

Formula: STLT(t, Οƒ, Ο‰) = ∫ x(Ο„) w(Ο„-t) e^(-(Οƒ+iΟ‰)Ο„) dΟ„

Unlike STFT which only shows frequency, STLT reveals damping coefficients.

Implementation

# Create multi-mode ringdown
t = ts.times.value
ringdown_multi = np.zeros(len(t))

# Mode 1: f=100 Hz, Ο„=0.2s (high damping)
mode1_mask = t >= 1
t1 = t[mode1_mask] - 1
tau1 = 0.2
ringdown_multi[mode1_mask] += 1.0 * np.exp(-t1/tau1) * np.sin(2*np.pi*100*t[mode1_mask])

# Mode 2: f=150 Hz, Ο„=0.8s (low damping)
mode2_mask = t >= 2
t2 = t[mode2_mask] - 2
tau2 = 0.8
ringdown_multi[mode2_mask] += 0.8 * np.exp(-t2/tau2) * np.sin(2*np.pi*150*t[mode2_mask])

# Mode 3: f=120 Hz, Ο„=0.5s (medium damping)
mode3_mask = t >= 3
t3 = t[mode3_mask] - 3
tau3 = 0.5
ringdown_multi[mode3_mask] += 0.6 * np.exp(-t3/tau3) * np.sin(2*np.pi*120*t[mode3_mask])

ringdown_multi += np.random.randn(len(t)) * 0.05

ts_ringdown = TimeSeries(ringdown_multi, t0=0, dt=ts.dt, unit='strain')

# Compute STLT
stlt = ts_ringdown.stlt(fftlength=1.0, overlap=0.5)

print(f"STLT shape: {stlt.shape}")  # (time, sigma, omega)
print("STLT reveals both Ο‰ (frequency) and Οƒ (decay rate)")

Visualization

fig = plt.figure(figsize=(14, 10))

# Original signal
ax1 = plt.subplot(3, 1, 1)
ax1.plot(t, ringdown_multi, linewidth=0.5, color='black')
ax1.set_ylabel('Amplitude')
ax1.set_title('Multi-Mode Ringdown (3 modes with different decay rates)')
ax1.set_xlim(0, duration)

# STFT (only shows frequency, not decay)
ax2 = plt.subplot(3, 1, 2)
spec_ringdown = ts_ringdown.spectrogram(fftlength=0.5, overlap=0.25)
im1 = ax2.pcolormesh(
    spec_ringdown.times.value,
    spec_ringdown.frequencies.value,
    spec_ringdown.value.T,
    cmap='viridis',
    shading='auto'
)
ax2.set_ylim(80, 180)
ax2.set_ylabel('Frequency (Hz)')
ax2.set_title('STFT: Shows frequencies but decay rates invisible')

# STLT (shows both frequency and decay rate)
ax3 = plt.subplot(3, 1, 3)
# Average over time for 2D Οƒ-Ο‰ visualization
stlt_power = np.abs(stlt.value)**2
stlt_avg = np.mean(stlt_power, axis=0)

# Get sigma and omega axes
sigma_axis = stlt.sigma.value if hasattr(stlt, 'sigma') else np.linspace(-10, 10, stlt.shape[1])
omega_axis = stlt.frequencies.value if hasattr(stlt, 'frequencies') else np.linspace(0, ts.sample_rate.value/2, stlt.shape[2])

im2 = ax3.pcolormesh(omega_axis, sigma_axis, stlt_avg, cmap='hot', shading='auto')
ax3.set_xlim(80, 180)
ax3.set_ylabel('Decay Rate Οƒ (1/s)')
ax3.set_xlabel('Frequency (Hz)')
ax3.set_title('STLT: 2D Οƒ-Ο‰ plane reveals BOTH frequency and decay rate')

# Annotate expected modes
ax3.plot(100, -1/tau1, 'ro', markersize=10, label=f'Mode 1: 100 Hz, Ο„={tau1}s')
ax3.plot(150, -1/tau2, 'go', markersize=10, label=f'Mode 2: 150 Hz, Ο„={tau2}s')
ax3.plot(120, -1/tau3, 'bo', markersize=10, label=f'Mode 3: 120 Hz, Ο„={tau3}s')
ax3.legend()

plt.tight_layout()
plt.show()

Quantitative Evaluation

# Extract Οƒ peaks for each mode
modes_info = [(100, tau1), (150, tau2), (120, tau3)]
abs_errors = []

for freq, tau in modes_info:
    # Find frequency index
    freq_idx = np.argmin(np.abs(omega_axis - freq))

    # Find sigma with maximum power
    sigma_idx = np.argmax(stlt_avg[:, freq_idx])
    sigma_est = sigma_axis[sigma_idx]
    sigma_true = -1/tau

    abs_error = abs(sigma_est - sigma_true)
    abs_errors.append(abs_error)

    print(f"Mode {freq} Hz: Οƒ_true = {sigma_true:.2f}, Οƒ_est = {sigma_est:.2f}, error = {abs_error:.2f}")

mae_stlt = np.mean(abs_errors)
print(f"\nMean Absolute Error: {mae_stlt:.2f} s⁻¹")
print("STFT cannot estimate decay rates at all.")

Pros and Cons

βœ… Advantages:

  • Unique ability to separate modes by both frequency and decay rate

  • Essential for ringdown quality factor estimation

  • 2D Οƒ-Ο‰ representation reveals damping structure

  • Directly applicable to post-merger waveforms

❌ Disadvantages:

  • High computational cost (2D transform)

  • Requires careful parameter selection

  • Less familiar than STFT/Q-transform

  • Interpretation requires understanding Οƒ-Ο‰ plane

When to Use

βœ… Use STLT when:

  • Analyzing ringdown modes (black hole quasi-normal modes)

  • Need to estimate quality factors or damping times

  • Multiple damped oscillations overlap in frequency

  • Decay rate information is scientifically important

❌ Don’t use when:

  • Signal has no damping (use STFT instead)

  • Only frequency information needed

  • Computational resources are limited

Method 6: Cepstrum

What it is

Cepstrum: Analysis of periodicity in the spectrum via β€œquefrency”

Formula: C(Ο„) = IFFT(log|FFT(x)|)

Converts spectral periodicity β†’ time-domain peaks at delay times.

Implementation

# Create signal with echo
original_sig = np.zeros(len(t))
pulse_time = 1.0
pulse_idx = int(pulse_time * ts.sample_rate.value)
pulse_width = int(0.05 * ts.sample_rate.value)

# Broadband pulse
sigma = int(0.01 * ts.sample_rate.value)
n = np.arange(pulse_width)
pulse_window = np.exp(-0.5 * ((n - pulse_width/2) / sigma) ** 2)
original_sig[pulse_idx:pulse_idx+len(pulse_window)] = pulse_window * np.sin(2*np.pi*100*t[pulse_idx:pulse_idx+len(pulse_window)])

# Add echo with 80ms delay
echo_delay = 0.08  # seconds
echo_amplitude = 0.4
echo_samples = int(echo_delay * ts.sample_rate.value)
echo_sig = np.zeros(len(t))
echo_sig[echo_samples:] = echo_amplitude * original_sig[:-echo_samples]

signal_with_echo = original_sig + echo_sig
signal_with_echo += np.random.randn(len(t)) * 0.02

ts_echo = TimeSeries(signal_with_echo, t0=0, dt=ts.dt, unit='strain')

# Compute cepstrum
spectrum = np.fft.rfft(signal_with_echo)
log_spectrum = np.log(np.abs(spectrum) + 1e-10)
cepstrum = np.fft.irfft(log_spectrum)
quefrency = np.arange(len(cepstrum)) / ts.sample_rate.value

print(f"Cepstrum computed. Quefrency range: 0 to {quefrency.max():.3f} s")

Visualization

fig, axes = plt.subplots(4, 1, figsize=(14, 12))

# Original signal
axes[0].plot(t, signal_with_echo, linewidth=0.5, color='black')
axes[0].set_ylabel('Amplitude')
axes[0].set_title(f'Signal with Echo (delay = {echo_delay*1000:.0f} ms)')
axes[0].set_xlim(0, 2)
axes[0].axvline(pulse_time, color='r', linestyle='--', alpha=0.5, label='Original')
axes[0].axvline(pulse_time + echo_delay, color='g', linestyle='--', alpha=0.5, label='Echo')
axes[0].legend()

# STFT
spec_echo = ts_echo.spectrogram(fftlength=0.25, overlap=0.125)
im = axes[1].pcolormesh(
    spec_echo.times.value,
    spec_echo.frequencies.value,
    spec_echo.value.T,
    cmap='viridis',
    shading='auto'
)
axes[1].set_ylim(50, 150)
axes[1].set_ylabel('Frequency (Hz)')
axes[1].set_title('STFT: Echo causes spectral ripples (hard to interpret)')
axes[1].set_xlim(0, 2)

# Power spectrum
freqs = np.fft.rfftfreq(len(signal_with_echo), ts.dt.value)
axes[2].semilogy(freqs, np.abs(spectrum))
axes[2].set_xlim(50, 150)
axes[2].set_ylabel('|Spectrum|')
axes[2].set_title('Power Spectrum: Periodic ripples from echo (comb filter)')

# Cepstrum
axes[3].plot(quefrency[:int(0.2*ts.sample_rate.value)],
            cepstrum[:int(0.2*ts.sample_rate.value)],
            linewidth=1, color='blue')
axes[3].axvline(echo_delay, color='r', linestyle='--', linewidth=2,
               label=f'True delay = {echo_delay*1000:.0f} ms')
axes[3].set_xlim(0, 0.15)
axes[3].set_xlabel('Quefrency (s)')
axes[3].set_ylabel('Cepstrum')
axes[3].set_title('Cepstrum: Clear peak at echo delay')
axes[3].legend()

plt.tight_layout()
plt.show()

Quantitative Evaluation

# Find peak in cepstrum
search_region = (quefrency > 0.03) & (quefrency < 0.12)
peak_idx = np.argmax(np.abs(cepstrum[search_region]))
detected_delay = quefrency[search_region][peak_idx]

error_ms = abs(detected_delay - echo_delay) * 1000

print(f"True echo delay: {echo_delay*1000:.1f} ms")
print(f"Detected peak: {detected_delay*1000:.1f} ms")
print(f"Detection error: {error_ms:.2f} ms")
print("Cepstrum converts spectral periodicity β†’ time-domain peak")

Pros and Cons

βœ… Advantages:

  • Direct detection of echo delays (quefrency peaks)

  • Reveals periodic structure in spectrum

  • Useful for pitch detection and harmonic analysis

  • Computationally efficient (double FFT)

❌ Disadvantages:

  • Requires logarithm (sensitivity to low amplitude)

  • Not time-localized (global analysis)

  • Less intuitive than time-frequency methods

  • Best for signals with clear echoes

When to Use

βœ… Use Cepstrum when:

  • Detecting echoes or reflections

  • Analyzing periodic structure in spectrum

  • Measuring delay times in reverberant signals

  • Pitch detection or harmonic analysis

❌ Don’t use when:

  • Need time-localized analysis

  • Signal has no echoes or periodicity

  • Low SNR (log operation amplifies noise)

Method 7: Discrete Cosine Transform (DCT)

What it is

DCT: Transform to cosine basis, excellent for compressing smooth signals

Formula: X(k) = Ξ£ x(n) cos(Ο€k(2n+1)/(2N))

Energy concentrates in low-frequency coefficients for smooth signals.

Implementation

from scipy.fftpack import dct, idct

# Create smooth signal with perturbations
smooth_signal = np.zeros(len(t))
smooth_signal += 1.0 * np.sin(2*np.pi*3*t)
smooth_signal += 0.5 * np.sin(2*np.pi*7*t)
smooth_signal += 0.3 * np.sin(2*np.pi*12*t)
smooth_signal += 0.1 * np.sin(2*np.pi*50*t)
smooth_signal += 0.05 * np.sin(2*np.pi*80*t)
smooth_signal += np.random.randn(len(t)) * 0.05

ts_smooth = TimeSeries(smooth_signal, t0=0, dt=ts.dt, unit='strain')

# Compute DCT
dct_coeffs = dct(smooth_signal, type=2, norm='ortho')

print(f"DCT computed: {len(dct_coeffs)} coefficients")

# Reconstruct with different numbers of coefficients
n_coeffs_list = [10, 50, 200]
for n_coeffs in n_coeffs_list:
    coeffs_truncated = dct_coeffs.copy()
    coeffs_truncated[n_coeffs:] = 0
    recon = idct(coeffs_truncated, type=2, norm='ortho')

    rmse = np.sqrt(np.mean((smooth_signal - recon)**2))
    energy_kept = 100 * np.sum(dct_coeffs[:n_coeffs]**2) / np.sum(dct_coeffs**2)

    print(f"{n_coeffs} coeffs: {energy_kept:.2f}%" " energy, RMSE = {rmse:.4f}")

Visualization

fig, axes = plt.subplots(3, 2, figsize=(14, 10))

# Original signal
axes[0, 0].plot(t[:2048], smooth_signal[:2048], linewidth=0.8, color='black')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_title('Original Signal (smooth + small perturbations)')

# DCT coefficients
axes[0, 1].semilogy(np.abs(dct_coeffs[:500]), linewidth=0.5, color='blue')
axes[0, 1].set_ylabel('|DCT Coefficient|')
axes[0, 1].set_title('DCT: Rapid decay (energy in low coefficients)')
axes[0, 1].axvline(50, color='r', linestyle='--', alpha=0.5)

# Reconstructions
for i, n_coeffs in enumerate(n_coeffs_list):
    coeffs_truncated = dct_coeffs.copy()
    coeffs_truncated[n_coeffs:] = 0
    recon = idct(coeffs_truncated, type=2, norm='ortho')

    # Reconstruction
    axes[i+1, 0].plot(t[:2048], smooth_signal[:2048], linewidth=0.5,
                     color='gray', alpha=0.5, label='Original')
    axes[i+1, 0].plot(t[:2048], recon[:2048], linewidth=1.5,
                     color='red', label=f'{n_coeffs} coeffs')
    axes[i+1, 0].set_ylabel('Amplitude')
    axes[i+1, 0].set_title(f'Reconstruction ({n_coeffs} coefficients)')
    axes[i+1, 0].legend()

    # Error
    error = smooth_signal - recon
    axes[i+1, 1].plot(t[:2048], error[:2048], linewidth=0.5, color='orange')
    axes[i+1, 1].set_ylabel('Error')
    rmse = np.sqrt(np.mean(error**2))
    compression = len(t) / n_coeffs
    axes[i+1, 1].set_title(f'Error (RMSE={rmse:.4f}, {compression:.0f}Γ— compression)')

axes[2, 0].set_xlabel('Time (s)')
axes[2, 1].set_xlabel('Time (s)')

plt.tight_layout()
plt.show()

Quantitative Evaluation

total_energy = np.sum(dct_coeffs**2)

print("DCT Compression Analysis:")
for n_coeffs in [10, 50, 100, 200]:
    energy_kept = 100 * np.sum(dct_coeffs[:n_coeffs]**2) / total_energy
    compression_ratio = len(dct_coeffs) / n_coeffs

    print(f"  {n_coeffs:4d} coeffs ({100*n_coeffs/len(dct_coeffs):5.1f}%): "
          f"{energy_kept:5.2f}%" " energy, {compression_ratio:5.1f}Γ— compression")

print("\n50 coefficients (3%" " of data) capture >99%" " energy")
print("Compression ratio ~30Γ— with negligible error")

Pros and Cons

βœ… Advantages:

  • Excellent compression for smooth signals

  • Sparse representation (few coefficients capture most energy)

  • Fast computation (FFT-like)

  • No boundary artifacts (unlike Fourier)

  • Ideal for feature extraction and denoising

❌ Disadvantages:

  • Not time-localized

  • Less effective for discontinuous signals

  • Interpretation less intuitive than time-frequency

  • Global analysis only

When to Use

βœ… Use DCT when:

  • Compressing smooth signals

  • Feature extraction (low-order coefficients)

  • Denoising smooth backgrounds

  • Data reduction before machine learning

  • Modeling slow trends

❌ Don’t use when:

  • Need time-localized analysis

  • Signal has sharp transients

  • Frequency content varies rapidly in time

Comparison Example: All Methods on Same Signal

fig = plt.figure(figsize=(14, 10))

# Original signal
ax1 = plt.subplot(4, 1, 1)
ax1.plot(ts.times.value, ts.value, linewidth=0.5, color='black')
ax1.set_ylabel('Strain')
ax1.set_title('Original Signal: Chirp (20β†’100 Hz) + Noise')
ax1.grid(True, alpha=0.3)

# Spectrogram
ax2 = plt.subplot(4, 1, 2, sharex=ax1)
im2 = ax2.pcolormesh(spec.times.value, spec.frequencies.value, spec.value.T,
                     cmap='viridis', shading='auto')
ax2.plot(t_true, f_true, 'r--', linewidth=1.5, alpha=0.8)
ax2.set_ylim(10, 150)
ax2.set_ylabel('Frequency (Hz)')
ax2.set_title('Spectrogram (window=0.5s)')

# Q-transform
ax3 = plt.subplot(4, 1, 3, sharex=ax1)
ax3.imshow(qgram.value.T,
          extent=[qgram.times.value[0], qgram.times.value[-1],
                  qgram.frequencies.value[0], qgram.frequencies.value[-1]],
          aspect='auto', origin='lower', cmap='viridis', interpolation='bilinear')
ax3.plot(t_true, f_true, 'r--', linewidth=1.5, alpha=0.8)
ax3.set_ylim(10, 150)
ax3.set_ylabel('Frequency (Hz)')
ax3.set_title('Q-transform')

# HHT instantaneous frequency
ax4 = plt.subplot(4, 1, 4, sharex=ax1)
ax4.plot(inst_freq.times.value, inst_freq.value, linewidth=1, label='HHT Inst. Freq')
ax4.plot(t_true, f_true, 'r--', linewidth=2, label='True Frequency')
ax4.set_ylim(0, 150)
ax4.set_ylabel('Frequency (Hz)')
ax4.set_xlabel('Time (s)')
ax4.set_title('HHT Instantaneous Frequency')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Decision Tree

What information do you need?

1. General time-frequency view?
   β†’ Use **Spectrogram (STFT)**

2. Transient detection with adaptive resolution?
   β†’ Use **Q-transform**

3. Chirp frequency tracking (multi-scale)?
   β†’ Use **Wavelet (CWT)**

4. Instantaneous frequency as single-valued curve?
   β†’ Use **HHT**

5. Decay rates AND frequencies (ringdown)?
   β†’ Use **STLT**

6. Echo detection or spectral periodicity?
   β†’ Use **Cepstrum**

7. Signal compression or smooth feature extraction?
   β†’ Use **DCT**

Quick decision path:
β”œβ”€ Stationary signal? β†’ Spectrogram
β”œβ”€ Transient/burst? β†’ Q-transform
β”œβ”€ Chirp spanning octaves? β†’ Wavelet
β”œβ”€ Need inst. frequency? β†’ HHT
β”œβ”€ Damped oscillations? β†’ STLT
β”œβ”€ Has echoes? β†’ Cepstrum
└─ Smooth, need compression? β†’ DCT

Performance Comparison

Method

Computation Time*

Memory

Best For

Spectrogram (STFT)

1Γ— (baseline)

Low

General purpose

DCT

1-2Γ—

Low

Compression

Q-transform

5-10Γ—

Medium

Transients

Wavelet (CWT)

8-15Γ—

Medium-High

Chirps

Cepstrum

3-5Γ—

Medium

Echo detection

STLT

15-25Γ—

High

Ringdown modes

HHT

20-50Γ—

High

Instantaneous freq.

*Approximate, depends on parameters and signal length

Summary Table

Feature

STFT

Q-transform

Wavelet

HHT

STLT

Cepstrum

DCT

Resolution

Fixed

Adaptive Q

Scale-adaptive

Data-adaptive

Fixed

Quefrency

N/A

Best signal

Stationary

Chirps

Multi-scale

AM/FM

Damped

Echoes

Smooth

Comp. cost

Low

Medium

Med-High

Very High

High

Medium

Low

Output

2D (t,f)

2D (t,f)

2D (t,f)

f(t) curve

3D (t,Οƒ,Ο‰)

Ο„ spectrum

Coeffs

Time-localized

Yes

Yes

Yes

Yes

Yes

No

No

Unique info

β€”

Adaptive res.

Multi-scale

Inst. freq.

Decay rate Οƒ

Delays

Compression

GW use

General

Transients

Chirps

Glitches

Ringdown

Reflections

Features

Practical Recommendations

For Routine Analysis

Start with Spectrogram (STFT) - it’s fast, well-understood, and sufficient for most cases.

For Transient Detection

Use Q-transform - it’s the standard for gravitational wave burst searches.

For Chirp Analysis

Use Wavelet (CWT) for precise frequency trajectory tracking across multiple scales.

For Instantaneous Frequency

Use HHT when you need single-valued instantaneous frequency (no time-frequency uncertainty).

For Ringdown Analysis

Use STLT to estimate both frequencies and decay rates (quality factors) simultaneously.

For Echo Detection

Use Cepstrum to identify delay times and periodic structure in the spectrum.

For Data Reduction

Use DCT for efficient compression and feature extraction of smooth signals.

For Publication

  • Include Spectrogram (familiar baseline for readers)

  • Add specialized method(s) that reveal unique physics (Wavelet for chirps, HHT for inst. freq., STLT for ringdown, etc.)

  • Show robustness: demonstrate that key results hold across multiple methods

General Strategy

  1. Start simple: Always begin with Spectrogram

  2. Identify limitations: Where does STFT fail to show important features?

  3. Choose specialized method: Pick the method whose β€œunique info” matches your needs (see Summary Table)

  4. Validate: Cross-check with another method if possible

See Also: