%
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; clear; clear functions; 

phi_1 = 1.2728;
phi_2 = -0.81;

addpath('../../../Poularikas/AdaptiveFilteringCode');

% generate samples from this timeseries: 
% 
N = 500;

% the random driving input noise 
% 
v  = randn(1,N+1);                              

% our input signal is given by x(n) ... 
% ... we want to predict this signal from the previous two values:
% 
a_coefs = [ 1, -phi_1, -phi_2 ];
b_coefs = [ 1 ]; 
x = filter( b_coefs, a_coefs, v ); 
  
% check that this input has the correct autocorrelation matrix:
s_index = round( 0.1*N );
[C,lags] = xcorr( x(s_index:N), x(s_index:N), 5, 'coeff' );
r_x_0 = C( find(lags==0 ) );
r_x_1 = C( find(lags==+1) ); 
R_x = [ r_x_0, r_x_1; r_x_1, r_x_0 ]
r_x_tr = trace(R_x);
eta_bound = 1/r_x_tr;
eta = eta_bound/20.;

dn = x(2:end); % predict one step ahead using the previous two x values; dn(1) = x(3) = x_2; dn(2) = x(3) = x_2, ... dn(end) = x(end) = x_N 
x  = x(1:end-1); 

if 0
  % the LMS algorithm:
  [w,y,e,J,w1] = aalms1( x, dn, eta, 2 ); % the version with a *1* returns an iteration history of the coefficients weights 

  ph = plot( w1, '-o' );
  legend( ph, {'phi_1 estimate', 'phi_2 estimate'}, 'location', 'best' );
  xlabel('iteration number'); ylabel('weight value');
  hold on; 
  plot( phi_1 * ones(1,length(w1)), '-b' );
  plot( phi_2 * ones(1,length(w1)), '-g' );
  %saveas(gcf,'../../WriteUp/Graphics/Chapter8/prob_8_1_LMS_convergence.eps', 'epsc');
else
  % the normalize LMS algorithm: 
  c = 1.;
  [w,y,e,J] = aanlms( x, dn, eta, 2, c );
end