Averaged periodograms: Bartlett and Welch

Cut one long record of noise into shorter segments, and average their periodograms. Here 4096 samples are cut into 1, 4, 16 and then 64 segments. Watch the estimate smooth out, while its peak sinks below the dashed truth.

Average the pieces: less scatter, less detail

The periodogram (25.1) left us with two faults. In “More data, same scatter” its values scattered about the true PSD by as much as the PSD itself, however long the record was. In “Short records blur the peak” a short record smeared a sharp peak and lifted the valleys.

This page cures the first fault with an average, and pays for it with the second. Here is the recipe.

Take a record of $N$ samples and cut it into $K$ pieces, called segments, each $L$ samples long, so $K=N/L$. Earlier chapters used $K$ for a harmonic and for a gain, and $L$ for an interpolation factor; on this page $K$ counts segments and $L$ is a length.

Take each segment’s periodogram, add the $K$ of them, and divide by $K$. That average is the estimate $\hat S_x(e^{j\Omega})$; the hat means “estimated”, as in 25.1. Segment $i$ is the stretch $x[iL]$ to $x[iL+L-1]$, counting from $i=0$, and its periodogram is 25.1’s with $L$ samples:

This average is Bartlett’s method. Why should it help? One segment’s periodogram scatters about its own average by about that average (25.1).

Suppose the segments are nearly independent of each other. Then “A longer average: less noise, more delay” in Simple smoothing filters (18.2) gives the answer: averaging $K$ unrelated values divides their spread by $\sqrt K$. Four segments halve the scatter, and sixteen leave a quarter.

I assume this on the whole page. It is fair here: for this page’s noise, samples more than 57 apart have a correlation coefficient under 0.05, and every segment is at least 64 long.

The price is detail. Each periodogram sees only $L$ samples, so it blurs like 25.1’s short record. In “Five windows, one trade” of Window functions compared (15.2), the rectangular window’s main lobe is 2 bins wide, and one bin is $2\pi/L$. So the estimate cannot show details much narrower than about $2\pi/L$.

It is like several quick glances at a scene against one long look: steadier, but less detailed. With $N$ fixed, more segments means shorter ones, and you cannot have both.

To measure the scatter alone, I compare the estimate with its own average, bin by bin. The spread is the standard deviation, over the bins, of $10\log_{10}\bigl(\hat S_x/\mathbb{E}\{\hat S_x\}\bigr)$, in dB. For each $L$ the average $\mathbb{E}\{\hat S_x\}$ is 25.1’s average periodogram with $L$ in place of $N$, so it already contains the blur, and the spread shows only the scatter.

For one periodogram, $K=1$, the exponential law of 25.1 gives a spread of 5.57 dB. Statistics derives that number; I only use it.

The record is 4096 samples of 25.1’s AR(2) noise,

with $v[n]$ white Gaussian noise of variance 1 (seed 252; measured 0.994). I run the recursion for 500 samples first and drop them, so the record starts settled. Its true PSD has a sharp peak of 22.1 dB at $\Omega=0.2998\pi$ and falls to $-9.6$ dB at $\pi$.

The picture at the top of the page cuts this record into $K$ segments of length $L$, with a rectangular window and no overlap. Its spread reads 5.51, 2.48, 1.15 and 0.50 dB for $K$ = 1, 4, 16 and 64. Its average peak is the average periodogram at the true peak, so it shows the blur without the scatter.

After the clip, drag the slider to 128 and watch the two readouts that matter. $K$ becomes 32, the spread 0.86 dB, and the average peak 21.3 dB: steadier than at 256, and blurrier.

Look at the end frame. The estimate is smooth, but its peak is only 20.6 dB on average, 1.5 dB under the truth. That is the blur. One bin at $L=64$ is $2\pi/64=0.098$ rad, about the width of the true peak itself, 0.103 rad at half its height.

At $L=256$ a bin is $2\pi/256=0.025$ rad, a quarter of the peak’s width. There the average peak, 21.7 dB, is within 0.4 dB of the truth. The valley fills too: at $L=64$ the rectangular window’s leakage lifts the average at $\pi$ from $-9.6$ dB to $-7.1$ dB.

Now the scatter. In 25.1’s terms, spread ÷ mean of $\hat S_x/\mathbb{E}\{\hat S_x\}$ is 0.967, 0.515, 0.273 and 0.112 for $K$ = 1, 4, 16 and 64. That is close to 18.2’s 1, 1/2, 1/4 and 1/8.

In dB the rule bends a little. A small change by a fraction $e$ moves a level by about $4.34e$ dB; 10 % is 0.41 dB. So once $K$ is large the spread is about $4.34/\sqrt K$ dB: 1.09 dB at $K=16$ and 0.54 dB at $K=64$, against the measured 1.15 and 0.50.

So choose $L$ from the detail you need. Make one bin, $2\pi/L$, a few times narrower than the narrowest peak you want to see; then $K=N/L$ is what the record gives you. In hertz, one bin is $f_s/L$.

Overlap buys segments

Welch’s method makes two changes to Bartlett’s recipe. First, it multiplies each segment by a window $w[n]$ before the DFT. A Hann window has side lobes from −31.5 dB instead of the rectangular window’s −13.3 dB, so it leaks far less. Its main lobe is 4 bins instead of 2, so it blurs a little more (15.2).

That window costs data, though. Hann falls to 0 at both ends, so the samples near each join between segments barely count. “Hop and overlap” in Spectrograms & the STFT (15.5) found the same for frames.

So the second change is to overlap the segments. Each one starts a hop after the last one, with $\text{hop}=L(1-\text{overlap})$. The last segment must fit inside the record, so the number of segments is

where $\lfloor\cdot\rfloor$ rounds down. With no overlap the hop is $L$, and $K=N/L$ as before.

Welch’s method averages $K$ windowed periodograms, as Bartlett’s averages plain ones. Segment $i$ now starts at sample $i\,\text{hop}$, and its windowed periodogram is

The window takes away some of the noise’s power, and dividing by $\sum w[n]^2$ instead of $L$ puts it back: for white noise the average comes out right. With $w[n]=1$ the sum is $L$, and Bartlett’s recipe returns.

SciPy computes exactly this: scipy.signal.welch(x, window='hann', nperseg=256, noverlap=128, detrend=False, return_onesided=False) gives the instrument’s estimate at 50 % overlap.

The instrument keeps the same 4096 samples and segments of 256, and changes only the overlap. Picture a relay, where each runner starts before the last one finishes.

Watch the spread as the windows move closer: 1.08 dB, then 0.76 dB, then 0.75 dB.

Notice the last step. Going from 50 % to 75 % doubles the segments, but the spread only falls from 0.76 dB to 0.75 dB.

The reason is that overlapping segments share samples, so their periodograms are no longer independent. For white noise, two Hann segments that overlap by half give periodograms with a correlation coefficient of only 0.03: nearly independent. At 75 % overlap, the periodograms of next-door segments correlate by 0.43.

Welch worked out what that does to the scatter. At 50 % overlap the 31 segments scatter like 29.4 independent ones; at 75 % the 61 scatter like only 32.0. So with Hann, 50 % overlap already gets most of what overlap can give, and more overlap adds little.

The same trade from the lag side

There is a second road to the same place, through the autocorrelation. Random processes (24.2) estimated $\hat R_x[\ell]$ from one recording, and Power spectral density (24.4) showed that the PSD is the DTFT of $R_x[\ell]$. So take the DTFT of the estimate.

Over all lags, that DTFT is called the correlogram. With the division by $N$ that 24.2 mentions in place of $N-\ell$, it equals the periodogram exactly, and it scatters just as much.

The trouble is in the far lags. 24.2’s estimate at lag $\ell$ averages only $N-\ell$ products, so the far lags are noisy. The Blackman–Tukey method multiplies the estimate by a lag window $w[\ell]$, which is 1 at $\ell=0$, falls to 0 by a chosen largest lag and stays 0 beyond it. Then it takes the DTFT:

A short lag window gives a smooth spectrum and a long one a sharp, noisy spectrum: the same trade. Bartlett’s method sits on the same road. By 25.1, its average is the DTFT of $R_x[\ell]$ tapered by the triangle $1-\lvert\ell\rvert/L$, a lag window that reaches lag $L-1$.

Multitaper estimates are a third relative: they average the periodograms of the whole record seen through several different, carefully chosen windows.

Worked example

1. Bartlett. Cutting 4096 samples into segments of 4096, 1024, 256 and 64 gives $K=4096/L$ = 1, 4, 16 and 64. The spread is 5.51 dB, 2.48 dB, 1.15 dB and 0.50 dB: each 4× in $K$ multiplies it by 0.45, 0.46 and 0.44 in turn.

2. Welch. Segments of 256 with hops of 256, 128 and 64 give

The spread is 1.08 dB, 0.76 dB and 0.75 dB.

3. Choosing $L$ for this record. The true peak is 0.103 rad wide at half its height. At $L=64$ one bin, 0.098 rad, is about as wide as the peak, and the peak loses 1.5 dB.

At $L=256$ a bin is about a quarter of the peak’s width, the average peak is 21.7 dB, and Hann with 50 % overlap gives a spread of 0.76 dB. At $f_s=1$ kHz, $L=256$ resolves about $f_s/L=3.9$ Hz.

Where you’ll meet this

Welch’s method is the usual PSD estimate in measurement software, vibration analysis and EEG, and Hann with 50 % overlap is the usual setting. SciPy’s welch uses it by default: a Hann window, segments of 256 and an overlap of half a segment.

Two of its defaults differ from this page. It removes each segment’s mean first. It also returns the positive frequencies only, doubled except at 0 and $\pi$: the one-sided density of Reading a spectrum: scaling and units (15.4). Its fs argument turns the result into units per hertz.

No average beats the limit of about $2\pi/L$. When that limit is too coarse for the record you have, other estimates help. Parametric models and linear prediction (25.3) fits a model to the data, and High-resolution frequency estimation (25.4) separates tones closer than one bin.

For more, see P. D. Welch, “The use of fast Fourier transform for the estimation of power spectra”, IEEE Trans. Audio and Electroacoustics (1967); 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.

Reference card