# calculating gain with levinson-durbin.

Started by March 17, 2006
```Hi,

I'm using the levinson-durbin algorithm to calculate an all-pole filter.
This all-pole filter is used to generate background noise to overcome
silence when the telephone line is muted.

I believe a description of the all-pole filter is : H(z) = G / (a0+a1*z-1 +
a2*z-2....).

Calculating the denomenator in Matlab is not a problem, using the following
line of matlab code:

[ap e] = levinson(autocorrelation matrix, nr of poles);

ap is the denominator, e is the power prediction error.

Can someone tell me how to calculate the gain G of the transfer function
using the power prediction error?. I can not sort this out.

With kind regards,

Johan Kleuskens
The Netherlands

```
```Levinson assumes that the gain is 1. Check doc levinson in Matlab.

```
```GAIN=SQRT(E)

-------------
"Johan Kleuskens" <j.kleuskens@opentsp.com> skrev i en meddelelse
news:441a8d69\$0\$11070\$e4fe514c@news.xs4all.nl...
> Hi,
>
> I'm using the levinson-durbin algorithm to calculate an all-pole filter.
> This all-pole filter is used to generate background noise to overcome
> silence when the telephone line is muted.
>
> I believe a description of the all-pole filter is : H(z) = G / (a0+a1*z-1
> + a2*z-2....).
>
> Calculating the denomenator in Matlab is not a problem, using the
> following line of matlab code:
>
>    [ap e] = levinson(autocorrelation matrix, nr of poles);
>
> ap is the denominator, e is the power prediction error.
>
> Can someone tell me how to calculate the gain G of the transfer function
> using the power prediction error?. I can not sort this out.
>
> With kind regards,
>
> Johan Kleuskens
> The Netherlands
>

```
```I found that answer also in the description of the levinson in simulink.
However simulation results show that this in not correct.

news:dvetok\$2i93\$1@newsbin.cybercity.dk...
> GAIN=SQRT(E)
>
> -------------
> "Johan Kleuskens" <j.kleuskens@opentsp.com> skrev i en meddelelse
> news:441a8d69\$0\$11070\$e4fe514c@news.xs4all.nl...
>> Hi,
>>
>> I'm using the levinson-durbin algorithm to calculate an all-pole filter.
>> This all-pole filter is used to generate background noise to overcome
>> silence when the telephone line is muted.
>>
>> I believe a description of the all-pole filter is : H(z) = G / (a0+a1*z-1
>> + a2*z-2....).
>>
>> Calculating the denomenator in Matlab is not a problem, using the
>> following line of matlab code:
>>
>>    [ap e] = levinson(autocorrelation matrix, nr of poles);
>>
>> ap is the denominator, e is the power prediction error.
>>
>> Can someone tell me how to calculate the gain G of the transfer function
>> using the power prediction error?. I can not sort this out.
>>
>> With kind regards,
>>
>> Johan Kleuskens
>> The Netherlands
>>
>
>

```
```>I found that answer also in the description of the levinson in simulink.
>However simulation results show that this in not correct.

if you explain your simulation results (how you obtained them) then I might
be able to help

```
```Hello Jack,

What i did is the following:

f=firlpnorm(10,[0 0.5 0.6 1],[0 0.5 0.6 1],[1 0 0 0]); %generate the filter
to be modeled

fvtool(f,1); % show filter

x=(rand(500000,1)-0.5)*2; %generate white noise

xf=filter(f,1,x); %filtered noise

ac=xcorr(xf(1:1024)); %generate autocorrelation matrix

[ap e] = levinson(ac(1024:1024+10),10); %model filter

fvtool(sqrt(e),ap); %show filter

fvtool shows that the filter shapes match but the gain does not. The gain of
the modeled filter is about 27 dB. When changing the number of zeros in the
fir filter f to from 10 to 30 zeros, but the gain of the modeled filter
changes to 25dB. One should expect the gain would stay the same.

f=firlpnorm(30,[0 0.5 0.6 1],[0 0.5 0.6 1],[1 0 0 0]); %generate the filter
to be modeled

fvtool(f,1); % show filter

x=(rand(500000,1)-0.5)*2; %generate white noise

xf=filter(f,1,x); %filtered noise

ac=xcorr(xf(1:1024)); %generate autocorrelation matrix

[ap e] = levinson(ac(1024:1024+10),10); %model filter

fvtool(sqrt(e),ap); %show filter

Johan Kleuskens

news:dvh611\$joj\$1@newsbin.cybercity.dk...
> >I found that answer also in the description of the levinson in simulink.
> >However simulation results show that this in not correct.
>
> if you explain your simulation results (how you obtained them) then I
> might be able to help
>

```
```Here is my guess of what you are doing wrong:

1) You are using FIR filter coefficients
2) You are filtering random noise with a FIR filter

LPC analyzes data which is assumed to be the output
of an IIR filter.

Try this:

clc
close all
clear

% Number of points
N=100;

% IIR denominator order
P=10;

% Design Pth order IIR filter
[num,den] = iirlpnormc(0,P,[0 .15 .4 .5 1],[0 .4 .5 1],[1 1.6 1 0 0],[1 1 1
10 10])

% Random Noise Input
x=randn(1,N);

% Filter random noise
y=filter(num,den,x);

% Find filter impulse response
y_impulse=filter(num,den,[1 zeros(1,N-1)]);

% Calc autocorrelation (random noise scenario)
ac=xcorr(y);
ac=ac(N:N+P);

% Estimate filter coefficients
[ap e] = levinson(ac,P)

% Compare spectra
Syy=abs(fft(y)).^2;
Syy_predicted=abs(freqz(sqrt(e),ap,N,'whole')).^2;

plot(Syy)
hold
plot(Syy_predicted,'r--')
hold

% Calc autocorrelation (impulse response scenario)
ac=xcorr(y_impulse);
ac=ac(N:N+P);

% Estimate
[ap e] = levinson(ac,P)

% Compare
Syy_impulse=abs(fft(y_impulse)).^2;
Syy_predicted=abs(freqz(sqrt(e),ap,N,'whole')).^2;

figure
plot(Syy_impulse)
hold
plot(Syy_predicted,'r--')
hold

```
```I expanded your code a little bit, and i think now i know what is really
wrong:

I calculated the power ratio of the output signals of the reference and
predicted filters. In the "random noise section" the power ratio is around
88, meaning that the power of the output signal of the predicted filter is
much more than it should be.

In the "impuls respons section" the ratio is 1, as it should be.

Doesn't this mean that gain calculation works best if the auto correlation
is taken from an impuls respons, and when using noise for auto correlation a
very (probably very )inaccurate result is achieved?

By the way, i tried this also for a fir-filter, and the conclusion is the
same, as expected.

The modified code is listed below

Clc

close all

clear

% Number of points

N=100;

% IIR denominator order

P=10;

% Design Pth order IIR filter

[num,den] = iirlpnormc(0,P,[0 .15 .4 .5 1],[0 .4 .5 1],[1 1.6 1 0 0],[1 1 1
10 10]);

% Random Noise Input

x=randn(1,N);

% Filter random noise

y=filter(num,den,x);

% Find filter impulse response

y_impulse=filter(num,den,[1 zeros(1,N-1)]);

% Calc autocorrelation (random noise scenario)

ac=xcorr(y);

ac=ac(N:N+P);

% Estimate filter coefficients

[ap e] = levinson(ac,P);

% Compare spectra

Syy=abs(fft(y)).^2;

Syy_predicted=abs(freqz(sqrt(e),ap,N,'whole')).^2;

y_predicted = filter(sqrt(e),ap,x);

% Calculate power of output signals of reference and predicted filter

pow_ref = y * y';

pow_predicted = y_predicted * y_predicted';

pow_predicted/pow_ref

figure(1)

plot(Syy)

hold

plot(Syy_predicted,'r--')

hold

% Calc autocorrelation (impulse response scenario)

ac=xcorr(y_impulse);

ac=ac(N:N+P);

% Estimate

[ap e] = levinson(ac,P);

% Compare

Syy_impulse=abs(fft(y_impulse)).^2;

Syy_predicted=abs(freqz(sqrt(e),ap,N,'whole')).^2;

y_predicted = filter(sqrt(e),ap,x);

% Calculate power of output signals of reference and predicted filter

pow_ref = y * y';

pow_predicted = y_predicted * y_predicted';

pow_predicted/pow_ref

figure(2)

plot(Syy_impulse)

hold

plot(Syy_predicted,'r--')

hold

When using random noise
news:dvm6lp\$2l4n\$1@newsbin.cybercity.dk...
> Here is my guess of what you are doing wrong:
>
> 1) You are using FIR filter coefficients
> 2) You are filtering random noise with a FIR filter
>
> LPC analyzes data which is assumed to be the output
> of an IIR filter.
>
> Try this:
>
> clc
> close all
> clear
>
> % Number of points
> N=100;
>
> % IIR denominator order
> P=10;
>
> % Design Pth order IIR filter
> [num,den] = iirlpnormc(0,P,[0 .15 .4 .5 1],[0 .4 .5 1],[1 1.6 1 0 0],[1 1
> 1 10 10])
>
> % Random Noise Input
> x=randn(1,N);
>
> % Filter random noise
> y=filter(num,den,x);
>
> % Find filter impulse response
> y_impulse=filter(num,den,[1 zeros(1,N-1)]);
>
> % Calc autocorrelation (random noise scenario)
> ac=xcorr(y);
> ac=ac(N:N+P);
>
> % Estimate filter coefficients
> [ap e] = levinson(ac,P)
>
> % Compare spectra
> Syy=abs(fft(y)).^2;
> Syy_predicted=abs(freqz(sqrt(e),ap,N,'whole')).^2;
>
> plot(Syy)
> hold
> plot(Syy_predicted,'r--')
> hold
>
> % Calc autocorrelation (impulse response scenario)
> ac=xcorr(y_impulse);
> ac=ac(N:N+P);
>
> % Estimate
> [ap e] = levinson(ac,P)
>
> % Compare
> Syy_impulse=abs(fft(y_impulse)).^2;
> Syy_predicted=abs(freqz(sqrt(e),ap,N,'whole')).^2;
>
> figure
> plot(Syy_impulse)
> hold
> plot(Syy_predicted,'r--')
> hold
>
>
>
>

```