%  adlnc.m
%            7 March 2005
% Adaptive Noise Cancelation

clear, 
figure(1), clf, figure(2), clf

% The input signal u(t) is a frequency and 
%                  amplitude  modulated sinusoid

f  = 4e3 ;       %  signal frequency 
fm = 300 ;      % frequency modulation
fa  = 200 ;     % amplitude modulation

ts = 2e-5 ;  % sampling time
N =  400 ;   % number of sampling points
t = (0:N-1)*ts ;  % 0  to 10 msec

ut = (1+0.2*sin(2*pi*fa*t)).*sin(2*pi*f*(1+0.2*cos(2*pi*fm*t)).*t) ;

% The noise is 
xt = sawtooth(2*pi*1e3*t, 0.7) ;

% the filtered noise
b = [ 1  -0.6 -0.3]  ;
vt = filter(b, 1, xt) ;

% noisy signal
dt = ut+vt ;

figure(1), % myfigA(1, 11,10)
subplot(2,1,1)

plot(1e3*t, ut, 1e3*t, dt), grid on
axis([0 8 -1.3 1.3])
title('Input  u(t)  and noisy input signal  d(t)'), 

subplot(2,1,2)

plot(1e3*t, xt, 1e3*t, vt), grid on
title('Noise  x(t)  and coloured noise  v(t)'), 
xlabel('time -- msec')
drawnow
%  print -f1 -depsc2 cAdlnc2

p = 4 ; % dimensionality of the input space

% formation of the input matrix X of size  p by N

X = convmtx(xt, p) ; X = X(:, 1:N) ;

 y  = zeros(1,N) ;   % memory allocation for y
eps = zeros(1,N) ;   % memory allocation for uh = eps

eta = 0.02  ; % learning rate/gain
w = 2*(rand(1, p)  -0.5) ; % Initialisation of the weight vector

for c = 1:4

for n = 1:N  %  learning loop
 y(n) = w*X(:,n) ;      % predicted output signal
 eps(n) = dt(n) - y(n) ; % error signal
 w = w + eta*eps(n)*X(:,n)' ;
end

eta = 0.8*eta ;

end

figure(2), %  myfigA(2,11,10)
subplot(2,1,1)
plot(1e3*t, ut, 1e3*t, eps), grid on
title('Input signal  u(t)  and estimated signal  uh(t)')
xlabel('time -- msec')
subplot(2,1,2)
plot(1e3*t(p:N), ut(p:N)-eps(p:N)), grid on
title('estimation error'), xlabel('time --[msec]')
drawnow

% print -f2 -depsc2 cAdlnc3

% w = 0.7933   -0.0646   -0.7231    0.0714