% 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.80,0.85);
    %[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)
    subplot(2,1,1);plot(wx/pi,(Wx/max(Wx))); title('input PSD'); grid
    subplot(2,1,2);plot(wa/pi,(H_2/max(H_2)),'r'); title ('parmetrized output PSD');grid
    pause;
    subplot(2,1,1);plot(wx/pi,log10((Wx/max(Wx)))); title('input PSD'); grid
    subplot(2,1,2);plot(wa/pi,log10((H_2/max(H_2)),'r')); title ('parmetrized output PSD');grid