Levinson-Durbin Algorithm and Signal Modeling

Afshin Sepehri

Department of Electrical and Computer Engineering

University of Maryland

Email: afshin@eng.umd.edu

Introduction

The Levinson-Durbin algorithm uses the autocorrelation method to estimate the linear prediction parameters for a segment of a random signal. This project is about implementing the Levinson-Durbin algorithm and modeling random signals using it.

1.1. Implementation of the Levinson-Durbin Algorithm

In this part, an optimal all-pole model for a process is designed provided the autocorrelation values. The Levinson-Durbin was implemented according to the algorithm depicted in Table 5.1 in [1]. This is an iterative method, which is executed for p number of times where p is the order of the desired filter passed to the program. The outputs of this function are a_p, the optimal all-pole parameter vector, b[0], the numerator, e, a vector of squared errors and g, a vector of reflection coefficients. The implemented function in MATLAB follows:

function [ap, b0, e, g] = LevinsonDurbin(r, p);

ap = 0;

g = [];

aj(1) = 1;

ej = r(1);

e = [ej];

for j=1:p,

aj1 = zeros(j+1, 1);

aj1(1) = 1;

gammaj = r(j+1);

for i=2:j,

gammaj = gammaj + aj(i)*r(j-i+2);

end

lambdaj1 = -gammaj/ej;

g = [g ; lambdaj1];

for i=2:j,

aj1(i) = aj(i)+lambdaj1*(aj(j-i+2)');

end

aj1(j+1) = lambdaj1;

ej1 = ej*(1-abs(lambdaj1)^2);

e = [e ; ej1];

aj = aj1;

ap = aj1;

ej = ej1;

end

b0 = sqrt(ej1);

1.2. Debugging The Algorithm

To verify the written function, a sample process model is provided. Now by evaluating the autocorrelation and calling the above function, we try to find the same model. The spectrum of the process is:

P_x[z] = 4Q(z)Q^*(1/z^*)

Where

Q(z) = 1/(1-0.5z^-1)(1+0.25z^-1)

So,

Q(z) = (2/3)/(1-0.5z^-1) + (1/3)/(1+0.25z^-1)

By taking inverse Z transform we get:

q(n) = (2/3)0.5ⁿu(n) + (1/3)(-0.25)ⁿu(n)

and since autocorrelation r_x(n) is inverse Z transform of P_x[z], then:

r_x(n) = 4q(n)*q(-n)

That is the convolution of two sequences q(n) and q(-n). Now, by evaluating r_x(n) for 0£ n £N, we can find the required autocorrelation sequence for Levinson-Durbin algorithm implemented in previous section. By increasing N, we can find more and more accurate autocorrelation and therefore finding coefficients more accurately. However, in this case, since q(n) is a second order all-pole model, we can find exact position of the poles just by using a relatively small N. Table 1 shows how the position of the poles change by N. Note that, in this example, we know the correct order of the model which is p=2. So, by choosing the right order and high enough N, we should be able to find the right position of the poles, which is shown in the table.

N	First Pole (z = 0.5)	Second Pole (z = -0.25)
2	0.4674	-0.2254
3	0.4889	-0.2398
4	0.4977	-0.2480
5	0.4994	-0.2494
6	0.4998	-0.2499
7	0.5000	-0.2500

Table 1: Position of the model poles found using different N’s

1.3. Signal Modeling

In this section, we have an unknown random process to model using Levinson-Durbin method. To find the model, first we need to find the autocorrelation function as accurate as possible. To find this autocorrelation, we generate M different sample paths of the process:

x₀(n), x₁(n), … , x_M-1(n)

For each x_i(n) we find r_xi(n) as

r_xi(k) = S_n=k^N x_i(n) x_i(n - k) k = 0, …., N

and then by averaging:

r_x(k) = (1/M) S_i=0^M-1 r_xi(k)

we can find the autocorrelation pattern. Obviously, by increasing the two parameters N and M we can get more accurate autocorrelation values. Figure 1 shows power spectrum of x (Fourier transform of autocorrelation). WE can see sensitivity of the autocorrelation with respect to N and M. By changing N and M from 100 to 1000 we can get a better view of the spectrum.

Figure 1: Power spectrum of random process x(n) using different values of N and M

Top left: N=100 and M=100 Top right: N=1000 and M=100

Bottom left: N=100 and M=1000 Bottom right: N=1000 and M=1000

The larger view of power spectrum can be seen in figure 2.

Figure 2: Power spectrum of BlackBoxA using N=1000 and M=2000

Now we pass the autocorrelation to the implemented Levinson-Durbin algorithm of section 1.1 to get the parameters of the model. The important problem here is the order of the model we are trying to find. WE know that by increasing the order, the error would be decrease. However it is not a good idea to have very high-order. So we try to find the optimal order by looking at the modeling error. One advantage of Levinson-durbin method is that we can get e as a vector of errors of modeling for different orders. These errors are shown in figure 3 for orders up to 25.

Figure 3: Errors of modeling with respect to the order of the all-pole model

As it can be seen clearly, error is decreased drastically by increasing the order up to order 4 and afterwards, no real change is made. Therefore, we can pick order four as the order of the optimal model. (i.e. p=4)

By applying this, we find the model denominator coefficients (a_p) as:

a_p = [1.0000 0.0998 0.2618 0.0568 -0.1688] ^T

and the numerator coefficient is

b₀ = 31.5882

the reflection Coefficients G as returned in vector g are:

g = [0.0682 0.3082 0.0759 -0.1688] ^T

We can find the model poles as roots of the denominator polynomial:

Poles = [-0.6006 0.0004 - 0.7498j 0.0004 + 0.7498j 0.4999]

So, by looking at the found poles, I personally guess that the real poles used in the model were

Guessed Poles = [ -0.6 -0.75j 0.75j 0.5 ]

Which were found using our method with a reasonable accuracy. However, these were guessed by the assumption that the original model is an all-pole model. If it were an ARMA or MA model, so the poles would not be justified.