a = [1 0.97];
na = length(a)-1;
b = [1 2];
% condizioni iniziali :
u_m1=0;
y_m1=0;
y_m2=0;
% 
lambda = 0.94;  % 0.94
N=64;
omega_nu = (2*pi/N) * 0.3*(N/2);
figure(2); clf, hold on, title('ingresso (verde) uscita (blu)')
a_est = zeros(na,N);
a_hist = zeros(na,N);
w = randn(1,N);
A=[]; d=[];
for n=1:N
  % legge di variazione dei parametri:
  if n>20 a(2)=-0.97; end;
  a_hist(:,n) = a(2:end);
  % calcolo l'ingresso
  %  u(n) = cos(omega_nu * n);    plot([n-1 n],[u_m1 u(n)],'g-');
  if n==1  u(n)=1; else u(n)=0; end;
  % aggiorno l'uscita del sistema
  y(n) = 1/a(1) * ( -a(2)*y_m1 + b(1)*u(n) + b(2)*u_m1 );    
  if 1
    y(n) = y(n) + 0.3*w(n);   % 0.3
  end;
  plot([n-1 n],[y_m1 y(n)],'b-');
  % calcolo le stime
  %keyboard,
  if n > na+1  rigadati = y(n-na:n-1)'; end;
  if n > 2*na+1
    A = [A; rigadati];
    d = [d; -y(n)];
    [Q, R] = qr(A);
    a_est(:,n) = R \ (Q'*d);
  elseif n == 2*na+1
    A = [A; rigadati];
    d = [d; -y(n)];
    [Q, R] = qr(A);
    a_est(:,n) = R \ (Q'*d);
  elseif n > na+1
    A = [A; rigadati]; 
    d = [d; -y(n)];
  end;
  if 1 & n > na+2  A = lambda * A;   d = lambda * d;  end;
  % aggiorno gli elementi di ritardo (NB: da quello più lontano nel tempo !!!)
  u_m1 = u(n);
  y_m2 = y_m1;
  y_m1 = y(n);
end;
figure(1); %clf, 
           hold on,
           plot(a_hist,'b-');
           plot(a_est,'k-');
