function [b,a] = prony(h, nb ,na)
%PRONY	Prony's method for time-domain IIR filter design.
%	[B,A] = PRONY(H, NB, NA) finds a filter with numerator order
%	NB, denominator order NA, and having the impulse response in
%	vector H.   The IIR filter coefficients are returned in
%	length NB+1 and NA+1 row vectors B and A, ordered in
%	descending powers of Z.
%
%	See also: BUTTER, CHEBY1, CHEBY2, ELLIP, INVFREZ.

%	L. Shure 5-17-88
%	Revised LS 17-Dec-90
%	(c) Copyright 1988, by The MathWorks, Inc.

%	References:
%	  [1] T.W. Parks and C.S. Burrus, Digital Filter Design,
%	      John Wiley and Sons, 1987, p226.

K = length(h) - 1;
M = nb; N = na;
c = h(1);
H = tril(toeplitz(h/c));	
% K+1 by N+1
if (K > N)
	H(:,(N+2):(K+1)) = [];
end
% Partition H matrix
H1 = H(1:(M+1),:);	% M+1 by N+1
h1 = H((M+2):(K+1),1);	% K-M by 1
H2 = H((M+2):(K+1),2:(N+1));	% K-M by N
a = [1; -H2\h1].';
b = c*a*H1.';

