%  adsid.m
%            7 March 2005
% Adaptive System Identification

clear, 
figure(1), clf

% Input signal x(t)

f = 0.8 ;      %  Hz
ts = 0.005 ; % 5 msec -- sampling time
N1 = 800 ; N2 = 400 ; N = N1 + N2 ;

t1 = (0:N1-1)*ts ;  % 0 to 4 sec
t2 = (N1:N-1)*ts ;  % 4 to 6 sec
t = [t1 t2] ;       % 0 to 6 sec

xt = sin(3*t.*sin(2*pi*f*t)) ;
% figure(1)
% plot(t, xt), grid, title('Input Signal, x(t)'), xlabel('time sec')

p = 3 ; % Dimensionality of the system

b1 = [ 1  -0.6 0.4] ; % unknown system paremeters during  t1
b2 = [0.9 -0.5 0.7] ; % unknown system psremeters during  t2

 [d1, stt] = filter(b1, 1, xt(1:N1) ) ;
 d2 =  filter(b2, 1, xt(N1+1:N), stt) ;
 dd = [d1 d2] ; %  target signal

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

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

% Alternatively, we could calculated D as

d = [b1*X(:,1:N1) b2*X(:,N1+1:N)] ;
% plot(t, d-dd), grid,

% plot(t, d), grid, title('Output Signal from the System, d(n)'), ... 
% xlabel('time sec')

 y  = zeros(size(d)) ;   % memory allocation for y
eps = zeros(size(d)) ;   % memory allocation for eps

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

for n = 1:N  %  learning loop

 y(n) = w*X(:,n) ;      % predicted output signal
 eps(n) = d(n) - y(n) ; % error signal
 w = w + eta*eps(n)*X(:,n)' ;
 if n == N1-1, w1 = w ; end
end

figure(1)
myfigA(1,11,10)
subplot(2,1,1)
plot(t, xt), grid on, 
title('Input Signal, x(t)'), % xlabel('time sec')
subplot(2,1,2)
plot(t, d, 'b', t, y, '-r'), grid on 
axis([0 6 -1.2 1.2])
title('Target d(t) and predicted y(t) signals'), 
xlabel('time [sec]')
drawnow
  print -f1 -deps cAdsid2

figure(2)
myfigA(2,11,10)
plot(t, eps), grid on
title(['prediction error for eta = ', num2str(eta)])
xlabel('time [sec]')
drawnow
  print -f2 -deps cAdsid3

[b1; w1]
[b2;  w]