% ar_model.m
%
%  �ƒ�@AR_model identification  �„
%
Fs = 200;   T = 1;   N = Fs*T;
t = (0:1/Fs:T-1/Fs);
rn = 0.02*randn(size(t));
%
% ----------------------signal---------------------------
x = sin(2*pi*2*t) + sin(2*pi*10*t) + sin(2*pi*25*t) + rn;
% -------------FFT spectrum (Welch's method)-------------
[Pxx,f] = psd(x,256,Fs,hanning(256),128,'none'); 
% -------------------------------------------------------
%
%  AR model estimation 
%
M = 4;                                  % order of AR model
[a,e] = aryule(x,M);                    % Yule-Walker method
% [a,e] = arburg(x,M);                  % Burg method
Pxx_ar = e*abs(freqz(1,a,N/2,Fs)).^2;   % spectrum (AR model)
f_ar = (0:1/T:Fs/2-1/T);
%
subplot(311);plot(f,10*log10(Pxx),'-b',f_ar,10*log10(Pxx_ar),':r');
legend('Pxx (FFT spectrum)','Pxx (AR spectrum)')
ylabel('Magnitude (dB)'); xlabel('Frequency (Hz)');
%
%  prediction by the AR model
%
y = zeros(1,N);
a(1) = [];
for I=M+1:N
   for J=1:M
      y(I)=y(I)-a(J)*x(I-J);
   end
end
% ------------------------------
xest = y(M+1:N);
t(1:M) = [];
x(1:M) = [];
subplot(312);plot(t,x,'-b',t,xest,':r');ylabel('Amplitude')
axis([0,T,-3,3]);legend('x (obserbed)','x (predicted)')
subplot(313);plot(t,x-xest);
axis([0,T,-0.5,0.5]);legend('x-xest')
ylabel('Prediction Error'); xlabel('Time (sec)');
% --------------------------------------------------------
%                                  2001.7.09  by K.Tsukada