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
calculating gain with levinson-durbin.
Started by ●March 17, 2006
Reply by ●March 17, 20062006-03-17
Reply by ●March 17, 20062006-03-17
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 >
Reply by ●March 17, 20062006-03-17
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 >> > >
Reply by ●March 18, 20062006-03-18
>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
Reply by ●March 20, 20062006-03-20
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 >
Reply by ●March 20, 20062006-03-20
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
Reply by ●March 20, 20062006-03-20
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 > > > >