Zi =filteric(B,A,X,Y)

function Zi =filteric(B,A,X,Y)
%
% The function Y=FILTER(B,A,X,Zi) is defined by the coefficients A and B of the
% filter and the vector Zi of the initial delays, of length
% N = max(length(A), length(B))-1.
%
% Given B,A,X,Y, this function produces the vector Zi such that the output
% of filter(B,A,X,Zi) agrees with Y in the first N terms.
%
% WARNING: The vectors X and Y must be at least of length N.
%

% test that X and Y contain enough terms
lB = length(B);
lA = length(A);
nst  = max(lA,lB) - 1;
if ((length(X) < nst) | (length(Y) < nst)),
  error('Not enough initial conditions');
end;

% check validity of the coefficients. A(1) must be non zero
  if(A(1)==0),
    error('First denominator coefficient must be non-zero.');   
  end; % if(A(1)==0)

% normalize filter coefficients by A(1)
factor=A(1);
A=A/factor;
B=B/factor;

% Truncate X and Y to length nst and make them into column vextors
XIc =X(1:nst);   XIc=XIc(:);
YIc =Y(1:nst);   YIc=YIc(:);

% Build auxiliary variables for the computation of Zi and truncate them
% to length nst

B     =B(:);
cnvbx = conv(B,XIc);
cnvbx = cnvbx([1:nst]');

A(1)=0;
A=A(:);
cnvay = conv(A,YIc);
cnvay = cnvay([1:nst]');
% Compute Zi
Zi    = YIc -cnvbx + cnvay;