function [Bm,Am,FIR]=sos2delparf_ls(SOS,G,NLS);

%SOS2DELPARF_LS converts the series SOS form to delayed parallel SOS form.
%   [Bm,Am,FIR]=SOS2DELPARF_LS(B,A,NLS) converts the series second-order form
%   defined in the matrix SOS and gain G 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 parameters are computed by a least-squares fit to match the impulse response
%   of the original filter, instead of partial fraction expansion. This procedure
%   is robust even for large (>100) filter orders.
%   
%   The NLS parameter sets the length of the impulse fit.
%
%   The SOS matrix is defined as used by the TF2SOS and ZP2SOS commands:
%   SOS = [ b01 b11 b21  1 a11 a21 
%           b02 b12 b22  1 a12 a22
%                ...
%           b0L b1L b2L  1 a1L a2L ]
%   And G is the optional series gain coefficient.
%
%   The Bm and Am matrices contain 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. 
%
%   Also, remember that converting a transfer function to parallel second-order
%   is not possible if there is pole mulitplicity.
%  
%   For the details, see:
%
%   Balázs Bank, "Converting IIR filters to parallel form," IEEE Signal Processing
%   Magazine, Tips and tricks, 2018.
%
%   http://www.mit.bme.hu/~bank/parconv
%
%   C. Balazs Bank, 2017-2018.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computing the impulse response of the direct-form filter 

if nargin < 3, % if NLS is not specified
    NLS=2*length(B);
end;

if nargin < 2, % if G is not specified
    G=1;
end;

% computing the impulse response of the series implementation
imp=zeros(1,NLS); 
imp(1)=1;
serresp=imp*G;
NSEC=size(SOS,1); % numer of second-order sections
for k=1:NSEC,
    serresp=filter(SOS(k,1:3),SOS(k,4:6),serresp);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% preparing the target for the LS fit

NFIR=1; % FIR part order is always 1 in the series to parallel case
parfresp=serresp(NFIR+1:end); % the part modeled by the parallel filter
L=length(parfresp);

imp=zeros(L,1);
imp(1)=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%constructing the modeling signal matrix

for k=1:NSEC,
    resp=filter(1,SOS(k,4:6),imp); %impluse response of the denominator
    M(:,2*k-1)=resp;
    M(:,2*k)=[0 ;resp(1:L-1)]; %the response delayed by one sample
end;

y=parfresp(:); %making the target a column vector

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Looking for min(||y-M*par||) as a function of par by the LS fit

par=M\y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% constructing the Bm and Am matrices 

for k=1:NSEC,
    Am(:,k)=SOS(k,4:6).'; % denominators are the same as for the series case
    Bm(:,k)=par(2*k-1:2*k); % numerators from the LS fit
end;

FIR=0; % the FIR part is simply the first NFIR samples of the impulse response
if NFIR>0,
    FIR=serresp(1:NFIR);
end;