Parametric models and linear prediction

Here 4096 samples of noise, made by a filter with four coefficients, are fitted one order at a time: each order predicts every sample from one more sample before it. Watch the error power left over: it falls until order 4, then stays near 1.

A few numbers instead of a whole spectrum

In The periodogram (25.1), “More data, same scatter” showed the trouble with estimating a PSD bin by bin. Each bin scatters about the truth by about the truth itself, and more data only gives more bins.

But much of the noise we meet is white noise through a filter with a handful of coefficients. If I know that, I only need to estimate those few numbers. Then the filter hands me the whole spectrum, smooth, through 24.4’s rule.

Here is the model. An autoregressive model of order $p$, AR($p$) for short, makes each sample from the $p$ before it plus fresh white noise:

The number $p$ is the model order. The coefficients $a_k$ are the output-side coefficients of a difference equation, as in Difference equations (6.1), with $a_0=1$. So $x[n]$ is white noise $v[n]$, of variance $\sigma_v^2$, through the all-pole filter $1/A(z)$, where $A(z)=1+a_1z^{-1}+\dots+a_pz^{-p}$.

The AR(1) noise of “White noise forgets, coloured noise remembers” in Random processes (24.2), $x[n]=a\,x[n-1]+v[n]$, is the case $p=1$ with $a_1=-a$. The minus sign is the difference-equation habit; this page keeps it.

What is its spectrum? In Power spectral density (24.4), “Memory in time, narrowness in frequency” gave two facts. White noise has a flat PSD, $\sigma_v^2$, and a filter multiplies a PSD by $\lvert H\rvert^2$. Here $H=1/A$, so

Every pole pair of $1/A(z)$ near the unit circle makes a peak at its angle, as “Where the pair sits is how h[n] moves” showed in Transfer functions, poles & zeros (16.3). A pair takes two coefficients, so $p$ numbers can draw up to $p/2$ such peaks.

There are two other models. A moving-average model, MA, is white noise through an FIR filter: $x[n]=b_0v[n]+b_1v[n-1]+\dots$, with zeros instead of poles. An ARMA model has both, and its PSD is $\sigma_v^2\lvert B(e^{j\Omega})\rvert^2/\lvert A(e^{j\Omega})\rvert^2$. This page fits AR models, because their coefficients come from linear equations.

Predict the next sample

How do I find the $a_k$ from data? Turn the model round. If $x[n]$ is made from its past plus something new, I can guess it from its past. The prediction from the $p$ samples before is

and the prediction error is what is left:

That is $x[n]$ through the FIR filter $A(z)$. If $x$ really is AR($p$) with these $a_k$, the error is exactly $v[n]$: the part nobody could have guessed.

I choose the $a_k$ that make the mean square of the error as small as possible. As in “Least squares or least worst” in Optimal FIR design (19.3), the error is linear in the coefficients, so its mean square is a quadratic. Its lowest point comes from $p$ linear equations.

Here is a way to see those equations without calculus. Suppose the error were still correlated with one of the past samples, say ${x[n-i]}$. Then I could add a little more of ${x[n-i]}$ to the prediction and shrink the error. So at the best coefficients the error is uncorrelated with each of the $p$ past samples:

Averages add (24.2), and each product averages to an autocorrelation: $\mathbb{E}\{x[n]\,x[n-i]\}=R_x[i]$ and $\mathbb{E}\{x[n-k]\,x[n-i]\}=R_x[i-k]$. So for $i=1,\dots,p$,

These are the Yule–Walker equations. For $p=3$ they read

where I used $R_x[-\ell]=R_x[\ell]$ from 24.2. In general they are $\mathbf{R}\mathbf{a}=-\mathbf{r}$. The $p\times p$ matrix $\mathbf{R}$ has $R_x[i-k]$ in row $i$, column $k$; $\mathbf{a}$ holds $a_1$ to $a_p$, and $\mathbf{r}$ holds $R_x[1]$ to $R_x[p]$.

Look at the matrix. Every diagonal holds one value: $R_x[0]$ on the main one, $R_x[1]$ beside it, and so on. A matrix that is constant along each diagonal is called Toeplitz.

The error power left at the best coefficients is worth knowing too. The error is uncorrelated with the past samples, so its mean square equals the mean of the error times $x[n]$ alone:

With data, I replace each $R_x[\ell]$ by an estimate from $N$ samples. 24.2 divided the $N-\ell$ products at lag $\ell$ by their number. Here I divide by $N$:

The reason is stability. With $\frac1N$, $\mathbf{R}$ is positive definite: every weighted sum of samples gets a positive power from it, as a real power must. That keeps the error powers of the next section positive, and the fitted model stable. Dividing by $N-\ell$ does not promise it.

Solving it one order at a time

Elimination solves $p$ equations in about $p^3$ steps. The Levinson–Durbin recursion uses the Toeplitz pattern to do it in about $p^2$. It solves the order-1 problem, then grows the answer to order 2, 3, and on up to $p$.

It starts from no prediction at all, with error power $E_0=R_x[0]$, the whole power. Going from order $m-1$ to order $m$ takes three lines. First the reflection coefficient

Then the new coefficients, for $k=1,\dots,m-1$, all computed from the old ones at once, and one more:

Then the new error power:

The top of $k_m$ has a meaning. It is the average of the order-$(m-1)$ error times ${x[n-m]}$, the one sample the shorter predictor did not use. If the error has nothing in common with it, $k_m=0$, the coefficients stay as they were, and so does the error power.

The name comes from Filter structures (21.1). In “Three more ways to build it”, the lattice of $1/A(z)$ had $k_2=a_2$ and $k_1=a_1/(1+a_2)$. That is this recursion at order 2. It sets $a_2=k_2$ and turns $a_1=k_1$ into $k_1+k_2k_1=k_1(1+k_2)$, so 21.1’s formulas are the recursion read from the $a$‘s back to the $k$‘s.

21.1 also showed that the biquad’s stability triangle of Stability and causality (16.4), “The triangle is the inside of the circle”, is $\lvert k_1\rvert<1$ and $\lvert k_2\rvert<1$. The rule holds at every order: $1/A_p(z)$ is stable when every $\lvert k_m\rvert<1$. A positive definite $\mathbf{R}$ keeps every $E_m$ positive, and $E_m=E_{m-1}(1-k_m^2)$ then forces each $\lvert k_m\rvert<1$.

Try it on 24.2’s AR(1) noise with $a=0.9$, scaled to variance 1, so $R_x[\ell]=0.9^{\lvert\ell\rvert}$. Order 1 gives $k_1=-0.9/1=-0.9$ and $E_1=1\cdot(1-0.81)=0.19$, the driving noise’s power. Order 2 gives $k_2=-(0.81+(-0.9)(0.9))/0.19=0$. The second sample back adds nothing, because the noise is AR(1).

The picture at the top of the page runs the recursion on a known AR(4) process. Its poles are two pairs, $0.95$ at angles $\pm0.2\pi$ and $0.9$ at $\pm0.45\pi$, and $\sigma_v^2=1$. Multiplying out the two biquads gives the coefficients $a$ = 1, −1.8187, 2.1453, −1.4992, 0.7310. Its true PSD has peaks at 0.20π and 0.44π. The second peak sits a little below its pole angle, because it rides on the downward slope of the first.

One order at a time

That picture feeds the $\frac1N$ estimates $\hat R_x[0]$ to $\hat R_x[8]$ of 4096 samples to the recursion, and shows each order as it is added: a stem for $k_m$ and a bar for $E_m$. Its dotted level marks the driving noise’s power, 1, which measures 1.003 for the draws used.

Watch the bars. Each order shrinks the error power, from 13.350 to 1.041 at order 4. Then it stops: from order 4 on, $E_m$ stays between 1.040 and 1.041, close to the driving noise’s power. Past the true order the new $k_m$ are close to 0, the “nothing in common” case.

Notice also that order 3 helps little, from 2.607 to 2.273, and order 4 a lot. So look for where the error power stops falling for good, not for its first pause.

The estimated order-4 model is $a$ = 1, −1.8052, 2.1355, −1.4928, 0.7362, against the true 1, −1.8187, 2.1453, −1.4992, 0.7310. SciPy’s solve_toeplitz, which solves $\mathbf{R}\mathbf{a}=-\mathbf{r}$ directly, gives the same numbers to about $10^{-15}$, which is rounding error.

Every $\lvert k_m\rvert$ is at most 0.778, so the model is stable: 16.4’s triangle, at order 4. Its poles are 0.952 at $\pm0.202\pi$ and 0.902 at $\pm0.452\pi$, close to the true 0.95 at $\pm0.2\pi$ and 0.9 at $\pm0.45\pi$.

A smooth model against a ragged periodogram

Now the payoff: the spectrum. Once I have the $a_k$ of order $p$, the AR estimate of the PSD is the model’s PSD, with the final error power $E_p$ standing in for $\sigma_v^2$:

Here $A_p$ is the order-$p$ polynomial, and the hat means “estimated”, as in 24.1. The next instrument draws it against the periodogram of 25.1, $P_N(e^{j\Omega})=\frac1N\lvert X(e^{j\Omega})\rvert^2$, from the same samples: the first 1024 samples of the record in “One order at a time”, with the fit made from their own autocorrelation.

A peak here is a local maximum that stands at least 2 dB above its surroundings. The distance from the truth is the RMS, over Ω from 0 to π, of the difference in dB between the AR estimate and the true PSD.

Watch the solid curve as the order goes from 2 to 4 to 20: one peak between the true two, then both, then both with a bump the truth does not have. More order is not better.

When the clip ends, use the slider to try a few orders. At order 1 there is no peak at all, and the distance is 12.1 dB. Order 3 still finds one peak, at 0.30π, 6.0 dB away. From order 4 the two peaks are there: orders 8 and 12 sit 1.5 dB away and order 16 1.6 dB, a little worse than order 4’s 1.3 dB.

So too low an order merges the two peaks into one, between them. The right order finds both. A much higher order still finds them, but spends its spare coefficients on small bumps the truth does not have.

The periodogram of the same 1024 samples is 6.1 dB away from the truth, worse than any AR estimate from order 4 to 20. Its dots carry no assumption about the shape, and they pay for it in scatter. The AR curve assumes a few poles, and gains a smooth line when that assumption is right.

How do I choose the order with real data, where there is no truth to compare with? I look where the error power of “One order at a time” stops falling. Criteria such as Akaike’s AIC make this a rule: they add a penalty for each coefficient to the logarithm of the error power, and pick the order where the sum is smallest.

The AR estimate needs enough data too, because it is built from estimated correlations. With only the first 256 samples, the order-4 fit’s second peak, near 0.44π, stands just 0.5 dB above its surroundings.

Worked example

1. Levinson by hand, two orders. Take $R_x[0]=1$, $R_x[1]=0.5$, $R_x[2]=0.1$. Order 1: $k_1=-0.5/1=-0.50$, $E_1=1\cdot(1-0.25)=0.75$ and $a_1=-0.50$.

Order 2: $k_2=-(0.1+(-0.5)(0.5))/0.75$, which is $0.15/0.75=0.20$, and $E_2=0.75\cdot(1-0.04)=0.72$. The new $a_1=-0.5+0.2\cdot(-0.5)=-0.60$, and $a_2=0.20$.

Check both Yule–Walker equations: $-0.6\cdot1+0.2\cdot0.5=-0.5=-R_x[1]$ and $-0.6\cdot0.5+0.2\cdot1=-0.1=-R_x[2]$. The error power checks too: $1+(-0.6)(0.5)+(0.2)(0.1)=0.72$.

2. The clip’s data. In “One order at a time” the error power falls from 13.350 to 1.041 at order 4 and then stays flat, 1.041 and 1.040. That is near the driving noise’s measured power, 1.003, so the order is 4.

Where you’ll meet this

Speech is the classic case. A voice is a buzz or a hiss shaped by the throat and mouth, and over a short frame that shaping is close to an all-pole filter. Speech coders at 8 kHz fit an AR model of order 8 to 10 to each frame of 10 to 25 ms, and send its coefficients instead of the samples. The source–filter model of speech (30.1) builds this.

Some coders send the reflection coefficients, or numbers made from them. A $k_m$ is easy to keep inside $(-1,1)$ after rounding, so the decoded filter stays stable.

AR estimates also suit short records, where a periodogram cannot separate close peaks. High-resolution frequency estimation (25.4) takes that further. The same Toeplitz system, with a cross-correlation on the right, returns as the Wiener–Hopf equations of The Wiener filter (26.2).

For more, see J. Makhoul, “Linear prediction: a tutorial review” (Proc. IEEE, 1975); J. G. Proakis and D. G. Manolakis, Digital Signal Processing, ch. 12 and 14; and P. Stoica and R. Moses, Spectral Analysis of Signals (2005), ch. 3. In SciPy, solve_toeplitz(r[:p], -r[1:p+1]) solves Yule–Walker, with r holding the estimates $\hat R_x[0]$ to $\hat R_x[p]$.

Reference card