## Tuesday, 4 February 2014

### LPCs using the Modified Covariance Method in Julia

This post will discuss the modified covariance method (a.k.a the forward-backward method) of computing LPCs. Other methods include the Autocorrelation method, the Covariance method and the Burg method.

Linear prediction coefficients (LPCs) are used a lot in signal processing, and the theory around it is covered in other places (see here for an intuitive explanation). In this post I will go over how to implement LPCs using the modified covariance method in Julia. Given a signal $$x$$, the job of LPCs is to predict the current sample i.e. $$x(n)$$ using a weighted sum of previous samples:

$x(n) = \sum_{i=1}^L a(i) x(n-i) + e(n)$

Ideally, $$e(n)$$, the error, will be zero but usually it is not. As a result we want values of $$a(i)$$ that make $$e(n)$$ as small as possible. The autocorrelation method resulted in the Yule-Walker equations, which have a nice structure (the resulting matrix is toeplitz). The modified covariance method results in a similar set of equations, but the resulting matrix is not toeplitz, so we can't use the Levinson-Durbin algorithm to solve it. This means we have to rely on good old matrix inversion. The covariance equations we have to solve look like this:

$\sum_{i=1}^L \phi_{i,k} a(i) = -\phi_{k,0} \qquad 1 \leq k \leq L$

where the $$\phi$$s are defined as:

$\phi_{i,k} = \dfrac{1}{2(N-L)} \sum_{n=0}^{N-L-1} x(n+k)x(n+i) + x(n+L-k)x(n+L-i)$

where N is the length of the signal we are working with, and L is the order of the model we want. The matrix equation that results looks like this:

$\left[\begin{array}{cccc} \phi_{1,1} & \phi_{2,1} & \ldots & \phi_{L,1} \\ \phi_{1,2} & \phi_{2,2} & \ldots & \phi_{L,2} \\ \vdots & \vdots & \ddots & \vdots\\ \phi_{1,L} & \phi_{2,L} & \ldots & \phi_{L,L} \end{array}\right] \left[\begin{array}{c} a(1) \\ a(2) \\ \vdots \\ a(L) \end{array}\right] = -\left[\begin{array}{c} \phi_{1,0} \\ \phi_{2,0} \\ \vdots \\ \phi_{L,0} \end{array}\right]$

This can also be written as:

$Pa = -p$

Note that the matrix $$P$$ is symmetric but not toeplitz. To solve it we use:

$a = P^{-1}({-p})$

Once we have the values of $$a$$ we are done. Julia code for this procedure follows.