levinson

Levinson-Durbin recursion

Syntax

a = levinson(r)
a = levinson(r,n)
[a,e] = levinson(r,n)
[a,e,k] = levinson(r,n)


Description

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.

a = levinson(r) finds the coefficients of a length(r)-1 order autoregressive linear process which has r as its autocorrelation sequence. r is a real or complex deterministic autocorrelation sequence. If r is a matrix, levinson finds the coefficients for each column of r and returns them in the rows of a. n=length(r)-1 is the default order of the denominator polynomial A(z); that is, a = [1 a(2) ... a(n+1)]. 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}}$$

a = levinson(r,n) returns the coefficients for an autoregressive model of order n.

[a,e] = levinson(r,n) returns the prediction error, e, of order n.

[a,e,k] = levinson(r,n) returns the reflection coefficients k as a column vector of length n.

Note

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

Examples

Autoregressive Process Coefficients

Estimate the coefficients of an autoregressive process given by

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

a = [1 0.1 -0.8];

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

rng('default')
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.

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

ar = levinson(r,numel(a)-1)

ar = 1×3

    1.0000    0.0957   -0.8026

Algorithms

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 & \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.