abbiamo_control_system_toolbox=1;
clear i;
%
% stima dei parametri di un modello ARMA dalla risposta al campione unitario
%
%
% definiamo il modello ARMA in modo arbitrario:
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;
%
% calcolo la risposta al campione unitario:
N = 64;
u=zeros(1,N);  u(1)=1;
h = filter(b,a,u)';
figure(1); plot(0:length(h)-1,h);
freqresp(h,5);
disp('campione unitario:')
%
% ora stimo i coefficienti "a" e "b" a partire da "h":
%A = hankel(h(2:na), h(na:N-1))';
A = hankel(h(2:N-na), h(N-na:N-1));
d = - h(na+2:N);
[Q, R] = qr(A);
B = hankel([zeros(1,na), h(1)], h(1:na+1));
a_est = R \ (Q'*d);
a_est = [a_est; 1];
disp(['a = ' mat2str(a,4)]);
disp(['a_est = ' mat2str(a_est(end:-1:1)',4)]);
b_est = B * a_est;
b_est = b_est(1:nb+1);
disp(['b = ' mat2str(b,4)]);
disp(['b_est = ' mat2str(b_est',4)]);
disp(' ');
disp(' ');
disp(' ');
h_est = filter(b_est,a_est(end:-1:1),u)';
freqresp(h_est,5,'r--',1);

