Introduction
I’ve recently been playing with this algorithm and wanted to do a write-up of my understanding on how it works. Better references certainly exist, but I wanted to write down how my brain thinks about the algorithm if for no other reason than for my own future reference.
I’m going to focus on the algorithm and explain how to build it - it won’t be the most efficient version, but will hopefully demonstrate the concepts.
Toeplitz systems
The Levinson-Durbin algorithm finds the solution to a Toeplitz system which has a non-zero main diagonal:
Toeplitz systems arise frequently in signal processing problems and, whenever I see them, I think “convolution”. For example, given a discrete convolution:
Where
Toeplitz matrix multiplications can be implemented using convolutions!
The algorithm
Start by solving a 2x2 Toeplitz system using Cramer’s rule:
To move to higher order systems, let’s define some new quantities:
is the -by- Topelitz matrix that is taken as the first rows and columns of a larger Toeplitz matrix. In the above 2x2 example:
is the -row column vector that when multiplied by produces (i.e. the column vector which has a 1 in the top row and 0 in all other rows). See , think upper - it is orthogonal to all but the top row of .
is the -row column vector that when multiplied by produces (i.e. the column vector with a 1 in the bottom row and 0 in all other rows). See , think lower - it is orthogonal to all but the bottom row of .
is the top rows of the column vector being solved. is the top rows of the column vector that satisfies = .
Let’s go ahead and define a
Notice that
From here we derive three new quantities
It is a straight-forward to show that for
These values are useless independently, but because of linearity in matrix multiplication, we can use
To make the RHS
To make the RHS
Because
That is, use
The algorithm continues to replace the
C Implementation
The following implementation follows the naming of quantities in the above working. l
, u
and x
are first initialised to the quantities u
and x
get a new value appended to the end for each iteration. l
gets a new element pre-pended per iteration which removes the need for a reverse-access read.
/* Solves Tx = y for x where x and y are vectors and T is Toeplitz.
*
* `ord` is the number of rows and columns in the matrix.
* `t` contains the toeplitz diagonal coefficients from the bottom left to the
* top right (it must contain 2*ord-1 elements). t[ord-1] must be non-zero.
* `y` contains the desired output (it must contain ord elements).
* `x` will contain the computed solution. */
static void levindurb(double *x, const double *t, const double *y, unsigned ord) {
const int diag_offset = ord-1;
double u[16];
double l[16];
double tmp;
int i, j;
(ord < 16);
assert
= 1.0 / (t[diag_offset] * t[diag_offset] - t[diag_offset+1] * t[diag_offset-1]);
tmp [0] = tmp * t[diag_offset];
u[1] = tmp * -t[diag_offset-1];
u[ord-2] = tmp * -t[diag_offset+1];
l[ord-1] = u[0];
l[0] = y[0] * u[0] + y[1] * l[ord-2];
x[1] = y[0] * u[1] + y[1] * l[ord-1];
xfor (i = 2; i < ord; i++) {
double alpha, beta, r, scale;
double *p_l = &(l[ord-i]);
for (j = 0, alpha = 0.0, beta = 0.0, r = 0.0; j < i; j++) {
+= x[j] * t[diag_offset-i+j];
r += u[j] * t[diag_offset-i+j];
alpha += p_l[j] * t[diag_offset+1+j];
beta }
= 1.0 / (1.0 - beta * alpha);
scale for (j = 0; j <= i; j++) {
double ut = (j < i) ? u[j] : 0.0;
double lt = (j > 0) ? p_l[j-1] : 0.0;
double xt = (j < i) ? x[j] : 0.0;
[j] = scale * (ut - alpha * lt);
u[j-1] = scale * (lt - beta * ut);
p_l[j] = xt + (y[i] - r) * p_l[j-1];
x}
}
}