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, misleading, and no greener for more data.

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 independence lost to windowing, letting you average more segments from the same data. You are, once again, spending resolution to buy stability: the exchange rate you pick depends on whether you need to resolve two close resonances or measure 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. The cross-spectrum also carries phase, \(\angle S_{xy}(f)\), and for a pure delay \(y(t) = x(t-\tau)\) that phase grows linearly with frequency, \(\angle S_{xy}(f) = -2\pi f \tau\), so dividing the unwrapped phase by \(-2\pi f\) recovers \(\tau\) directly: a phase slope, not a single phase value, is the quantity you fit. That slope is trustworthy only where coherence is high, because phase is meaningless once noise dominates a bin (it becomes an average of arbitrary angles), and because a phase that changes quickly with frequency can wrap past \(\pm\pi\) and alias into a shorter, wrong apparent delay unless you unwrap across frequency first. Coherence-weighted, unwrapped phase 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 for a roughly stationary signal when you want a stable, reproducible view of power across frequency: a sensor's noise floor (Chapter 2), a vibration environment, or band-power features for the classifiers of Section 7.7 and Chapter 8. Use coherence when the question is coupling between two channels: EEG functional connectivity (Chapter 31), structural transfer functions, leak localization, or sensor-fusion checks that redundant sensors actually agree (Chapter 48). Skip both when the signal is strongly nonstationary within your averaging window, since averaging then blurs distinct regimes into a meaningless mean; 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.