% PROGRAM PRONY.M
% Prony's method
% Find the poles for the FIXED POLE method
%   inputs:
%    y = example systems
%    np = number of poles that we need to find
%
%   internal variables:
%    M = number of example systems
%    S = number of samples
%    P = are the poles

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

H = zeros(M*(S-np),np);
Y = zeros(M*(S-np),1);

% Set-up and solve the least squares problem

for m=1:M,
 m
 for s=1:S-np,
  H((m-1)*(S-np)+s,:) = y(s+np-1:-1:s,m)';
  Y((m-1)*(S-np)+s,1) = y(s+np,m);
 end
end

A = inv(H'*H)*H'*Y;
pp = roots([1 -A']);
clear H m s S T M Y A
