Power spectral density and coherence

"A single FFT tells you what one window happened to contain. A power spectral density tells you what the process keeps doing."

A Patient AI Agent

This section assumes the sampling, aliasing, and windowing ideas from Chapter 3 and the leakage-versus-resolution tradeoff of Section 7.1. It also leans on the notion of a stationary random process and its autocorrelation from the probability primer in Chapter 4.

What a PSD is, and why the periodogram alone fails

For a wide-sense-stationary process \(x(t)\) with autocorrelation \(r_x(\tau) = \mathbb{E}[x(t)x(t+\tau)]\), the power spectral density is its Fourier transform (the Wiener-Khinchin theorem):

The units are power per hertz: \((\text{units of }x)^2/\text{Hz}\), so an accelerometer PSD reads in \(g^2/\text{Hz}\) and integrating \(S_x(f)\) over a band returns the mean-square (the power) the signal carries in that band. That integral property is what makes a PSD a measurement rather than a picture: the area under a resonance peak is the energy that will fatigue a bracket or heat a coil.

The naive estimator is the periodogram, \(\hat{S}(f) = \tfrac{1}{N}\lvert X(f)\rvert^2\), the squared magnitude of the DFT. It is asymptotically unbiased, but it has a fatal defect: its variance does not decrease as the record length \(N\) grows. Doubling your data doubles the number of frequency bins but leaves each bin just as noisy, roughly 100% relative error, because each bin is essentially a single sample of a chi-squared random variable. A periodogram of a smooth process looks like grass: spiky, irreproducible, and misleading.

Welch's method: trading resolution for a variance you can control

The cure is averaging, and the standard recipe is Welch's method. Split the record into \(K\) segments, overlapping (typically 50%), apply a window \(w[n]\) to each to control leakage, compute a periodogram per segment, and average them:

where \(U = \sum_n w[n]^2\) normalizes for the window's power so the PSD stays correctly scaled. Averaging \(K\) roughly independent periodograms cuts the variance by about \(K\). The catch is the classic bias-variance dial: shorter segments give more of them (lower variance) but coarser frequency resolution \(\Delta f = f_s / L\) for segment length \(L\), and a fatter window main lobe (more bias). Overlap recovers some independence lost to windowing, letting you average more segments from the same data. You are, once again, spending resolution to buy stability, and the right operating point depends on whether you are resolving two close resonances or estimating a broadband floor.

Cross-spectra and coherence: coupling between two channels

Sensor systems rarely have one channel. Given two stationary signals \(x\) and \(y\), the cross-spectral density \(S_{xy}(f)\) is the Fourier transform of their cross-correlation, and estimated exactly like a Welch PSD but multiplying one segment's DFT by the complex conjugate of the other's before averaging. From it comes the magnitude-squared coherence:

Coherence is the frequency-domain analogue of a squared correlation coefficient: at each frequency it reports how much of \(y\)'s power is linearly predictable from \(x\). A value near 1 means the two channels are tightly, linearly coupled at that frequency (a common source, or one driving the other through a linear path); a value near 0 means they are unrelated there. Its phase, \(\angle S_{xy}(f)\), gives the frequency-dependent time delay of the coupling, which is how you locate a leak by the lag between two acoustic sensors or estimate a transfer function \(H(f) = S_{xy}(f)/S_x(f)\) for a structure under vibration.

Computing it in practice

The code below estimates a PSD, a cross-spectrum, and coherence for two coupled signals: \(y\) is a filtered, delayed copy of a shared band-limited source plus independent noise. Welch handles the segmentation, windowing, overlap, and averaging; coherence peaks in the shared band and stays near the \(1/K\) floor elsewhere.

import numpy as np
from scipy import signal

fs = 8192
n = 200_000
rng = np.random.default_rng(0)
src = signal.lfilter(*signal.butter(4, [30, 45], "bandpass", fs=fs),
                     x=rng.standard_normal(n))          # shared 30-45 Hz source
x = src + 0.5 * rng.standard_normal(n)                  # channel 1 + noise
y = np.roll(src, 12) + 0.5 * rng.standard_normal(n)     # delayed copy + noise

f, Sxx = signal.welch(x, fs, nperseg=4096)              # PSD, units^2/Hz
f, Sxy = signal.csd(x, y, fs, nperseg=4096)             # cross-spectral density
f, Cxy = signal.coherence(x, y, fs, nperseg=4096)       # magnitude-sq coherence

band = (f > 30) & (f < 45)
print(f"mean coherence in 30-45 Hz band: {Cxy[band].mean():.2f}")
print(f"mean coherence outside band:     {Cxy[~band].mean():.2f}")

When to reach for PSD and coherence

Use a PSD whenever you care about how power distributes across frequency for a roughly stationary signal and you want a stable, reproducible estimate: characterizing a sensor's noise floor (as in Chapter 2), specifying a vibration environment, or building band-power features for the classifiers of Section 7.7 and Chapter 8. Use coherence whenever the question is about coupling between two channels: EEG functional connectivity between electrodes (Chapter 31), structural transfer functions, leak localization, or sensor-fusion sanity checks that two supposedly redundant sensors actually agree (Chapter 48). Do not use either when the signal is strongly nonstationary within your averaging window: averaging then blurs distinct regimes into a meaningless mean, and you should segment first or move to the time-frequency tools of Section 7.2. And remember coherence measures only linear coupling; two channels can be perfectly dependent through a nonlinear law yet show low coherence.

What's Next

In Section 7.7, we step back from individual transforms and ask the design question that governs this whole chapter: given a modality and a task, which representation earns its keep. We will weigh the FFT, spectrogram, wavelet, Hilbert-Huang, cepstral, and PSD tools you have now met against the signal's stationarity, the compute budget, and what a downstream model can actually use, turning a menu of methods into a decision.