function y=delparfilt(Bm,Am,FIR,x);
%DELPARFILT filtering by parallel second-order sections.
% Y=DELPARFILT(Bm,Am,FIRcoeff,X) filters the signal X by the parallel
% second-order sections given in the matrices Bm and Am. Additionally, X
% is filtered by an FIR filter given by FIRcoeff, and this is also added
% to the output. If there is no parallel FIR part (not even a simple
% coefficient), then use FIRcoeff=0 as a parameter.
%
% The Bm and Am matrices are containing the [b0 b1]' and [1 a0 a1]'
% coefficients for the different sections in their columns. For example,
% Bm(:,3) gives the [b0 b1]' parameters of the third second-order
% section.
%
% An important difference from PARFILT is that if there is an FIR part,
% then the IIR part is delayed so that there is no overlap between the two.
% This leads to improved numerical properties.
%
% The delayed parallel filter can be designed by the DELPARFILTDES command,
% or converted from a normal parallel filter by PARF2DELPARF.
%
% For the details, see:
%
% Balazs Bank and Julius O. Smith, III, "A delayed parallel filter structure
% with an FIR part having improved numerical properties", 136th AES
% Convention, Preprint No. 9084, Berlin, April 2014.
%
% http://www.mit.bme.hu/~bank/delparf
%
% C. Balazs Bank, 2014.
NFIR=length(FIR);
if NFIR==1, % if the FIR coeff is zero then there is no FIR part
if FIR==0,
NFIR=0;
end;
end;
y=zeros(size(x));
s=size(Am);
for k=1:s(2),
y(NFIR+1:end)=y(NFIR+1:end)+filter(Bm(:,k),Am(:,k),x(1:end-NFIR));
end;
if NFIR>0,
y=y+filter(FIR,1,x);
end;