Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % E5.m PSD parametric estimation for stationary signals
- % -----------------------------------------------------------
- clear;
- T=2000; % number of the input signal samples
- K=50; % number of the auto-correlation samples
- N=50; % Levinson filter order
- %x=gencomp(T,1,T,1,T,1,T,0.1,0.0,0.0);
- %x=gencomp(T,1,T,1,T,1,T,0.1,0.2,0.0);
- %x=gencomp(T,1,T,1,T,1,T,0.1,0.2,0.3);
- %x=gencomp(T,1,T,1,T,1,T,0.1,0.5,0.9);
- x = gen(T,0.1,0.15,0.20,0.25);
- %[x,fp]=wavread('s1.wav'); % reading a speech sample s1.wav (sampling frequency: fp=12000 Hz)
- %T=length(x);
- %y=decimate(x,2); % decimation to fp=6000 Hz
- %wavwrite(y,fp/2,16,'filename.wav'); % writing a decimated speech sample s1.wav (fp=6000 Hz, quantization 16 bits)
- xcentr = cntr(x);
- xstand = stnd(xcentr);
- [cs,ccs,nccs]=korel(xstand,K);
- [rho,A,e] = lev(cs,N);
- for i=1:N a(i)=A(i,N); end
- for i=N+1:T a(i)=0; end
- [Ha,wa] = ampl(a);
- H2=Ha.*Ha;
- [Wx,wx] = wgmx(xstand);
- We=(Wx/max(Wx)).*(H2/max(H2));
- for i=1:length(H2) H_2(i)=1/H2(i); end % squared inverse of the Levinson filter magnitude characteristic
- y=conv(xstand,A(:,N));
- [Wy,wy] = wgmx(y);
- %subplot(1,1,1); plot(wa/pi,We/max(We)); title('wgm of output sig'); grid;
- %pause
- %subplot(1,1,1); plot(nccs,ccs); title('Input signal auto-correlation'); grid
- %pause
- %subplot(1,1,1); stem(A(:,N)); title('Lev. filter imp. resp.'); grid
- %pause
- %subplot(1,1,1); stem(rho(2:N)); title('Schur coeff.'); grid
- %pause
- %subplot(1,1,1); plot(e); title('Least-squares error'); grid;
- %pause
- subplot(3,1,1); plot(wx/pi,Wx/max(Wx)); title('Input signal PSD'); grid;
- %subplot(3,1,2); plot(wa/pi,Ha/max(Ha)); title('Lev. filter magnitude'); grid;
- %subplot(3,1,3); plot(wy/pi,Wy/max(Wy)); title('Output signal PSD'); grid;
- %pause
- %subplot(2,1,1); plot(wx/pi,log10(Wx/max(Wx))); title('Input signal PSD'); grid;
- %subplot(2,1,1); plot(wa/pi,log10(H_2/max(H_2)),'r');
- %subplot(2,1,2); plot(wa/pi,log10(H_2/max(H_2))); title('Parametric output signal PSD estimate'); grid;
- subplot(2,1,2); plot(wa/pi,H_2/max(H_2)); title('Parametric output signal PSD estimate'); grid;
- %subplot(1,1,1)
- %plot(wx/pi,(Wx/max(Wx)));
- %plot(wa/pi,(H_2/max(H_2)),'r');
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement