The periodogram

The obvious way to get a spectrum from a recording is to square its DFT and divide by the number of samples. Here that is tried on white noise, whose true spectrum is flat at 0 dB, with 64, 256 and then 1024 samples. Watch the dots: they get denser, but no tighter.

A spectrum from one record

Power spectral density (24.4) gave a random process a spectrum, the PSD $S_x(e^{j\Omega})$. In “Memory in time, narrowness in frequency” it came from a formula, because the process was known. In practice you have a recording, not a formula, and you want its spectrum anyway.

A recording is one realisation, in the words of Random processes (24.2). As in “Stationary, and ergodic”, I assume the process is ergodic, so that one long recording can stand for the whole ensemble. This page asks how well it does.

Take the $N$ samples $x[0]$ to $x[N-1]$ and their DTFT, the sum of turned samples of “Samples, turned and added” in The DTFT (12.2):

Square its size and divide by $N$. That is the periodogram of the $N$ samples:

Arthur Schuster gave it that name in 1898, when he used it to hunt for hidden periods in meteorological records. The letter $P$ is the average power of How big is a signal (1.3); the subscript $N$ and the argument mark the periodogram.

Why divide by $N$? At the bins, $\Omega=2\pi k/N$, the DTFT is the DFT $X[k]$, as Sampling the spectrum (13.1) showed. By “Energy in time and in bins” in The DFT (13.2), $\sum_n x[n]^2=\frac1N\sum_k\lvert X[k]\rvert^2$. Divide both sides by $N$: the $N$ values of $\lvert X[k]\rvert^2/N$ average to the mean square of the samples, their power.

So the periodogram spreads the record’s power over frequency, as a PSD does. Its area works out too: $\frac1{2\pi}\int_{-\pi}^{\pi}P_N\,d\Omega$ is the same mean square, just as 24.4’s PSD integrates to $R_x[0]$.

To read it in volts squared per hertz, one-sided, as in “Why the floor depends on N” of Reading a spectrum: scaling and units (15.4), multiply it by $2/f_s$. Bins 0 and $N/2$ are not doubled. That is 15.4’s PSD formula with no window, where every $w[n]$ is 1 and $\sum_nw[n]^2=N$.

Bias, variance and consistency

The periodogram is a guess, and a guess needs words. An estimator turns data into a guess at something fixed. I write any estimate of the PSD as $\hat S_x(e^{j\Omega})$, with the hat of Random variables for signals (24.1). The periodogram is one such estimate.

A new record gives a new guess, so an estimate is itself a random variable, with a mean and a spread. Its bias is how far its average lies from the truth. Its variance is how widely it scatters about its own average. It is consistent if both shrink to 0 as $N$ grows: then a long enough record pins the truth down.

Think of a bathroom scale. One that always reads 2 kg heavy is biased. One whose needle wobbles has variance. A good scale has neither, and standing on it longer should only help.

I measure scatter by the spread, the standard deviation of “Draws pile up into a density” in 24.1. Dividing it by the mean gives the readout “spread ÷ mean”: 0.1 means the guesses wander by a tenth of their size.

Let’s test the periodogram where the truth is known. White Gaussian noise of variance 1 has a flat PSD (24.4): $S_x=1$, which is 0 dB as a power ($10\log_{10}$, 1.3), at every $\Omega$. Its periodogram averages exactly to the truth, as the next section shows. So whatever you see is scatter, and you would hope it shrinks as $N$ grows.

More data, same scatter

The picture at the top of the page draws the periodogram of white noise for the first 64, 256 and 1024 samples of one record, one dot per bin. Its shaded band holds 90 % of values; the next section says where the band’s edges come from. Spread ÷ mean reads 1.099, 0.944 and 0.963 for the three lengths.

Its readouts are taken over bins 1 to N/2 − 1, strictly between 0 and π. For a real signal the bins above $N/2$ mirror those below (“Bins in hertz, and the mirror”, 13.2), and bins 0 and $N/2$ are real numbers that follow a different law.

Notice that the dots get denser but no tighter. Spread ÷ mean stays near 1 at every $N$: each dot is wrong by about the size of the truth itself. It is like asking more people a question that each of them answers at random. You get more answers, not better ones.

The mean readout stays near 1. It sits close to each record’s own mean square, 1.128 for the first 64 samples and 1.026 for all 1024, because by 13.2’s energy rule all $N$ values average to exactly that.

Why the scatter stays

In 12.2 each sample became an arrow, turned back by $\Omega n$, and the arrows were added nose to tail. With noise, the arrows have random lengths and signs. Their sum is the end point of a random walk of $N$ steps.

How far does it get on average? Multiply out $\lvert X\rvert^2=X\,X^*$ and take the average term by term, since averages add (24.2). Each product $x[n]\,x[m]$ averages to $R_x[n-m]$, the autocorrelation of 24.2. For white noise that is $\sigma_x^2$ when $m=n$ and 0 otherwise:

So the squared length grows like $N$ on average, and dividing by $N$ makes the average right: $\mathbb{E}\{P_N\}=\sigma_x^2$ at every $\Omega$ and every $N$. For white noise the periodogram has no bias.

But one walk’s end is random. Sometimes the steps nearly cancel and sometimes they line up. The squared length scatters about as widely as its own average, and both grow like $N$. Dividing by $N$ scales them together, so spread ÷ mean stays near 1 however long the record is.

For Gaussian noise the law of that scatter is known. I assume it here; statistics derives it. At every bin strictly between 0 and $\pi$, $P_N$ follows the exponential law: the share of values above $c$ times the truth is $e^{-c}$,

Its spread equals its mean, so the law’s spread ÷ mean is 1. The measured 1.099, 0.944 and 0.963 scatter about it, because 31 to 511 dots only estimate the law.

Now read the band from it. A share 0.05 lies above $c=-\ln0.05=2.996$ times the truth, which is 4.8 dB. A share 0.05 lies below $c=-\ln0.95=0.0513$ times the truth, where $1-e^{-c}=0.05$, which is $-12.9$ dB. So 90 % of values lie between $-12.9$ dB and 4.8 dB of the truth, at any $N$, and the 511 dots put 90.4 % there.

The band is lopsided: it reaches 12.9 dB down but only 4.8 dB up. On a dB axis the low tail is long, and more dots reach further into it: the lowest is $-16.5$ dB at $N=64$ and $-35.2$ dB at 1024.

Even the middle sits low. Half of all values lie below $\ln2=0.693$ of the truth, $-1.6$ dB, and 63.2 % below the truth itself, since $1-e^{-1}=0.632$. Of the 511 dots, 313 do, 61.3 %.

So the periodogram fails the last test. Its bias is 0 for white noise, but its spread stays at the size of $S_x$ itself whatever $N$ is, so its variance never shrinks: it is not consistent. A longer record buys more bins, closer together, each as unreliable as before.

Short records blur the peak

White noise is the easy case, because its PSD is flat. A PSD with a sharp peak shows the second fault. First I need the periodogram’s average for any process, not only white noise.

Multiply out $\lvert X\rvert^2$ as before, but now group the products $x[n]\,x[m]$ by their lag $\ell=n-m$. There are $N-\lvert\ell\rvert$ of them at lag $\ell$.

Added up, they are $N-\lvert\ell\rvert$ times 24.2’s estimate $\hat R_x[\ell]$ from “White noise forgets, coloured noise remembers”, which divides by the number of products. A negative lag uses the same products as the positive one, so $\hat R_x[-\ell]=\hat R_x[\ell]$. Dividing by $N$,

The periodogram is the DTFT of the estimated autocorrelation, with its lags weighted by a triangle.

Each product at lag $\ell$ averages to the true $R_x[\ell]$, so $\hat R_x[\ell]$ does too, and averages add:

Without the triangle this sum is the DTFT of $R_x[\ell]$, which 24.4 showed is the true PSD. With it the far lags are cut down, so the average is the truth seen through a window. For white noise only $\ell=0$ is left, and the triangle is 1 there: no bias, as found above.

Which window? Put the record window itself into the lag form: $N$ ones, the cut of “Cut a tone short, and its arrows spread” in Properties of the DTFT (12.3). Every product is 1, so every $\hat R$ is 1, and the right side is the DTFT of the triangle alone. The left side is $\lvert W(e^{j\Omega})\rvert^2/N$, with 12.3’s $W$:

A product of two sequences in the lag domain is a periodic convolution of their transforms, by “Multiplication in time is periodic convolution in frequency” (12.3). So the average periodogram is the true PSD smoothed by this kernel:

The kernel is never negative, and its own area, $\frac1{2\pi}\int$, is 1, the triangle’s value at $\ell=0$. So the average at each $\Omega$ is a weighted average of the true PSD nearby. An average can never exceed the largest value it averages, or fall below the smallest, so a peak can only drop and a valley can only fill.

How wide is the smoothing? $\lvert W\rvert^2$ has the zeros of $W$, so its main lobe is 12.3’s $4\pi/N$ wide. Its side lobes are the leakage of “A tone between bins leaks into all of them” in Windowing & spectral leakage (15.1).

A noise with a sharp peak

The test process is the ringing recursion of “Where the pair sits is how h[n] moves” in Transfer functions, poles & zeros (16.3), driven by white noise $v[n]$ of variance $\sigma_v^2=1$:

It is 16.3’s pair with radius $r=0.95$ at the angle $0.3\pi$: $2\cdot0.95\cos(0.3\pi)=1.1168$ and $0.95^2=0.9025$, so its poles are $0.95\angle\pm0.3\pi$. A sample depends on its last two, so this is an AR(2) process, 24.2’s AR(1) with one more past value.

Its PSD follows from $S_y=\lvert H\rvert^2S_x$ of 24.4, with white noise going in. Here $H(e^{j\Omega})=1/(1-1.1168\,e^{-j\Omega}+0.9025\,e^{-j2\Omega})$, 16.3’s transfer function on the unit circle, so

It peaks at 22.1 dB at $\Omega=0.2997\pi$, just below the poles’ angle, and falls to $-9.6$ dB at $\pi$. The peak is narrow. The rule of “Poles behind the zeros narrow the notch” in Resonators, notches and combs (17.4), a width of about $2(1-r)=0.1$ rad/sample, holds well, and the $-3$ dB points are 0.103 rad/sample apart. Its power is $R_x[0]=8.23$, and 90 % of it lies between $0.2\pi$ and $0.4\pi$ on either side of 0.

Think of a photograph taken through a small aperture. The bright spot spreads out, and its light reaches the dark corners.

Watch the solid curve close in on the dashed one as $N$ grows. It is the exact average of the sum above, not one noisy record.

Notice the first frame: at 32 samples the valley fills by 3.8 dB, more than the peak drops, 2.9 dB.

Both follow from the kernel’s width. At $N=32$ its main lobe is $4\pi/32=0.39$ rad/sample wide, nearly four times the peak’s 0.103, so the peak is smeared over its neighbours. At 128 the main lobe, 0.098 rad/sample, is about as wide as the peak, and the peak falls only 0.7 dB short. At 1024 it is 0.012 rad/sample, and the peak is within 0.1 dB.

The valley at $\pi$ sits 31.7 dB below the peak, so even the faint side lobes of the kernel carry enough of the peak’s power there to matter. At $N=32$, 52 % of the average at $\pi$ comes from the band around the peak, $0.2\pi$ to $0.4\pi$ on either side of 0. That is leakage, the same as in 15.1, only now it leaks noise power instead of a tone.

One root, two cures

Both faults come from the lag form. One record gives one estimate $\hat R_x[\ell]$ per lag, and the far lags rest on few products, down to a single one at $\ell=N-1$. Every one of them enters the sum, noise and all: that is the variance. The triangle tapers the far lags away: that is the bias.

The cures follow. A smoother window than the plain cut has lower side lobes, as in Window functions compared (15.2), and so leaks less. Averaging the periodograms of many pieces of the record cuts the variance; that is the subject of Averaged periodograms: Bartlett and Welch (25.2).

Worked example

1. A periodogram by hand. Take $x$ = 1, 2, 1, 0, so $N=4$. The DFT is $X[k]$ = 4, $-2j$, 0, $2j$, and the periodogram at the bins is $\lvert X[k]\rvert^2/4$ = 4, 1, 0, 1. These average to 1.5, the mean square $(1+4+1+0)/4$.

Now by lags. 24.2’s estimates are $\hat R_x[0]=6/4=1.5$, $\hat R_x[1]=(2+2+0)/3=4/3$, $\hat R_x[2]=(1+0)/2=1/2$ and $\hat R_x[3]=0$. With the triangle’s weights 1, 3/4, 2/4 and 1/4, the lag form gives

At $\Omega$ = 0, $\pi/2$ and $\pi$ this is 4, 1 and 0, as before. Between bins it gives new values: at $\pi/3$ it is 2.25.

Scale it as a one-sided PSD for $f_s=1000$ Hz, with $x$ in volts. The bins are 0, 250 and 500 Hz. Bin 0 is $4/1000=0.004$ V²/Hz, bin 1 is doubled, $2\times1/1000=0.002$ V²/Hz, and bin 2 is 0. Times the bin width, 250 Hz, they add to $(0.004+0.002)\times250=1.5$ V², the mean square again.

SciPy’s periodogram(x, fs=1000, detrend=False) returns exactly these numbers. By default it first subtracts the mean, and then bin 0 reads 0.

2. White noise. For $N$ = 64, 256 and 1024 seeded samples, spread ÷ mean is 1.099, 0.944 and 0.963. The exponential law says 1 at every $N$, and puts 90 % of the values between $-12.9$ dB and 4.8 dB of the truth.

3. AR(2) bias. At 32 samples the average periodogram’s peak is 19.1 dB against a true 22.1 dB, 2.9 dB short. At $\pi$ it is $-5.8$ dB against a true $-9.6$ dB, 3.8 dB high. At 1024 samples both are within 0.3 dB.

Where you’ll meet this

The periodogram is the raw material of almost every spectral estimate. Averaged periodograms: Bartlett and Welch (25.2) averages the periodograms of segments. Each column of a spectrogram in Spectrograms & the STFT (15.5) is, up to a scale factor, the periodogram of one windowed frame.

Alone, it serves best for strong tones in weak noise. A tone of amplitude $A$ on a bin gives $\lvert X[k]\rvert=AN/2$ (13.2), so its periodogram value $NA^2/4$ grows with $N$, while the noise bins only scatter about their truth. At $N=1024$ a tone of amplitude 1 in white noise of variance 1 reads 24.1 dB, 19.3 dB above the top of the 90 % band.

For more, see A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, ch. 10; J. G. Proakis and D. G. Manolakis, Digital Signal Processing, ch. 14; and P. Stoica and R. Moses, Spectral Analysis of Signals (2005), ch. 2. In NumPy, np.abs(np.fft.fft(x))**2 / len(x) gives $P_N$ at the bins, and SciPy’s scipy.signal.periodogram gives the one-sided density. Models that sidestep the periodogram’s faults come in Parametric models and linear prediction (25.3).

Reference card