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.

	function armcov(x,L)
	N = length(x)
	phi = zeros(L+1,L+1)
	for k = 0:L, i=0:L, n=0:N-L-1
	phi[k+1,i+1] += x[n+1+k]*x[n+1+i] + x[n+1+L-k]*x[n+1+L-i]
	end
	phi = phi./(2*(N-L))
	a = phi[2:end,2:end]\-phi[2:end,1]
	P = phi[1,1] + phi[1,2:end]*a
	[1,a],P
	end
	x =[1:10,10:-1:1]
	a,P = armcov(x,5)
	print(a)
	print(P)

view raw armcov.jl hosted with ❤ by GitHub