% Model reduction using Hankel norm
% Create matrix H from Beliczynski, Kale and Cain article  (SM12.M)
[r,c]=size(y);
disp('Find poles for each example system....')
disp('Find eigenvalues...')
for j=1:c,
 disp(j)
 H = zeros(r-1);
 for i=1:r-1,			% skip c0
  H(i,:) = [y(i+1:r,j)' zeros(1,i-1)];
 end
 [v,e]=eig(H);
% plot the eigenvalues
% for k=1:r-1,
%  e(k) = e(k,k);
% end
% plot(abs(e(100:-1:40)));
% pause
% find the state space description
 ak=v(2:99,90:99)'*v(1:98,90:99);
 bk=v(1,90:99)';
 ck=y(2:100,j)'*v(1:99,90:99);
 [nume(j,:),dene(j,:)]=ss2tf(ak,bk,ck,0,1);
%plot(1:121,y(:,j),1:121,dimpulse(nume(j,:),dene(j,:),121))
%pause
end
clear H ak bk ck v e
%  Create a new system [a b c d] in the following way:
%        .  
%       |x1|   |A1  0 .. 0||x1|   |B1|
%       |x2| = | 0 A2    0||x2| + |B2|u
%       |..|   | ........ ||..|   |..|
%       |xn|   | 0  0 ..An||xn|   |Bn|
%
%       |y1|   |C1  0 .. 0||x1|   |D1|
%       |y2| = | 0 C2 .. 0||x2| + |D2|u
%       |..|   | ........ ||..|   |..|
%       |yn|   | 0  0 ..Cn||xn|   |Dn|
%
%
%   INPUTS: matrices NUM,DEN
%

[dr,dc]=size(dene);

[a1,b1,c1,d1] = tf2ss(nume(1,:),dene(1,:));
[a1,b1]=d2c(a1,b1,1);

for j=2:dr
 [a2,b2,c2,d2] = tf2ss(nume(j,:),dene(j,:));
 [a2,b2]=d2c(a2,b2,1);
 [a1,b1,c1,d1] = mappend(a1,b1,c1,d1,a2,b2,c2,d2); 
end
[ax,bx,cx,dx,ay,by,cy,dy,aug]=ohkapp(a1,b1,c1,d1,1,10);
[ad,bd]=c2d(ax,bx,1);
ph=eig(ad);
clear a1 b1 c1 d1 nume dene
