N = 16;
%
nx=6;
u = rand(1,N);
x0 = ones(nx,1);
A = rand(nx,nx);
b = rand(nx,1);
%c = rand(1,nx);
%
% vediamo per quali combinazioni di misura degli stati il sistema è osservabile :
%
% caso di 1 uscita:
for i=1:nx
  c = zeros(1,nx); c(1,i)=1;
  OSS = c;
  for k=1:nx-1
    OSS = [c; OSS*A];
  end;
  disp([num2str(c) ' :  indice di osservabilita'' = ' num2str(rank(OSS))]);
end;
disp(' ');

%
%
%
nx=6;
u = 1*rand(1,N);
x0 = ones(nx,1);
[Q,X]=qr(randn(nx));
D=diag([0.1*(1:nx)]);
A=Q*D*Q';
b = rand(nx,1);
c = rand(1,nx);
%y = simula_DLTI_StateSpace_discreto(A,b,c,[0],u,x0);
%
% calcola il guadagno dell'osservatore 
%
Ko = place(A', c', 0.1*(1:nx))';
Ao = A - Ko*c;
eig(Ao)
%
x0_est = zeros(nx,1);
x = x0;
x_est = x0_est;
figure(1); clf, hold on,
  plot(0,x,'b+');
  plot(0,x_est,'r.');
for i=1:N
  y(i) = c * x;
  x = A * x  +  b * u(i);
  x_est = Ao * x_est + b * u(i) + Ko * y(i);
  figure(1); 
    plot(i,x,'b+');
    plot(i,x_est,'r.');
  %  keyboard
end


