%
% stima dei parametri di un modello state-space dalla risposta al campione unitario:
%
N = 32;
% sequenza "campione unitario" di lunghezza N:
u=zeros(1,N);  u(1)=1;
if 0
  % calcolo la risposta del sistema meccanico precedentemente visto:
  param_sistema_meccanico_3gdl;
  build_sistema_meccanico_3gdl;
  load = u';
  simula_sistema_meccanico_3gdl;  % NB: calcola "response".
  h = response;  % in questo caso la risposta del sistema è la risposta al gradino unitario, "h".
  nx=6;
else
  nx=6;
  x0 = zeros(nx,1);
  A = rand(nx,nx);
  b = rand(nx,1);
  c = rand(1,nx);
  h = simula_DLTI_StateSpace_discreto(A,b,c,[0],u,x0);
end
freqresp(h,10); pause,
%
% uso la "h" per identificare un modello state-space di ordine nx=6:
H = hankel(h(2:nx+2), h(nx+2:2*nx+1));
%
% metodo di realizzazione in forma controller-Hessenberg:
[O,R] = qr(H);
O_1 = O(1:nx,1:nx);
O_2 = O(2:nx+1,1:nx);
% Dato che O_1 * A_c = O_2 e le matrici hanno rango pieno, si ha:
A_c = O_1 \ O_2;
% dopodichè:
b_c = R(1:nx,1);
c_c = O(1,1:nx);
% verifico la risposta in frequenza:
u=zeros(1,N);  u(1)=1;
x0 = zeros(nx,1);
h_c = simula_DLTI_StateSpace_discreto(A_c,b_c,c_c,[0],u,x0);
freqresp(h_c,10,'r--',1); pause,
%
% metodo di realizzazione in forma bilanciata:
% Nota: dato che calcola la SVD di H, è dei due il miglior metodo per
% determinare "nx"!
[U,S,V] = svd(H);
sq = sqrt(S(1:nx,1:nx));
O = U(1:nx+1,1:nx)*sq;
R = sq*V(1:nx,1:nx)';
A_b = O(1:nx,1:nx) \ O(2:nx+1,1:nx);
b_b = R(1:nx,1);
c_b = O(1,1:nx);
% verifico la risposta in frequenza:
u=zeros(1,N);  u(1)=1;
x0 = zeros(nx,1);
h_b = simula_DLTI_StateSpace_discreto(A_b,b_b,c_b,[0],u,x0);
freqresp(h_b,10,'g:+',1);

