Power spectral density

A random signal keeps its power somewhere along the frequency axis. Here AR(1) noise, whose samples remember a share $a$ of their neighbour, is drawn twice: its autocorrelation, and its power at each frequency. Watch both as $a$ grows from 0: the longer the memory, the narrower the spectrum.

Power per frequency

In Random processes (24.2), “White noise forgets, coloured noise remembers” measured how alike two samples $\ell$ apart are. That was the autocorrelation, $R_x[\ell]=\mathbb{E}\{x[n]\,x[n-\ell]\}$. White noise had none beyond lag 0, and AR(1) noise kept a share $a$ of its likeness per step. On this page I ask the same question in frequency: where along the frequency axis does a random signal keep its power?

Let’s start with a signal whose answer we know: a sine. In “A test sine goes in, a scaled and shifted sine comes out” of Frequency response of discrete-time systems (12.4), a sine of amplitude $A$ came out of a system $\lvert H(e^{j\Omega})\rvert$ times as big. Its average power, from How big is a signal (1.3), was $A^2/2$ going in. Coming out it is $\lvert H\rvert^2A^2/2$. So a system multiplies the power at each frequency by $\lvert H\rvert^2$.

A noise has no single frequency. Its power is spread over the band, as the quantizer’s error was in “Sharing the same error among more slices” of Oversampling and noise shaping (11.3). There I said that this error’s power sits evenly from 0 to $f_s/2$, and promised the proof for this page.

To describe a spread I need a density: the power per unit of frequency near $\Omega$. That is the power spectral density, or PSD, written $S_x(e^{j\Omega})$. Adding it up over the whole band gives the whole power. For a process of mean 0, like every one on this page, that power is $R_x[0]$, the average power $P$ of 1.3:

Now filter the noise. Each narrow slice of frequency goes through like a test sine, so its power is multiplied by $\lvert H\rvert^2$. The output’s PSD is

This needs a stable system and an input that has been running for a long time, 12.4’s steady state. The angle $\angle H$ does not appear: delaying a noise does not move its power to other frequencies.

Wiener–Khinchin: the PSD is the DTFT of the autocorrelation

Where does $S_x$ come from? Take $N$ samples of the process and their DFT $X[k]$. In “Energy in time and in bins” of The DFT (13.2), the energy of the samples was $\frac1N\sum_k\lvert X[k]\rvert^2$. Divide by $N$ once more and you have the power, the average over the $N$ bins of $\lvert X[k]\rvert^2/N$. So $\lvert X[k]\rvert^2/N$ is the power that sits at bin $k$.

What is that number on average? Write it at any $\Omega$, with $X_N(e^{j\Omega})=\sum_{n=0}^{N-1}x[n]e^{-j\Omega n}$, the DTFT of the $N$ samples. Multiply out the square and take the expected value of “Mean and variance” in Random variables for signals (24.1). Each product $x[n]\,x[m]$ becomes $R_x[n-m]$, and the $N-\lvert\ell\rvert$ pairs with the same lag $\ell=n-m$ give the same term:

As $N$ grows the weights $1-\lvert\ell\rvert/N$ tend to 1, and what is left is the DTFT of the autocorrelation, from The DTFT (12.2):

This is the Wiener–Khinchin theorem. Read backwards, the inverse DTFT gives $R_x[\ell]$ from $S_x$, and at $\ell=0$ it is the power formula above. An autocorrelation is even, $R_x[-\ell]=R_x[\ell]$, so the sum is $R_x[0]+2\sum_{\ell\ge1}R_x[\ell]\cos(\Omega\ell)$: a real curve, and an even one.

White noise $v[n]$ of variance $\sigma_v^2$ has $R_v[\ell]=\sigma_v^2\,\delta[\ell]$. Only the $\ell=0$ term survives, so

The PSD is flat: equal power at every frequency. That is the promise of 11.3 kept, for an error that behaves like white noise.

The AR(1) pair

Now the coloured noise of 24.2, AR(1) scaled to variance 1, $x[n]=a\,x[n-1]+\sqrt{1-a^2}\,v[n]$ with $\sigma_v^2=1$. Its autocorrelation is $R_x[\ell]=a^{\lvert\ell\rvert}$.

Split the sum at $\ell=0$. The lags $\ell\ge0$ give the geometric series of “One arrow, two curves” in 12.2, the DTFT of $a^nu[n]$, which settles because $\lvert a\rvert$ is less than 1. The negative lags give the same series in $e^{+j\Omega}$, starting one step later:

Put both over the common denominator $(1-ae^{-j\Omega})(1-ae^{j\Omega})=1-2a\cos\Omega+a^2$. The numerator is $1-ae^{j\Omega}+ae^{j\Omega}-a^2=1-a^2$, so

There is a second road to the same answer. AR(1) noise is white noise of variance $1-a^2$ sent through the recursion $1/(1-ae^{-j\Omega})$, so $S_y=\lvert H\rvert^2S_x$ gives $(1-a^2)/\lvert1-ae^{-j\Omega}\rvert^2$, the same curve. At $\Omega=0$ the formula is $(1+a)/(1-a)$, and at $\Omega=\pi$ it is $(1-a)/(1+a)$.

Memory in time, narrowness in frequency

The picture at the top of the page draws $R_x[\ell]$ and its DTFT together for AR(1) noise, as $a$ grows from 0. Its decibels are $10\log_{10}S_x$, as for any power in “Signal and noise, in decibels” (1.3). Think of a crowd: a crowd that sways slowly hums low, and a jittery one hisses.

Notice the end caption’s last sentence. The power is the same, 1, at every $a$: it only moves in frequency. The area is that of $S_x$ itself; in dB the curve is mirrored about 0 dB, because $S_x(e^{j0})\,S_x(e^{j\pi})=1$.

Where does the curve cross the white level? $S_x=1$ when $1-a^2=1-2a\cos\Omega+a^2$, that is, when $\cos\Omega=a$. For $a=0.5$ that is $\Omega=\pi/3$, and for $a=0.9$ it is $0.144\pi$. Below that frequency the power is raised, and above it the power is cut.

How narrow is narrow? At $a=0.9$, the slowest tenth of the band, $\lvert\Omega\rvert$ below $0.1\pi$, holds 79.6 % of the power; white noise keeps 10 % there.

The PSD falls to half its peak at $\Omega=0.105$ rad/sample. The autocorrelation falls to $1/e$ of its first value after $-1/\ln0.9=9.49$ lags. Their product, $9.49\times0.105$, is 1.00: the longer the memory, the narrower the spectrum.

After the clip, use the slider to try a negative $a$. At $-0.5$ the PSD is −4.8 dB at 0 and 4.8 dB at π: memory that alternates in sign moves the power up. Neighbours now tend to have opposite signs, so the power gathers near $\pi$. At either end of the slider the PSD reaches ±15.9 dB, still inside the panel.

Does a recording agree with the formula? I drew 400 realisations of 256 samples each, with $a=0.9$, from the site’s seeded generator (seed 244). For each I took $\lvert X[k]\rvert^2/256$ and then averaged the 400 results, bin by bin. At every bin from 1 to 127 the average is within 0.61 dB of $S_x$.

At $\Omega=0$ the average reads 0.900 of $S_x$. The peak there is so sharp that 256 samples blur it: the weights $1-\lvert\ell\rvert/N$ above predict 0.963. The rest is the scatter of 400 random results, about 5 % for one bin. Estimating a PSD from data is the subject of The periodogram (25.1).

Equal power in every octave

Our ears judge pitch by ratios. From 20 Hz to 40 Hz is one octave, a doubling, as in 11.3, and so is 10 kHz to 20 kHz. The range we hear, 20 Hz to 20 kHz, is just under 10 octaves.

White noise has the same power in every hertz. Its octave from 10 to 20 kHz is 10 000 Hz wide, and its octave from 20 to 40 Hz only 20 Hz wide. So the top octave holds 500 times the power, 27.0 dB more, and white noise sounds like a hiss.

Pink noise has the same power in every octave instead: 20 to 40 Hz holds as much as 10 to 20 kHz. Each octave is twice as wide as the one before it, so to hold the same power its PSD must be half as high. Half the power is $10\log_{10}\tfrac12=-3.01$ dB, so pink noise falls 3.01 dB per octave.

Brown noise is white noise added up, a random walk: $s[n]=s[n-1]+v[n]$. The first difference of 11.3, $1-e^{-j\Omega}$, has gain $2\lvert\sin(\Omega/2)\rvert$, as 12.4 worked out, and adding up undoes a first difference. So the running sum has gain $1/(2\lvert\sin(\Omega/2)\rvert)$, and

For small $\Omega$, $\sin(\Omega/2)$ is close to $\Omega/2$, so $S_s$ is close to $\sigma_v^2/\Omega^2$. Each doubling of frequency then divides the power by 4: $10\log_{10}\tfrac14=-6.02$ dB per octave.

A pure random walk wanders off: its variance grows with $n$, so it is not stationary, and its PSD is infinite at $\Omega=0$.

In practice I let the sum leak a little, $s[n]=0.995\,s[n-1]+v[n]$, the leaky integrator of Difference equations (6.1) with a different constant. That is AR(1) noise with $a=0.995$, a longer memory than the first instrument reaches. Its PSD is flat below a corner of about $(1-0.995)$ rad/sample and falls about 6 dB per octave above it. At $f_s=44.1$ kHz the corner is $0.005\cdot44\,100/2\pi=35.1$ Hz.

No simple recursion falls exactly 3 dB per octave. J. O. Smith’s pink-noise filter comes close with three poles and three zeros. They take turns along the real axis between 0.1 and 1. The falls of the poles and the rises of the zeros overlap, and on average the slope is near −3 dB per octave.

White, pink, brown

The second instrument passes white noise at 44.1 kHz through the two filters and draws the three PSDs. The frequency axis is logarithmic: each octave takes the same width, so a fixed number of dB per octave is a straight line. Rain on leaves sounds pink, a distant waterfall brown, and radio static white.

Notice that pink measures −2.94 dB per octave, not exactly −3.01. The filter only approximates the ideal slope, but over the band we hear it never strays from it by more than 0.67 dB.

Brown’s slope eases near the top, to −5.46 dB from 5 to 10 kHz and −3.60 dB from 10 to 20 kHz. $\lvert H\rvert$ is even and repeats every $2\pi$, so it is mirrored about $\Omega=\pi$, which is 22.05 kHz here, and every curve levels off as it nears that frequency.

After the clip, choose each colour and press “Hear it”: white hisses, pink sounds even, and brown rumbles.

Every colour here came from the same white noise. Only $\lvert H\rvert^2$ differed, and $S_y=\lvert H\rvert^2S_x$ turned it into the colour. To make noise of a given shape, design a filter whose gain squared has that shape, as closely as a filter can, and feed it white noise.

Worked example

1. AR(1) PSD. For $a=0.5$, $S_x$ at $\Omega=0$ is $1.5/0.5=3.00$, which is 4.77 dB, and at $\pi$ it is $0.5/1.5=0.3333$, which is $-4.77$ dB. For $a=0.9$ the two values are $1.9/0.1=19.00$ and $0.1/1.9=0.0526$, or 12.79 dB and $-12.79$ dB.

2. Through a filter. 12.2’s signal $0.8^nu[n]$ had $\lvert X\rvert=5$ at $\Omega=0$ and 0.556 at $\pi$. Used as a filter on white noise of variance $1-0.8^2=0.36$, it gives $S_y=25\times0.36=9$ at 0 and $0.36/1.8^2=0.111$ at $\pi$. The AR(1) formula with $a=0.8$ gives the same, $1.8/0.2=9$ and $0.2/1.8=0.111$, or 9.54 dB and $-9.54$ dB.

3. Octaves. Pink noise falls $10\log_{10}2=3.01$ dB per octave and brown noise $20\log_{10}2=6.02$ dB. Over the $\log_2 1000=9.97$ octaves from 20 Hz to 20 kHz, pink falls 30.0 dB, which is $10\log_{10}1000$, and brown falls 60.0 dB.

4. Units. The band from $-\pi$ to $\pi$ is $f_s$ hertz wide, so in hertz the two-sided density is $S_x(f)=S_x(e^{j\Omega})/f_s$, with $\Omega=2\pi f/f_s$. For white noise of variance 1 V² at $f_s=8$ kHz that is $1.25\times10^{-4}$ V²/Hz. The one-sided density of “One spectrum, four units” in Reading a spectrum: scaling and units (15.4) doubles it, to $2.5\times10^{-4}$ V²/Hz, or $-36.0$ dB re 1 V²/Hz.

Where you’ll meet this

Sound engineers play pink noise through loudspeakers and read the room with an analyser whose bands are an octave or a third of an octave wide. Pink noise puts equal power in every such band, so a flat reading means a flat system.

Random walks model slow drift: a gyroscope’s bias, a clock’s error, the position of an object pushed by random forces. The Kalman filter of RLS and the Kalman filter (26.4) builds such coloured noise into its model. Noise shaping, in “What a first difference does to slow and fast wiggles” of 11.3, is a designed $\lvert H\rvert^2$. The quantizer’s white error goes through $1-e^{-j\Omega}$, and its PSD becomes $4\sin^2(\Omega/2)$ times the flat one.

Estimating a PSD from a recording comes next, in The periodogram (25.1) and Averaged periodograms: Bartlett and Welch (25.2). Parametric models and linear prediction (25.3) fits AR models and reads their PSD from the formula. The Wiener filter (26.2) builds a filter from the PSDs of a signal and of its noise.

Reference card