ipf.m

solution for the Statistical modelling for digital signal processing by hayes

function [a,b,err] = ipf(x,p,q,n,a)

%IPF	Pole-zero signal modeling using iterative prefiltering.

%---

%USAGE	[a,b,err] = ipf(x,p,q,n,a) 

%

%	The sequence x is modeled as the unit sample response of

%	a filter having a system function of the form

%		H(z) = B(z)/A(z) 

%	The polynomials B(z) and A(z) are formed from the vectors

%		b=[b(0), b(1), ... b(q)]

%		a=[1   , a(1), ... a(p)]

%	The inputs are

%          x   :  input sequence to be modeled

%          p   :  number of poles in the model

%          q   :  number of zeros in the model

%          n   :  number of iterations to be used to find A(z) and B(z)	

%          a   :  initial estimate for the denominator (default is zero).

%	and the outputs are

%          a   :  vector of coefficients for A(z)

%          b   :  vector of coefficients for B(z)

%          err :  model error

%

%---------------------------------------------------------------

% copyright 1996, by M.H. Hayes.  For use with the book 

% "Statistical Digital Signal Processing and Modeling"

% (John Wiley & Sons, 1996).

%---------------------------------------------------------------



x   = x(:);

N   = length(x);

if p+q>=length(x), error('Model order too large'), end

if nargin < 5 

    a   = prony(x,p,q);

    end

delta = [1; zeros(N-1,1)];



for i=1:n

    f   = filter(1,a,x);

    g   = filter(1,a,delta);

    u   = convm(f,p+1); 

    v   = convm(g,q+1);

    ab  = -[u(1:N,2:p+1) -v(1:N,:) ]\u(1:N,1);

    a   = [1; ab(1:p)];

    b   = ab(p+1:p+q+1);

    err = norm( u(1:N,1) + [u(1:N,2:p+1) -v(1:N,:)]*ab);

    end;