close all;
clear all;
% Definir o sistema

%Num = [0 1]
%Den = [1 -0.9625]
%Num = [0 0 1]
%Den = [1 -1.2115 .2406]
Num = [0 1 -.2 .01]
Den = [1 -1.7115 .8464 -.1203]
%Encontra a resposta impulsiva com n pontos

n = 300;
y1 = dimpulse(Num, Den, n);

%Imprime a resposta
figure, stem(y1);

%Cria o sistema
MatrizA = diag(ones(n-2,1),-1);
MatrizB = zeros(n-1,1); MatrizB(1) = 1;
MatrizC = y1(2:n)';
MatrizD = y1(1);

y2 = dimpulse(MatrizA,MatrizB,MatrizC,MatrizD,1, n);    

%figure, stem(y2);

H = zeros(n-1);
%H(1,1:n-1) = MatrizC(1:n-1);
%H(2,1:n-2) = MatrizC(2:n-1);
%H(3,1:n-3) = MatrizC(3:n-1);
%H(4,1:n-4) = MatrizC(4:n-1);
for indice=1:n-1
   H(indice,1:n-indice) = MatrizC(indice:n-1);
end

[u,s,v] = svd(H);

k=3
Ak = v(2:n-1,1:k)'*v(1:n-2,1:k);
Bk = v(1,1:k)';
Ck = MatrizC*v(1:n-1,1:k);

[Num2, Den2] = ss2tf(Ak, Bk, Ck, 0,1)	 

[Z,P,K] = ss2zp(Ak,Bk,Ck,0,1)

sing = diag(s);
dsing = diff(sing);

figure, stem(sing);
figure, stem(dsing);

ordemdiff = find(dsing==min(dsing))
ordemP = find(sing>((max(sing)-min(sing))*0.001))