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.



Many thanks for your input.



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
"Jack" <jack@nospammers.please> wrote in message 
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
>
>
>
> 

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

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



This is why i think your answer is not correct.



Johan Kleuskens








"Jack" <jack@nospammers.please> wrote in message 
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
> 

>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 

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

"Jack" <jack@nospammers.please> schreef in bericht 
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
>>
>
> 

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
> 

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

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