% 
% For illustrative purposes this duplicates the MATLAB results on epage 124 from the book.
% 
% Written by:
% -- 
% John L. Weatherwax                2008-11-24
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

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

clear functions; 
clear; 
%clc; 
close all; drawnow; 

randn('seed',0); rand('seed',0); 

N  = 4000; 
n  = 0:(N-1); 

% some random Gaussian driving input noise ... 
% 
v_std = 1.0;  % <- add more or less noise ... 
v     = v_std*randn(1,N);    

% our desired signal s(n) is ( this is unknown to the algorithm ) ... 
% 
s = sin( 0.2*pi*n ); 

% the signal plus noise:
% 
spv = s + v; 

% our unknown system will be a MA(4) system and have these coefficients
% --we denote our system by x(n):
w_system = [1, 0.2, -0.5, +0.75]; 
x        = filter( w_system, 1, spv ); 

% our desired signal is the noised input s(n): 
sd = spv; 

% plot the output of our noise modified signal:
% 
fo=figure; hs=plot( n, spv, '-og', 'markersize', 6 ); hold on;
xlabel( 'n index' ); ylabel( 'observed values' ); title( 's(n)+v(n) with y(n)' ); 

% 
% for one learning rate and value of M: 
% 
M1 = 2; 
mu = 0.01; 
[w,y,e,J,w1] = aalms1(x,sd,mu,M1); 
fprintf('standard deviation of the error when M=%10d = %10.6f\n',M1,std(e)); 

figure(fo); hy1=plot( n, y, '-xr' ); % <- the predicted signal at each timestep ... should match spv(n) 

fe = figure; eh1 = plot( n, e, '-k' ); hold on; % <- the error as a function of n 

fw = figure; ws1 = plot( w1, '-' ); hold on; % <- the convergence of the weights ... 

fJ = figure; plot( n, J, '-x' ); 

% 
% for another learning rate and another value of M: 
% 
M2 = 6; 
mu = 0.001; 
[w,y,e,J,w2] = aalms1(x,sd,mu,M2); 
fprintf('standard deviation of the error when M=%10d = %10.6f\n',M2,std(e)); 

figure(fo); hy2=plot( n, y, '-xk' ); 
legend( [hs,hy1,hy2], {'signal plus noise', 'LMS prediction #1', 'LMS prediction #2'} ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter7/inv_system_w_lms_signal_plus_pred.eps', 'epsc' ); 

figure(fe); eh2 = plot( n, e, '-r' ); axis tight; 
xlabel( 'n' ); ylabel( 'prediction error' ); 
legend( [eh1,eh2], {['error under M=',num2str(M1)], ['error under M=',num2str(M2)]} ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter7/inv_system_w_lms_error.eps', 'epsc' ); 

figure(fw); ws2 = plot( w2, '-' ); hold on; 
xlabel( 'n' ); ylabel( 'weight convergence' ); title( 'convergence of the weights' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter7/inv_system_w_lms_weight_conv.eps', 'epsc' );