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

%TF2DELPARF_LS converts the transfer function form to delayed parallel SOS form.
%   [Bm,Am,FIR]=TF2DELPARF_LS(B,A,NLS) 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 difference from TF2DELPARF is that now 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
%   suggested for large (>100) filter orders due to its better numerical
%   preformance. The NLS parameter sets the length of the impulse response fit.
%
%   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:
%
%   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;

imp=zeros(1,NLS); 
imp(1)=1;
impresp=filter(B,A,imp); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computing the poles of the filter

p=roots(A); % computing the poles of the filter
% if there was an error in root finding giving unstable poles:
% we make the filter stable by flipping the poles into the unit cirle
for k=1:length(p), 
   if abs(p(k))>1			
      p(k)=1/conj(p(k));
  end;
end;

p=cplxpair(p); %order it to complex pole pairs + real ones afterwards

%in order to have second-order sections only (i.e., no first order)
pnum=length(p); %number of poles
ppnum=2*floor(pnum/2); %the even part of pnum
ODD=0;
if pnum>ppnum, %if the number of poles is odd
    ODD=1;
end;

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

NFIR=max(length(B)-length(A)+1,0); % FIR part order
parfresp=impresp(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:2:ppnum-1, %second-order sections
    resp=filter(1,poly(p(k:k+1)),imp); %impluse response of the two-pole filter
    M(:,k)=resp;
    M(:,k+1)=[0 ;resp(1:L-1)]; %the response delayed by one sample
end;

if ODD, %if the number of poles is odd, we have a first-order section
    resp=filter(1,poly(p(pnum)),imp);
    M(:,pnum)=resp;
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 from par

for k=1:ppnum/2, % second-order sections
    Am(:,k)=poly(p(2*k-1:2*k)).';
    Bm(:,k)=par(2*k-1:2*k);
end;
if ODD, %we extend the first-order section to a second-order one by adding zero coefficients
    Am(:,k+1)=[poly(p(pnum)).'; 0];
    Bm(:,k+1)=[par(pnum); 0];
end;

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