function [a,b,V,FPE]=pzls(y,x,p,q)
% function [a,b,V,FPE]=pzls(y,x,p,q)
% Pole-zero modeling using equation-error
% (=> linear) least-squares
% Model parameters [1 a1...ap 1 b1...bq V]
% V=input signal variance.
% FPE=Akaike's final prediction error criterion
%
% Programmed by: Dimitris Manolakis, 9/28/93
%
%-----------------------------------------------------------
% Copyright 2000, by Dimitris G. Manolakis, Vinay K. Ingle,
% and Stephen M. Kogon.  For use with the book
% "Statistical and Adaptive Signal Processing"
% McGraw-Hill Higher Education.
%-----------------------------------------------------------



Lx=length(x);
pqmax=max(p,q);
pq=p+q;
R=zeros(pq,pq);
d=zeros(pq,1);
for n=pqmax+1:Lx
   z=[-y(n-p:n-1)' x(n-q:n-1)']';
   R=R+z*z';
   d=d+y(n)*z;
end

c=R\d;
a=[1 flipud(c(1:p))']';
b=[1 flipud(c(p+1:p+q))']';

% Computes the equation error like pe.m

b(1)=0;
e=filter(a,1,y)-filter(b,1,x);
for n=1:pqmax e(n)=0; end;    
V=(e'*e)/(Lx-pqmax);
FPE=V*(1+pq/Lx)/(1-pq/Lx);

b(1)=1;