% solution to Examples Sheet - Spectrum Analysis - Question 17
% S. S. Singh, November 2008

clear all

% declaring the true AR parameters
a1=0;
a2=0.81;
b0=1;

N =1000; % no of data points

% now generate samples from the AR process
x=[ 0 0]';  % to store previous two values [x_{-1} x_{-2}]' which are assumed zero 
a=[a1 a2]';
X=[];
for i=1:(1000+N)
    xnew=-1*a'*x+b0*randn;
    X=[X;xnew]; % append sample to record
    x=[xnew x(1)]';  % update for iteration i+1
end


Y=X(end-N+1:end);   % these are now the generated 24 samples




%%%%%%%%%%%%%%%
% Estimate the autocorrelation 
[Rhat, lags] = xcorr(Y,2,'unbiased'); % check help file for functionxcorr
                                      % selecting unbiased estimate
                                      % selecting 2 means R_XX[-2],R_XX[-1],R_XX[0],R_XX[1],R_XX[2] will be computed
                                      % and stored in Rhat in the same order
                                    
                                      
%%%%%%%%%%%%%%%
% Estimate the AR coefficients using Yule Walker equations                                     
aHat=-1*inv([Rhat(3) Rhat(4);Rhat(4) Rhat(3) ])*[Rhat(4) Rhat(5)]'
b0Hat=sqrt([ Rhat(3) Rhat(4) Rhat(5)]*[1; aHat])



%%%%%%%%%%%%%%%
% Estimated PSD using the Periodogram
Periodo=fft(Y,1024);
Periodo=(1/N)*abs(Periodo).^2;


%%%%%%%%%%%%%%%
% Estimated PSD from estimated AR coefficients
w=2*pi*[0:1023]/1024;
clear j
Shat=abs(1+ aHat(1)*exp(-j*w)+aHat(2)*exp(-j*w*2)).^2;
Shat=b0Hat^2./Shat;



%%%%%%%%%%%%%%%
%true PSD
w=2*pi*[0:1023]/1024;
clear j
S=abs(1+ a(1)*exp(-j*w)+a(2)*exp(-j*w*2)).^2;
S=b0^2./S;

%%%%%%%%%%%%%%%
%plot
figure(101)
subplot(3,1,1)
plot(w, Periodo)
temp=axis;
title('Periodogram')
hold on
subplot(3,1,2)
plot(w,  Shat)
title('Parametric')
axis(temp);
hold on
subplot(3,1,3)
plot(w, S)
title('True PSD')
axis(temp);
hold on


