The Spectral Estimation Problem: Beyond Band Power

From band power to a random process

In Vol I you computed band power by filtering, squaring, and averaging — a perfectly reasonable recipe for a demo. Volume II starts from an uncomfortable fact: the power spectral density (PSD) you are after is a property of a stochastic process, not of your particular two-second segment. Your data are one finite, noisy realization; the PSD is a population parameter you must estimate. This whole track is about doing that estimation well.

S(f) = \sum_{\tau=-\infty}^{\infty} r(\tau)\, e^{-i2\pi f\tau}, \qquad r(\tau)=\mathbb{E}\!\left[x_n\,x_{n+\tau}\right]

Wiener–Khinchin: the PSD is the Fourier transform of the autocovariance.

The power spectrum — how much signal power sits at each frequency — is just the Fourier transform of the autocovariance, which measures how the signal correlates with a time-shifted copy of itself. Rhythm in that correlation shows up as peaks in the spectrum.

S(f): The power spectral density at frequency f.
r(\tau): Autocovariance — correlation of the signal with itself at lag \tau.
e^{-i2\pi f\tau}: The Fourier kernel, turning lags into frequencies.
\mathbb{E}: Expectation — an average over the random process.

The periodogram: seductive but biased

The obvious estimator is the periodogram — the squared magnitude of the DFT. For white noise it is asymptotically unbiased, but for colored spectra at finite N it is biased, because the finite window convolves the true spectrum with the window's spectral kernel. That smearing is spectral leakage: power from strong components (line noise, a large low-frequency rhythm) bleeds into neighbouring bins and can bury weak high-frequency activity.

\hat{S}_p(f) = \frac{1}{N}\left|\sum_{n=0}^{N-1} x_n\, e^{-i2\pi f n}\right|^2

The periodogram: squared magnitude of the DFT.

The most obvious spectral estimate: take the Fourier transform of your data and square it. It's easy and tempting, but as the next blocks show, it's statistically flawed — noisy no matter how much data you collect.

\hat{S}_p(f): The periodogram estimate of the spectrum.
x_n: The data samples.
\sum_n x_n e^{-i2\pi f n}: The discrete Fourier transform (DFT) of the data.
|\cdot|^{2}: Squared magnitude — i.e. power.

The expected periodogram is the true spectrum convolved with the Fejér kernel (the squared Dirichlet kernel) of the rectangular window. Its main lobe sets the frequency resolution (\approx 1/T); its side lobes, falling off only as \sim 1/f^2 (first side lobe near −13 dB), are the leakage. Tapering (Hann, Hamming, …) trades a wider main lobe for far lower side lobes.

Inconsistency: variance that never shrinks

Leakage is only half the trouble. The periodogram is inconsistent: its variance does not go to zero as you collect more data. Adding samples buys finer frequency resolution, not a smoother estimate — each new frequency bin is essentially one more noisy sample.

\operatorname{Var}\!\big[\hat{S}_p(f)\big] \;\approx\; S^2(f) \qquad \text{(independent of } N \text{)}

For a Gaussian process the periodogram's variance is ~S², regardless of record length.

The nasty surprise: the periodogram's error doesn't shrink as you gather more data. Its variance stays about the size of the spectrum itself, so a raw periodogram is always jagged. That's exactly why we need averaging methods like Welch and multitaper.

\operatorname{Var}[\hat{S}_p(f)]: The variance (uncertainty) of the periodogram estimate.
S^{2}(f): The squared true spectrum — the floor the variance never drops below.
N: The record length; the variance is independent of it.

Stationarity and the bias–variance dial

Averaging independent looks assumes the statistics don't change across those looks — stationarity. Brain signals are at best locally stationary: an event-related desynchronization can drop mu power within a few hundred milliseconds. So estimator design is really a three-way tension between frequency resolution, variance (statistical stability), and time resolution. You can have any two at the expense of the third.

Choose a segment short enough that the signal is ~stationary within it, but long enough for the frequency resolution you need (\Delta f \approx 1/T).
Choose how to taper each segment to control leakage.
Choose how many independent looks to average (segments, tapers, or trials) to hit your target variance.

Map of the track

With that framing, the rest of the track is a tour of principled ways to spend the bias–variance–resolution budget: Guide 2 builds the workhorse nonparametric estimators (Welch and multitaper); Guide 3 moves to time-resolved analysis (STFT, wavelets, the analytic signal) and the uncertainty floor; Guide 4 couples channels and frequencies (coherence and phase–amplitude coupling); Guide 5 confronts the aperiodic background and turns spectra into a working BCI via frequency tagging.