% Kaelin, Lindgren and Moschytz pole identification (SM6.M)
%
% np = number of poles
% M = number of example systems
% S = number of samples
% y = output of the system
% u = input
% alpha = denominator of transfer function

S = size(y); M = S(1,2); S = S(1,1); 

                        % set some initial alpha and e1
e1 = zeros(2*np+1,1); e1(1,1) = 1;
alpha = zeros(np+1,1); alpha(1,1) = 1; 
beta = zeros(np,1); beta(1,1) = 1;
                        % iteraction
for i = 1:10,
 i
 F = zeros(2*np+1);
 for k = 1:M,                        
                           % form matrix A
  for j = 1:S - np,
   A(:,j)=[zeros(j-1,1);alpha(np+1:-1:1,1);beta(np:-1:1,1);zeros(S-np-j,1)];
  end
                        % form matrix H
  for j = 1:S - np,
   H(j,:) = [y(j+np:-1:j,k)' u(j+np:-1:j+1,1)'];
  end
                         % find F
  F = H'*inv(A'*A)*H + F;
 end
 invF = inv(F);
 v = -invF*e1/invF(1,1);
 alpha = v(1:np+1,1); beta = v(np+2:2*np+1,1);
 roots(-alpha)
end 
 
pklm = roots(-alpha);       %find the poles
clear H F invF j k i A S M 




