function y=pefilter(b,a,x,ybeg);

% function y=pefilter(b,a,x,ybeg);
% y is the output of the filter H(z)=B(z)/A(z)
% The initial conditions are chosen so that
% y(1:L)=ybeg, L=length(ybeg)
%
% Programmed by: Dimitris Manolakis, 10/4/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.
%-----------------------------------------------------------



if nargin<4, ybeg=[]; end
p=length(a)-1;
q=length(b)-1;
M=max([p q length(ybeg)]);
ybeg=[ybeg' zeros(1,M)]';

y1=filter(a,1,ybeg(1:M));
y2=filter(b,1,x(1:M));
w=y1-y2;
b1=[b' zeros(1,M-q)]';
a1=[a' zeros(1,M-p)]';
y=filter(b1,a1,x,w);