% DSP Fall 2019
% LMS algorithm demo 
clear all; close all; hold off

% Channel system order
sysorder = 5;

% Number of system points
N = 1000;

x     = randn(N,1);
n1    = randn(N,1);
[b,a] = butter(2,0.25);

% Create a discrete-time transfer function model
Gz = tf(b,a,-1); 

% This function is submitted to make inverse Z-transform
% The first sysorder weight value
% h = ldiv(b,a,sysorder)';
% if you use ldiv this will give h :filter weights to be
h = [0.0976; 0.2873; 0.3360; 0.2210; 0.0964];
% Simulate time response of dynamic system to arbitrary inputs
y1 = lsim(Gz,x);

% Add noise
n1 = n1 * std(y1)/(10*std(n1));
d1 = y1 + n1;

total_L = size(d1,1);

% Begin the algorithm
w1 = zeros(sysorder,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set mu1, mu2
mu1 = 0.005;
mu2 = 0.15;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for n = sysorder : N 
	u    = x(n:-1:n-sysorder+1);
    y1(n) = w1' * u;
    e1(n) = d1(n) - y1(n);
	w1 = w1 + mu1 * u * e1(n);
end 

% Check results
for n1 = (N+1):total_L
	u     = x(n1:-1:n1-sysorder+1);
    y1(n1)= w1' * u;
    e1(n1) = d1(n1) - y1(n1);
end 

% Simulate time response of dynamic system to arbitrary inputs
y2 = lsim(Gz,x);
n2 = randn(N,1);

% Add noise
n2 = n2 * std(y2)/(10*std(n2));
d2 = y2 + n2;
w2 = zeros(sysorder,1);

for n = sysorder : N 
	u      = x(n:-1:n-sysorder+1);
    y2(n) = w2' * u;
    e2(n) = d2(n) - y2(n);
	w2     = w2 + mu2 * u * e2(n);
end 

% Check results
for n = N+1 : total_L
	u    = x(n:-1:n-sysorder+1) ;
    y2(n) = w2' * u ;
    e2(n) = d2(n) - y2(n) ;
end 

figure(1)
subplot(1,2,1)
hold on
plot(d1)
plot(y1,'r');
title_1 = strcat({'System Output, mu_1 = '}, num2str(mu1));
title(title_1);
xlabel('Samples')
ylabel('True and estimated output')
legend('Desired Output','Estimated Output')

subplot(1,2,2)
hold on
plot(d2)
plot(y2,'r');
title_2 = strcat({'System Output, mu_2 = '}, num2str(mu2));
title(title_2);
xlabel('Samples')
ylabel('Desired and estimated output')
legend('Desired Output','Estimated Output')

figure(2);
subplot(1,2,1)
plot(abs(e1))
title_3 = strcat({'Error Curve, mu_1 = '}, num2str(mu1));
title(title_3) ;
xlabel('Samples')
ylabel('Error value')
ylim([0 1.2])

subplot(1,2,2)
plot(abs(e2))
title_4 = strcat({'Error Curve, mu_2 = '}, num2str(mu2));
title(title_4)
xlabel('Samples')
ylabel('Error value')
ylim([0 1.2])


figure(3)
p1 = plot(h, 'k*'); p1.LineWidth = 10;
hold on
p2 = plot(w1, 'r+'); p2.LineWidth = 3;
p3 = plot(w2, 'b+'); p3.LineWidth = 3;
title('Comparison of the actual weights and the estimated weights');
legend('Actual weights', sprintf('Estimated Weights, mu_1 = %g', mu1), ...
    sprintf('Estimated Weights, mu_2 = %g', mu2))