Spectral Analysis

Overview

Extract frequency content from physics signals using FFT, power spectral density estimation, spectrograms, and wavelet transforms. Covers stationary and non-stationary signals.

When to Use

• Finding oscillation frequencies in time-series data
• Computing power spectra of turbulence, vibrations, or wave data
• Detecting transient features (spectrograms, wavelets)
• Filtering signals (low-pass, band-pass, notch)
• Cross-spectral analysis and coherence between signals

Core Workflows

1. FFT (Fast Fourier Transform)

import numpy as np
import matplotlib.pyplot as plt

# Signal: sum of two sinusoids + noise
dt = 0.001  # sampling interval [s]
fs = 1 / dt  # sampling frequency [Hz]
t = np.arange(0, 1, dt)
signal = 1.5 * np.sin(2*np.pi*50*t) + 0.8 * np.sin(2*np.pi*120*t)
signal += 0.5 * np.random.randn(len(t))

# Compute FFT
N = len(t)
fft_vals = np.fft.rfft(signal)
freqs = np.fft.rfftfreq(N, d=dt)
amplitude = 2.0 / N * np.abs(fft_vals)  # single-sided amplitude
phase = np.angle(fft_vals)

fig, axes = plt.subplots(2, 1, figsize=(10, 7))
axes[0].plot(t[:200], signal[:200], 'b-', linewidth=0.8)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Time Domain Signal')
axes[0].grid(True, alpha=0.3)

axes[1].plot(freqs, amplitude, 'r-', linewidth=0.8)
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('FFT Amplitude Spectrum')
axes[1].set_xlim(0, 200)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fft_analysis.png', dpi=150, bbox_inches='tight')

2. Power Spectral Density (Welch Method)

from scipy.signal import welch, periodogram

# Welch PSD (better noise averaging than raw FFT)
freqs_w, psd_w = welch(signal, fs=fs, nperseg=256, noverlap=128)

# Periodogram (raw, no averaging)
freqs_p, psd_p = periodogram(signal, fs=fs)

fig, ax = plt.subplots(figsize=(10, 5))
ax.semilogy(freqs_p, psd_p, 'gray', alpha=0.3, label='Periodogram')
ax.semilogy(freqs_w, psd_w, 'r-', linewidth=2, label='Welch PSD')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('PSD [V²/Hz]')
ax.set_title('Power Spectral Density')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 200)
plt.savefig('psd.png', dpi=150, bbox_inches='tight')

Welch parameters:

• nperseg: Segment length (frequency resolution = fs/nperseg)
• noverlap: Overlap between segments (typically 50% = nperseg//2)
• window: Window function ('hann' default, 'blackman' for better sidelobe suppression)
• Trade-off: longer segments → better frequency resolution, worse noise averaging

3. Spectrogram (Time-Frequency)

from scipy.signal import spectrogram

# Chirp signal (frequency sweeps from 10 to 200 Hz)
t_chirp = np.linspace(0, 2, 8000)
chirp = np.sin(2*np.pi * (10*t_chirp + 47.5*t_chirp**2))

f_spec, t_spec, Sxx = spectrogram(chirp, fs=4000, nperseg=256, noverlap=200)

fig, ax = plt.subplots(figsize=(10, 5))
pcm = ax.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx + 1e-20),
                     shading='gouraud', cmap='inferno')
ax.set_ylabel('Frequency [Hz]')
ax.set_xlabel('Time [s]')
ax.set_title('Spectrogram')
plt.colorbar(pcm, ax=ax, label='PSD [dB/Hz]')
plt.savefig('spectrogram.png', dpi=150, bbox_inches='tight')

4. Wavelet Transform (Continuous)

from scipy.signal import cwt, morlet2

# Continuous wavelet transform with Morlet wavelet
widths = np.geomspace(1, 128, num=100)  # scale parameters
cwtmatr = cwt(signal, morlet2, widths, w=6)

# Convert scales to frequencies: f = w * fs / (2*pi*scale)
frequencies = 6 * fs / (2 * np.pi * widths)

fig, ax = plt.subplots(figsize=(10, 5))
ax.pcolormesh(t, frequencies, np.abs(cwtmatr),
              shading='gouraud', cmap='viridis')
ax.set_ylabel('Frequency [Hz]')
ax.set_xlabel('Time [s]')
ax.set_title('Continuous Wavelet Transform (Morlet)')
ax.set_ylim(0, 200)
plt.colorbar(ax.collections[0], ax=ax, label='|CWT|')
plt.savefig('cwt.png', dpi=150, bbox_inches='tight')

5. Cross-Spectral Analysis and Coherence

from scipy.signal import coherence, csd

# Two related signals
signal2 = 1.2 * np.sin(2*np.pi*50*t + 0.3) + 0.6 * np.random.randn(len(t))

# Coherence (how correlated are the two signals at each frequency)
f_coh, Cxy = coherence(signal, signal2, fs=fs, nperseg=256)

# Cross-spectral density
f_csd, Pxy = csd(signal, signal2, fs=fs, nperseg=256)

fig, axes = plt.subplots(2, 1, figsize=(10, 7))
axes[0].plot(f_coh, Cxy, 'b-', linewidth=1.5)
axes[0].set_ylabel('Coherence')
axes[0].set_title('Coherence between signals')
axes[0].set_xlim(0, 200)
axes[0].grid(True, alpha=0.3)

axes[1].semilogy(f_csd, np.abs(Pxy), 'r-', linewidth=1.5)
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('|CSD| [V²/Hz]')
axes[1].set_title('Cross-Spectral Density')
axes[1].set_xlim(0, 200)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('coherence.png', dpi=150, bbox_inches='tight')

6. Filtering

from scipy.signal import butter, sosfilt, sosfiltfilt

def bandpass_filter(data, lowcut, highcut, fs, order=4):
    """Apply zero-phase Butterworth bandpass filter."""
    nyq = 0.5 * fs
    sos = butter(order, [lowcut/nyq, highcut/nyq], btype='band', output='sos')
    return sosfiltfilt(sos, data)

# Example: extract 45-55 Hz component
filtered = bandpass_filter(signal, 45, 55, fs, order=4)

Frequency Resolution vs Averaging Trade-off

Parameter	Effect on Resolution	Effect on Noise
Longer nperseg	Better frequency resolution	More noise (fewer averages)
Shorter nperseg	Worse frequency resolution	Less noise (more averages)
More noverlap	Same resolution	Slightly less noise
Windowing (Hann)	Wider main lobe	Better sidelobe suppression

Rule of thumb: Frequency resolution Δf = fs / nperseg

Common Pitfalls

Pitfall	Fix
Aliasing (fs too low)	Nyquist: fs ≥ 2 × f_max
Spectral leakage	Apply window function (Hann, Blackman)
Zero-padding confusion	Zero-padding interpolates FFT, doesn't improve resolution
Wrong units on PSD	Check: V²/Hz for continuous, V² for discrete
DC component dominates	Subtract mean before FFT
Non-uniform sampling	Use Lomb-Scargle periodogram (scipy.signal.lombscargle)