Advanced DSP

Levinson-Durbin Algorithm

Suppose we want to least-squares solve a convolution system Ax ≈ y. The columns of A are shifted variants of the same “signal” a. Thus, these column vectors, a₁, a₂, ...a_l have special properties. E.g. a₂ looks to a₁ as a₃ looks to a₂, as a₂ looks to a₃. The dot products between all these columns can be found by finding the autocorrelation of a, which can be found using the FFT. Similarly, the dot product of the columns with y is the cross correlation between a and y.

Suppose we already know the best prediction of y using [a₁, a₂, a₃], and we know the linear combiners which will achieve this. Thus, we know y ≈ β1a₁ + β2a₂ + β3a₃. To find the best prediction of y in terms of [a₁, a₂, a₃, a₄], we first orthogonalize a₄ with respect to [a₁, a₂, a₃], as follows.

Suppose we already know the best prediction of a₃ with respect to [a₁, a₂]. Suppose a₃ ≈ α₁a₁ + α₂a₂. We can also say a₄ ≈ α₁a₂ + α₂a₃. Immediately, a₄ can be orthogonalized with respect to [a₂, a₃]. Now if we orthogonalized a₁ with respect to [a₂, a₃] our task would be done. Ah, but a₁ ≈ α₁a₃ + α₂a₂. Thus, the orthogonalization of a₁ is q = a₁ − α₁a₃ − α₂a₂. We should use exactly ‹a₄,q›/‹q,q› scale of q. Out of this, ‹a₄,q› is found by substituting the value of q and using the autocorrelation values of a. ‹q,q› is given to us by the previous iteration. Multiply the coefficients in the q expression by this scale, and add it to the previous a₄ approximation expression, to get the coefficients which give the best approximation to a₄ in terms of [a₁, a₂, a₃]. Thus, we get the orthogonalization of a₄ with respect to [a₁, a₂, a₃] which we will call q₄ = a₄ − α₁'a₁ − α₂'a₂ − α₃'a₃.

This means my new best estimation of y is y ≈ β1a₁ + β2a₂ + β3a₃ + (‹y,q₄›/‹q₄,q₄›)q₄. Substituting the q₄ expression in the above expression will give us the required prediction of y in terms of [a₁, a₂, a₃, a₄].

There is a fast way to calculate ‹q₄,q₄›. ‹q₄,q₄› is the energy remaining in a₄ after the energy of the vector (‹a₄,q›/‹q,q›)q is taken away from the original remaining energy of ‹q,q›. Thus ‹q₄,q₄› = ‹q,q› − (‹a₄,q›)²/‹q,q›. Furthermore, ‹q₄,q₄› becomes the ‹q,q› for the next iteration.

Thus, we end up with the αs, βs and ‹q,q›, which are used in the next iteration. Each iteration takes time proportional to current order. If the final order is l, the algorithm runs in O(l²) time. The time to find the auto and cross correlations will O(ologo), where o is the dimensions of a and y. It is important to notice that once the correlations are found out, there is no more meddling in the signal space. We are directly in the space of x. All the equations above are evaluated only upto their coefficients in the a_•s, not till the final vector in the y space. The final βs are the required x.

The Levinson-Durbin algorithm is used to solve LS matrix inversion problems for Toeplitz matrices – LS deconvolution problems, (which themselves occur in many system identification algorithms) and prediction, Wiener filtering and autocorrelation estimation of stationary random processes. This algorithm is based on the idea of successive orthogonalization.