Spectral Analysis

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

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

Common Pitfalls