%  Linear Prediction Example: AR(2)
%  Kalman Filter Example 
%  State Model set up: 
a_1 = 2*0.99; a_2 = -0.99^2;
sigw = 0.05; sig_eps = 0.5; N = 500;
Phi = [a_1 a_2 ; 1 0]; C = [1 0]; 
Q = [sigw 0; 0 0]; R = sig_eps^2; 
%  Generate Observation process
w = awn(zeros(1,N),[0 sigw]);
d = filter(1,[1,-a_1,-a_2],w);
obs = awn(d,[0 sig_eps]); 
[X,K,MMSE] = dkf(obs,Phi,Q,R,C);
% Plot State Estimates
plot([0:1:N-1],X(1,:),[0:1:N-1],d,'r--')
set(gca,'FontName','Helvetica','FontSize',16)
xlabel('TIME   SAMPLE')
ylabel('d[n] & DKF state esitmate x[n]')
title('Prediction Example: \sigma_w = 0.05, \sigma_\epsilon = 0.5')
legend(gca,'EST','ACT',4)
print -deps dkf_est.eps
% Plot Kalman Gain vector
Kal = K.'; kal1 = Kal(:,1); kal2 = Kal(:,2);
plot([0:1:20-1],kal1(1:20),[0:1:20-1],kal2(1:20),'r--')
set(gca,'FontName','Helvetica','FontSize',16)
xlabel('TIME   SAMPLE')
ylabel('Kalman Gain: K[n]')
title('Prediction Example: \sigma_w = 0.05, \sigma_\epsilon = 0.5')
legend(gca,'K[1]','K[2]',4)
print -deps dkf_kal.eps
% Plot the MMSE
plot([0:1:20-1],MMSE(1:20))
set(gca,'FontName','Helvetica','FontSize',16)
xlabel('TIME   SAMPLE')
ylabel('MMSE (dB)')
title('Learning Curve: \sigma_w = 0.05, \sigma_\epsilon = 0.5')
print -deps dkf_mmse.eps
