function [Bm,Am,FIR]=tf2delparf(B,A);

%TF2DELPARF converts the transfer function form to delayed parallel SOS form.
%   [Bm,Am,FIR]=TF2DELPARF(B,A) converts the transfer function B(z)/A(z)
%   defined by vectors B and A to delayed parallel second-order sections
%   described by Bm and Am, plus a parallel FIR part, if the order of
%   B(z) is larger than the order of A(z).
%
%   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. These can be used by the filter command separatelly (e.g., by
%   y=filter(Bm(:,3),Am(:,3),x), or by the DELPARFILT command.
%   
%   Note that if there is an FIR part, then the IIR part must be delayed by
%   the length (order+1) of the FIR filter when this form is implemented. 
%   This leads to much better numerical properties compared to the standard
%   parallel filter computed by TF2PARF.
%
%   Also, remember that converting a transfer function to parallel second-order
%   is not possible if there is pole mulitplicity.
%  
%   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.
%
%   Balázs Bank, "Converting IIR filters to parallel form," IEEE Signal Processing
%   Magazine, Tips and tricks, 2018.
%
%   http://www.mit.bme.hu/~bank/parconv
%   http://www.mit.bme.hu/~bank/delparf
%
%   C. Balazs Bank, 2014-2018.

[r,p,k]=residued(B,A); % perform partial fraction expansion to the delayed form
[Bm,Am,FIR]=rpk2parf(r,p,k); % recombine to second-order sections