Skip to content

Power Spectral Density and the Wiener-Khinchin Theorem

Copy this iframe to your website:

1
<iframe src="https://dmccreary.github.io/signal-processing/sims/psd-wiener-khinchin/main.html" width="100%" height="786px" scrolling="no"></iframe>

Run MicroSim in Fullscreen

Description

The Wiener-Khinchin theorem says that the power spectral density (PSD) of a wide-sense stationary process is the Fourier transform of its autocorrelation, \(S(\omega) = \sum_k r[k]\,e^{-j\omega k}\). This MicroSim estimates the PSD of a random record in two ways and compares both with the exact \(S(\omega)\).

The Process list offers white noise, a first-order autoregressive process AR(1) \(x[n] = a\,x[n-1] + w[n]\), a moving average of \(q\) noise samples, and a sinusoid in noise. The first slider sets \(a\) (−0.95 to 0.95), \(q\) (1 to 32) or the frequency \(\omega_0\). Every process is driven by the same unit-variance white noise \(w[n]\) from a seeded generator: New realization draws new noise and Reset seed returns to seed 1, so comparisons are reproducible. Record length sets \(N\) (64 to 4096) and Averages sets the number of segments \(K\).

  • Top: the first 512 samples of \(x[n]\), with the analysis segments in gold.
  • Middle: the biased sample autocorrelation \(\hat r[k] = \frac{1}{N}\sum_n x[n]\,x[n+k]\) (blue), the theoretical \(r[k]\) (black rings) and the lag window \(c[k]\) (gold).
  • Bottom: the PSD in dB for \(0 \le \omega \le \pi\) rad/sample. Maroon: the averaged periodogram of \(K\) segments of \(M\) samples, with Hann windows and 50% overlap (Welch) or rectangular windows without overlap (Bartlett). Green: the Wiener-Khinchin estimate \(\sum_{|k|<M} c[k]\,\hat r[k]\,e^{-j\omega k}\). Dashed: the theoretical \(S(\omega)\).

The lag window is matched to the segment window, so the two estimates have nearly the same mean. With one rectangular segment, \(c[k] = 1\) at every lag and the curves coincide exactly, because \(|X(\omega)|^2/N\) is the Fourier transform of \(\hat r[k]\). The read-outs give the RMS dB difference between the estimates and each one's RMS error against theory. Point at a lag to read \(\hat r[k]\) and \(r[k]\), or at the PSD to read the three values at one frequency.

This MicroSim is used in Chapter 12: Stochastic Processes and Random Signals.

The Four Processes

With \(\sigma^2 = 1\):

Process Autocorrelation \(r[k]\) PSD \(S(\omega)\)
White noise \(x[n] = w[n]\) \(\delta[k]\) \(1\) (0 dB)
AR(1) \(x[n] = a\,x[n-1] + w[n]\) \(a^{\lvert k\rvert}/(1 - a^2)\) \(1/(1 - 2a\cos\omega + a^2) = 1/\lvert 1 - a e^{-j\omega}\rvert^2\)
MA(q) \(x[n] = (w[n] + \dots + w[n-q+1])/\sqrt{q}\) \(1 - \lvert k\rvert/q\) for \(\lvert k\rvert < q\), else 0 \(\sin^2(q\omega/2)\,/\,\big(q\sin^2(\omega/2)\big)\)
Sinusoid \(\sqrt{2}\cos(\omega_0 n + \phi) + w[n]\) \(\cos(\omega_0 k) + \delta[k]\) \(1 + \pi[\delta(\omega - \omega_0) + \delta(\omega + \omega_0)]\)

How the Two Estimates Are Computed

The averaged periodogram uses segments of \(M = \lfloor N/K \rfloor\) samples (Bartlett) or \(M = \lfloor 2N/(K+1) \rfloor\) samples with a hop of \(M/2\) (Welch), a window \(v[n]\) with \(U = \sum_n v^2[n]\), and

\[\hat S_{\text{avg}}(\omega) = \frac{1}{K}\sum_{i=0}^{K-1} \frac{1}{U}\left|\sum_{n=0}^{M-1} v[n]\,x[n + iD]\,e^{-j\omega n}\right|^2 .\]

Its mean is \(\sum_{|k|<M} c_v[k]\,r[k]\,e^{-j\omega k}\), where \(c_v[k] = \frac{1}{U}\sum_n v[n]\,v[n+|k|]\) is the autocorrelation of the window. Because \(E\{\hat r[k]\} = (1 - |k|/N)\,r[k]\), the Wiener-Khinchin estimate uses \(c[k] = c_v[k]\) for the Hann window and \(c[k] = (1 - |k|/M)/(1 - |k|/N)\) for the rectangular window. The Hann choice can never produce a negative PSD; the rectangular one can dip below zero at a few frequencies (this happens for the sinusoid with \(K\) = 2 to 5), and the read-out then says so. The RMS error against theory is measured where \(S(\omega)\) is within 40 dB of its maximum, which skips the exact zeros of the moving average.

Lesson Plan

Learning Objectives

After using this MicroSim, students will be able to:

  1. Explain the Wiener-Khinchin theorem and verify that the periodogram of one rectangular segment equals the Fourier transform of the biased sample autocorrelation
  2. Calculate \(r[k]\), \(S(0)\) and \(S(\pi)\) for AR(1) and moving-average processes from their closed forms
  3. Analyze how the number of averaged segments \(K\) trades variance against frequency resolution, using the RMS error against theory and the depth of the moving-average nulls
  4. Relate features of an autocorrelation (decay rate, alternating sign, finite length, a component that never decays) to features of the PSD
  5. Distinguish agreement between two estimates from accuracy against the true PSD

Target Audience

Undergraduate students in a digital signal processing, random signals or communications course (junior or senior year). Plan 30 minutes for the activities below.

Prerequisites

  • The DFT and decibels
  • Mean, variance and the autocorrelation of a wide-sense stationary process
  • Windowing a finite record (Chapter 7)

Activities

All numbers below are for seed 1, the seed on page load (press Reset seed to return to it). Point at a lag in the autocorrelation panel, or at a frequency in the PSD panel, to read the values there.

Activity 1: The theorem, exactly (6 min). Keep the defaults (AR(1), \(a = 0.90\), \(N = 1024\)), choose Bartlett and drag \(K\) to 1. The gold lag window becomes a flat line at 1, the green curve lies exactly on the maroon one, and the read-out reports an RMS difference of 0.00 dB: the periodogram is the Fourier transform of \(\hat r[k]\), term by term. Both curves are nevertheless far from theory (RMS error 6.35 dB). Now set \(N\) = 64, 256 and 4096: the difference stays 0.00 dB and the error stays large (6.06, 5.40 and 6.33 dB). A longer record does not make a single periodogram more accurate.

Activity 2: Averaging trades variance for resolution (10 min). Choose Welch, \(N = 4096\), and step \(K\) through 1, 3, 7, 15, 31 and 63 (so \(M\) = 4096, 2048, 1024, 512, 256, 128). The Welch error against theory falls from 6.18 to 2.85, 1.72, 1.15, 0.73 and 0.57 dB, and the two estimates converge on each other (RMS difference 4.68, 1.51, 0.75, 0.43, 0.28, 0.19 dB). Now choose the moving average with \(q = 8\) and repeat \(K\) = 15, 31, 63. The exact PSD is zero at \(\omega = \pi/4\), but the Welch estimate there reads −33.03, −26.12 and −19.40 dB: shorter segments smear power into the null. The error against theory is smallest at \(K = 31\) (0.83 dB) and rises again at \(K = 63\) (1.14 dB), where this bias outweighs the lower variance.

Activity 3: Read the spectrum from the autocorrelation (8 min). Use Welch, \(N = 1024\), \(K = 7\). With \(a = 0.90\) the autocorrelation decays slowly (theory \(r[0] = 5.26\), \(r[1] = 4.74\); the record gives \(\hat r[1] = 4.31\)) and the PSD is low-pass: theory 20.00 dB at \(\omega = 0\) and −5.58 dB at \(\pi\) (the estimates read 20.39 and 18.52 dB at \(\omega = 0\)). Set \(a = -0.90\): \(\hat r[k]\) alternates in sign (\(\hat r[1] = -4.31\)) and the PSD flips to high-pass, 20.00 dB at \(\omega = \pi\). Set \(a = 0\): \(\hat r[k]\) is a single spike and the PSD is flat at 0 dB. Choose the moving average with \(q = 8\) and \(K = 15\): \(\hat r[k]\) is a triangle that reaches zero at \(|k| = 8\), and the theoretical PSD (9.03 dB at \(\omega = 0\)) has nulls at \(\omega = \pi/4\), \(\pi/2\) and \(3\pi/4\), the multiples of \(2\pi/q\).

Activity 4: A correlation that never dies (6 min). Choose the sinusoid with \(\omega_0 = 0.25\pi\), Welch, \(N = 1024\). The autocorrelation is a cosine with a period of 8 lags that does not decay: \(\hat r[0] = 1.98\) (power 1 from the sinusoid plus 1 from the noise), \(\hat r[4] = -0.98\), \(\hat r[8] = 0.95\). The PSD shows a line on a 0 dB floor. Point at \(\omega = 0.250\pi\) and step \(K\) through 1, 3, 7, 15 and 31: the Welch peak reads 25.22, 22.21, 19.26, 16.26 and 13.35 dB, about 3 dB lower each time \(M\) is halved, while the floor stays at 0 dB. The line's power is fixed; the estimate spreads it over a main lobe of width proportional to \(1/M\).

Assessment

  1. An AR(1) process has \(a = 0.5\) and \(\sigma^2 = 1\). Find \(r[0]\), \(r[1]\), and \(S(\omega)\) at \(\omega = 0\) and \(\omega = \pi\) in dB. Answer: \(r[0] = 1/(1 - 0.25) = 1.33\) and \(r[1] = 0.67\); \(S(0) = 1/(1 - 0.5)^2 = 4\), or 6.02 dB; \(S(\pi) = 1/1.5^2 = 0.44\), or −3.52 dB. Pointing at \(\omega = 0\) in the sim shows a theory value of 6.02 dB.
  2. With Bartlett and \(K = 1\) the two curves coincide exactly. State the identity behind this. Answer: for a record of \(N\) samples, \(\frac{1}{N}|X(\omega)|^2 = \sum_{k=-(N-1)}^{N-1} \hat r[k]\,e^{-j\omega k}\) with the biased estimate \(\hat r[k] = \frac{1}{N}\sum_n x[n]\,x[n+|k|]\). It is the Wiener-Khinchin theorem applied to the data themselves.
  3. Why does the RMS error of a single periodogram stay near 6 dB when \(N\) grows from 64 to 4096, and what reduces it? Answer: a larger \(N\) gives more frequency samples but each still has a standard deviation about equal to its mean, so the estimate does not converge. Averaging \(K\) segments reduces the variance roughly \(K\)-fold, at the cost of resolution \(\propto 1/M\).
  4. For the moving average with \(q = 4\), give \(r[k]\) and the frequencies where \(S(\omega) = 0\). Why does the estimate not reach \(-\infty\) dB there? Answer: \(r[k]\) = 1, 0.75, 0.5, 0.25 for \(|k|\) = 0 to 3 and 0 beyond; nulls at \(\omega = \pi/2\) and \(\pi\). The window's main lobe and sidelobes leak power from neighbouring frequencies into the null (with \(N = 1024\) and Welch \(K = 15\) the estimate at \(\pi/2\) reads −27.40 dB).
  5. With the sinusoid, the Welch peak is 19.26 dB for \(M = 256\). Predict it for \(M = 128\), and explain why the peak height is not the sinusoid's power. Answer: about 3 dB lower; the sim shows 16.26 dB. A spectral line has a fixed power (its area); a window of \(M\) samples spreads that power over a main lobe of width about \(1/M\), so the peak height grows in proportion to \(M\).

Common Misconceptions

"A longer record gives a more accurate periodogram." In Activity 1 the single-segment error stays between 5.4 and 6.4 dB from \(N = 64\) to \(N = 4096\). Only averaging reduces it: Welch with \(N = 4096\) and \(K = 63\) reaches 0.57 dB.

"If two PSD estimates agree, they are correct." With Bartlett and \(K = 1\) the two estimates differ by 0.00 dB, yet each is 6.35 dB RMS from the true PSD. Their agreement demonstrates the Wiener-Khinchin identity, not accuracy: both are computed from the same noisy record and make the same errors.

References

  1. Wiener–Khinchin theorem - Wikipedia - the theorem that the PSD of a wide-sense stationary process is the Fourier transform of its autocorrelation
  2. Spectral density - Wikipedia - definition and properties of the power spectral density
  3. Periodogram - Wikipedia - the raw PSD estimate and why its variance does not fall as the record grows
  4. Welch's method - Wikipedia - averaging windowed, overlapping periodograms
  5. Bartlett's method - Wikipedia - averaging periodograms of non-overlapping segments, and its link to a triangular lag window
  6. Autoregressive model - Wikipedia - the AR(1) process, its autocorrelation \(a^{|k|}/(1-a^2)\) and its spectrum
  7. The Scientist and Engineer's Guide to Digital Signal Processing - Steven W. Smith - free online textbook covering the DFT, windows and noise spectra