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;