levinson

Levinson-Durbin recursion

Syntax

a = levinson(r,n)

[a,e,k] = levinson(___)

Description

a = levinson(r,n) returns the coefficients of an autoregressive linear process of order n that has r as its autocorrelation sequence.

[a,e,k] = levinson(___) also returns the prediction error e, and the reflection coefficients, k.

Examples

Autoregressive Process Coefficients

Estimate the coefficients of an autoregressive process given by

$$x(n)=0.1x(n-1)-0.8x(n-2)-0.27x(n-3)+w(n).$$

a = [1 0.1 -0.8 -0.27];

Generate a realization of the process by filtering white noise of variance 0.4.

v = 0.4;
w = sqrt(v)*randn(15000,1);
x = filter(1,a,w);

Estimate the correlation function. Discard the correlation values at negative lags. Use the Levinson-Durbin recursion to estimate the model coefficients. Verify that the prediction error corresponds to the variance of the input.

[r,lg] = xcorr(x,'biased');
r(lg<0) = [];

[ar,e] = levinson(r,numel(a)-1)

ar = 1×4

    1.0000    0.0772   -0.7954   -0.2493

e = 0.3909

Estimate the reflection coefficients for a 16th-order model. Verify that the only reflection coefficients that lie outside the 95% confidence bounds are the ones that correspond to the correct model order. See AR Order Selection with Partial Autocorrelation Sequence for more details.

[~,~,k] = levinson(r,16);
stem(k,'filled')

conf = sqrt(2)*erfinv(0.95)/sqrt(15000);
hold on
[X,Y] = ndgrid(xlim,conf*[-1 1]);
plot(X,Y,'--r')
hold off

Prediction Errors for Multiple Realizations

Generate the coefficients of an autoregressive process given by

$$x(n)=0.1x(n-1)-0.8x(n-2)-0.27x(n-3)+w(n).$$

a = [1 0.1 -0.8 -0.27];

Generate five realizations of the process by filtering white noise with different variances.

nr = 5;
v = rand(1,nr)

v = 1×5

    0.8147    0.9058    0.1270    0.9134    0.6324

w = sqrt(v).*randn(15000,nr);
x = filter(1,a,w);

Estimate the correlation function. Discard cross-correlation terms and correlation values at negative lags. Use the Levinson-Durbin recursion to estimate the prediction errors for the correct model order and verify that the prediction errors correspond to the variances of the input noise signals.

[r,lg] = xcorr(x,'biased');

[~,e] = levinson(r(lg>=0,1:nr+1:end),numel(a)-1)

e = 5×1

    0.7957
    0.9045
    0.1255
    0.9290
    0.6291

Input Arguments

r — Autocorrelation sequence

vector | matrix

Autocorrelation sequence, specified as a vector or matrix. If r is a matrix, the function finds the coefficients for each column of r and returns them in the rows of a.

Example: [r,lg] = xcorr(randn(1000,1),'biased'); r(lg<0) = [] estimates the autocorrelation sequence of a 1000-sample random signal for positive lags.

Data Types: single | double
Complex Number Support: Yes

n — Model order

length(r) – 1 (default) | positive integer scalar

Model order, specified as a positive integer scalar.

Data Types: single | double

Output Arguments

a — Autoregressive linear process coefficients

row vector | matrix

Autoregressive linear process coefficients, returned as a row vector or matrix. The filter coefficients are ordered in descending powers of z^–1:

$$H(z)=\frac{1}{A(z)}=\frac{1}{1+a(2)z^{-1}+\cdots +a(n+1)z^{-n}}.$$

If r is a matrix, then each row of a corresponds to a column of r.

e — Prediction error

scalar | column vector

Prediction error, returned as a scalar or column vector. If r is a matrix, then each element of e corresponds to a column of r.

k — Reflection coefficients

column vector | matrix

Reflection coefficients, returned as a column vector of length n. If r is a matrix, then each column of k corresponds to a column of r.

Note

k is computed internally while computing the a coefficients, so returning k simultaneously is more efficient than converting a to k with tf2latc.

Algorithms

The Levinson-Durbin recursion is an algorithm for finding an all-pole IIR filter with a prescribed deterministic autocorrelation sequence. It has applications in filter design, coding, and spectral estimation. The filter that levinson produces is minimum phase.

levinson solves the symmetric Toeplitz system of linear equations

$$\left[\begin{array}{cccc}r(1)& r{(2)}^{*}& \cdots & r{(n)}^{*}\\ r(2)& r(1)& \cdots & r{(n-1)}^{*}\\ ⋮& ⋮& \ddots & ⋮\\ r(n)& \cdots & r(2)& r(1)\end{array}\right]\left[\begin{array}{c}a(2)\\ a(3)\\ ⋮\\ a(n+1)\end{array}\right]=\left[\begin{array}{c}-r(2)\\ -r(3)\\ ⋮\\ -r(n+1)\end{array}\right],$$

where r = [r(1) ... r(n + 1)] is the input autocorrelation vector, and r(i)^* denotes the complex conjugate of r(i). The input r is typically a vector of autocorrelation coefficients where lag 0 is the first element, r(1).

Note

If r is not a valid autocorrelation sequence, the levinson function might return NaNs even if the solution exists.

The algorithm requires O(n²) flops and is thus much more efficient than the MATLAB^® backslash command for large n. However, the levinson function uses \ for low orders to provide the fastest possible execution.

References

[1] Ljung, Lennart. System Identification: Theory for the User. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• Code generation for this function requires the DSP System Toolbox™ software.

• Input n, when specified, must be a scalar.