% a synthetic noise cancellation problem
clear all

T=1000; %creating input signal u(n) for T time steps
t=1:T;

s=rand(size(t));    % - uncorrupted speech signal
                    % - we just use white noise as speech
                    
d=s+0.5*sin(t*pi/2+0.5);    % - the speech is corrupted by noise
                            % - we use a sinusoid for the corrupting noise
                            % - we have just created our desired signal
                            % - this is what Mic 1 measures

u=10*sin(t*pi/2);       % - this is the measurement of Mic 2 
                            % - Mic 2 measures just noise
                            % - this noise is a sinusoid too but with different amplitude and delay
                            % - this noise is correlated with the noise
                            %   measured at Mic 1
                            
% initialization of the LMS
M=2;
h=zeros(M,T);
h(:,1)=zeros(M,1);  % just setting the filter h to all zeros
mu=0.0001;
% initialization
for i=1:M
    uv(1,i)=u(1,M-i+1); % - initialising vector u(n) 
                        % - remember the input signal is appended with its
                        % past to form this vector
end
Y=[];   % to save the output of the filter
% and off we go! The LMS is being iterated
for n=M+1:T-1
    uv(1,:)=[u(1,n) uv(1,1:M-1)]; % appending the new data point
    e(n)=d(n)-uv*h(:,n);   % computing the error
    h(:,n+1)=h(:,n)+mu*e(n)*uv'; % updating the filter h
    Y=[Y,uv*h(:,n)];    % saving the output of the filter to plot
end

%%%
%% You should plot the different signals to understand how things are
%% working
%% Here is an example of what you can plot

figure
plot(e-s(1:T-1))    % - e should be close to s if the corrupting sinusoid has beed removed
                    % - e and s get closer as the filter converges
                    
title('error between true and estimated speech')
                    

figure
subplot(2,1,1)
plot(s(800:end),'r')
hold on
plot(d(800:end))
legend('speech recorded','speech orig')
                    
subplot(2,1,2)
plot(s(800:end),'r')
hold on
plot(e(800:end),'k-')
legend('speech orig', 'speech recon')

figure
plot(h')
title('Tap weights')

figure
plot(d(M+1:end)-s(M+1:end));
hold on;
plot(Y(1:end),'r')
title('Mic 1 noise (blue), Filter output (red)')


                  