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

%TF2PARF converts the transfer function form to parallel SOS form.
%   [Bm,Am,FIR]=TF2PARF(B,A) converts the transfer function B(z)/A(z)
%   defined by vectors B and A to 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 PARFILT command.
%
%   Note that the delayed parallel filter obtained by TF2DELPARF is
%   numerically preferred whenever there is a FIR part (NB>=NA).
%
%   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.
%
%   http://www.mit.bme.hu/~bank/delparf
%
%   C. Balazs Bank, 2014.

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