% Model reduction using Hankel norm
% Create matrix H from Beliczynski, Kale and Cain article  (SM11.M)
%    np = number of poles
[r,c]=size(y);
for j=1:c,
 j
 H = zeros(r-1);
 for i=1:r-1,			% skip c0
  H(i,:) = [y(i+1:r,j)' zeros(1,i-1)];
 end
 disp('Find eigenvalues...')
 [v,e]=eig(H);
% plot the eigenvalues
% for k=1:r-1,
%  e(k) = e(k,k);
% end
% plot(abs(e(r-1:-1:100)));
% disp('.')
% pause
% Find new system
 ak=v(2:r-1,r-np:r-1)'*v(1:r-2,r-np:r-1);
 bk=v(1,r-np:r-1)';
 ck=y(2:r,j)'*v(1:r-1,r-np:r-1);
 [nume(j,:),dene(j,:)]=ss2tf(ak,bk,ck,y(1,j),1);
 clear H v e k ak bk ck
 yi(:,j)=dimpulse(nume(j,:),dene(j,:),r);
 plot(1:r,y(:,j),'-',1:r,yi(:,j),':')
 disp('.')
 pause
end
  
