a = poly([0.2 0.4+0.87i 0.4-0.87i -0.1+0.1i -0.1-0.1i -0.89 0.88]); 
na = length(a)-1;
b = [0.0 2.4 4.4 8.4 7.4];
nb = length(b)-1;
%
% stima dei parametri da una coppia ingresso/uscita qualunque:
N = 64;
for stde = [0 1e-16 1e-12 1e-10 1e-9 1e-8 1e-6 1e-5 5e-5 1e-4 2e-4]
  figure(20); clf, hold on,
  figure(21); clf, hold on, plot(roots(a),'bo'); title('radici pol. "a"'); pause
  if 1
    omega_nu = (2*pi/N) * 0.1*(N/2);
    u = sin(omega_nu * (1:N));
    if 1
      disp(['ingresso sinusoidale con rumore (std=' num2str(stde) ') in uscita:']);
      y = filter(b, a, u) + stde*rand(1,N) ;  % questa routine risolve l'equazione alle differenze vista alla "ESE_1"
    else
    end
    figure(2); clf, hold on, plot(1:N,u(1:N),'g-', 1:N,y(1:N),'b-'); 
    figure(44); clf, hold on,
    if stde == 0
      for na_est=3:na+4
        A = [hankel([zeros(1,na_est-1), u(1:N-na_est+1)], u(N-na_est+1:N)) hankel([zeros(1,na_est-1), y(1:N-na_est+1)], y(N-na_est+1:N-1))];
        [U, S, V] = svd(A);
        figure(44); semilogy(log10(diag(S)),'bo'); 
        %disp(['valori singolari per na_est=' num2str(na_est) '  e rumore (std=' num2str(stde) '):']), diag(S),
        pause,
        figure(44); semilogy(log10(diag(S)),'go'); 
      end;
    end;
    A = [hankel([zeros((0)+1,(na-1)+1), u((0)+1:(N-1-na)+1)], u((N-1-na)+1:(N-1)+1)) ...
         hankel([zeros((0)+1,(na-1)+1), y((0)+1:(N-1-na)+1)], y((N-1-na)+1:(N-1-1)+1))];
    [U,S,V] = svd(A);   if stde>0  figure(20); semilogy(log10(diag(S)),'r*'); end;
    d = y((0)+1:(N-1)+1)';
    [Q, R] = qr(A);
    coeff_est = R \ (Q'*d);
    a_est = [1 -coeff_est(end:-1:na+2)'];
    disp(['a = ' mat2str(a,4)]);
    disp(['a_est = ' mat2str(a_est,4)]);
    %sum(abs(a)-abs([1 -coeff_est(end:-1:na+2)']))
    b_est = coeff_est(na+1:-1:na-nb+1)';
    disp(['b = ' mat2str(b,4)]);
    disp(['b_est = ' mat2str(b_est,4)]);
    %sum(abs(b)-abs(coeff_est(na+1:-1:na-nb+1)'))
    %disp('coeff_est :'), coeff_est
    y_est = filter(b_est, a_est, u) ;  
    figure(2); plot(1:N,y_est(1:N),'r--'); 
    figure(21); plot(roots(a_est),'r*'); pause
    pause
    %
    disp(' ');
    disp(' ');
    disp(' ');
  end
  disp(['ingresso impulsivo con rumore (std=' num2str(stde) ') in uscita:']);
  u=zeros(1,N);  u(1)=1;
  y = filter(b, a, u)  + stde*rand(1,N) ;  % questa routine risolve l'equazione alle differenze vista alla "ESE_1"
  figure(22); clf, hold on, plot(1:N,u(1:N),'g-', 1:N,y(1:N),'b-'); pause 
  if stde >= 0
    for na_est=3:na+4
      A = [hankel([zeros(1,na_est-1), u(1:N-na_est+1)], u(N-na_est+1:N)) hankel([zeros(1,na_est-1), y(1:N-na_est+1)], y(N-na_est+1:N-1))];
      [U, S, V] = svd(A);
      figure(44); semilogy(log10(diag(S)),'ro'); 
      disp(['valori singolari per na_est=' num2str(na_est) '  e rumore (std=' num2str(stde) '):']), diag(S),
      title(['std(err) = ' num2str(stde)]);
      pause,
      figure(44); semilogy(log10(diag(S)),'mo'); 
    end;
  end;
  A = [hankel([zeros((0)+1,(na-1)+1), u((0)+1:(N-1-na)+1)], u((N-1-na)+1:(N-1)+1)) ...
       hankel([zeros((0)+1,(na-1)+1), y((0)+1:(N-1-na)+1)], y((N-1-na)+1:(N-1-1)+1))];
  [U,S,V] = svd(A);   if stde>0 figure(20); semilogy(log10(diag(S)),'bo'); title(['std(err) = ' num2str(stde)]); pause, end;
  d = y';
  [Q, R] = qr(A);
  coeff_est = R \ (Q'*d);
  a_est = [1 -coeff_est(end:-1:na+2)'];
  disp(['a = ' mat2str(a,4)]);
  disp(['a_est = ' mat2str(a_est,4)]);
  %sum(abs(a)-abs([1 -coeff_est(end:-1:na+2)']))
  b_est = coeff_est(na+1:-1:na-nb+1)';
  disp(['b = ' mat2str(b,4)]);
  disp(['b_est = ' mat2str(b_est,4)]);
  %sum(abs(b)-abs(coeff_est(na+1:-1:na-nb+1)'))
  y_est = filter(b_est, a_est, u) ;  
  figure(22); plot(1:N,y_est(1:N),'r--'); pause 
  figure(21);  title(['std(err) = ' num2str(stde)]); plot(roots(a_est),'b+'); pause,
end;
%
% 



