levinson: Durbin-Levinson Recursion

In gjmvanboxtel/gsignal: Signal Processing

Durbin-Levinson Recursion

Description

Use the Durbin-Levinson algorithm to compute the coefficients of an autoregressive linear process.

Usage

levinson(acf, p = NROW(acf))

Arguments

acf	autocorrelation function for lags 0 to p, specified as a vector or matrix. If r is a matrix, the function finds the coefficients for each column of acf and returns them in the rows of a.
p	model order, specified as a positive integer. Default: NROW(acf) - 1.

Details

levinson uses the Durbin-Levinson algorithm to solve:

toeplitz(acf(1:p)) * x = -acf(2:p+1)

The solution c(1, x) is the denominator of an all pole filter approximation to the signal x which generated the autocorrelation function acf.

From ref [2]: Levinson recursion or Levinson–Durbin recursion is a procedure in linear algebra to recursively calculate the solution to an equation involving a Toeplitz matrix. Other methods to process data include Schur decomposition and Cholesky decomposition. In comparison to these, Levinson recursion (particularly split Levinson recursion) tends to be faster computationally, but more sensitive to computational inaccuracies like round-off errors.

Value

A list containing the following elements:

a

vector or matrix containing (p+1) autoregression coefficients. If x is a matrix, then each row of a corresponds to a column of x. a has p + 1 columns.

e

white noise input variance, returned as a vector. If x is a matrix, then each element of e corresponds to a column of x.

k

Reflection coefficients defining the lattice-filter embodiment of the model returned as vector or a matrix. If x is a matrix, then each column of k corresponds to a column of x. k has p rows.

Author(s)

Paul Kienzle, [e-mail],
Peter V. Lanspeary, [e-mail].
Conversion to R by Geert van Boxtel, [e-mail].

References

[1] Steven M. Kay and Stanley Lawrence Marple Jr. (1981). Spectrum analysis – a modern perspective. Proceedings of the IEEE, Vol 69, 1380-1419.
[2] https://en.wikipedia.org/wiki/Levinson_recursion

Examples

## 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 <- c(1, 0.1, -0.8, -0.27)
v <- 0.4
w <- sqrt(v) * rnorm(15000)
x <- filter(1, a, w)
xc <- xcorr(x, scale = 'biased')
acf <- xc$R[-which(xc$lags < 0)]
lev <- levinson(acf, length(a) - 1)

gjmvanboxtel/gsignal documentation built on Nov. 22, 2023, 8:19 p.m.