% Modified Kaelin, Lindgren and Moschytz pole identification
%
% np = number of poles
% M = number of example systems
% S = number of samples
% y = output of the system
% 
% alpha = denominator of transfer function

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

                        % set some initial alpha and e1
e1 = zeros(np+1,1); e1(1,1) = 1;
alpha = e1;
                        
                       % form matrix H
for i = 1:M,
 for j = 1:S - np,
  H(j+(i-1)*(S-np),:) = y(j+np:-1:j,i)';
 end
end

                        % iteraction
for i = 1:10,
 i
                           % form matrix A
 A = zeros(S,S-np);
 for j = 1:S - np,
  A(:,j) = [zeros(j-1,1) ; alpha(np+1:-1:1,1) ; zeros(S-np-j,1)];
 end
 A = [A ];

                         % find F
 F = H'*inv(A'*A)*H;
 invF = inv(F);
 alpha = -invF*e1/invF(1,1);
 roots(-alpha)
end 
 
pmklm = roots(-alpha);       %find the poles
%clear H F invF j k i A S M 




