%   Example matlab file for the paper: Balázs Bank, "Converting IIR
%   filters to parallel form," IEEE Signal Processing Magazine, Tips and
%   tricks, 2018.
%   Creating Fig. 3
%
%   http://www.mit.bme.hu/~bank/parconv/
%
%   C. Balázs Bank, 2018.

load roomresp; % loudspeaker-room response
Fs=44100;
N=6000;
impresp=impresp(225+1:N+225);


NB=200; % the order of the numerator
NA=200; % the order of the denominator

[B,A]=stmcb(impresp,NB,NA);  % designing direct-form filter by the Steighlitz-McBride method
imp=zeros(1,N);
imp(1)=1;
filtresp=filter(B,A,imp); % computing the filter response


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% delayed filter - partial fraction expansion

[Bmd,Amd,FIRd]=tf2delparf(B,A); % converting to the delayed parallel filter by partial fraction expansion
parfrespd=delparfilt(Bmd,Amd,FIRd,imp); % computing the parallel filter impulse reponse

% plotting the results
figure(1);
subplot(2,1,1);
[fr,h]=tfplot(filtresp); % plot the original filter response
semilogx(fr,db(h),'b','LineWidth',2);
hold on;
[fr,hp1]=tfplot(parfrespd); % plot the parallel filter response
semilogx(fr,db(hp1),'r--','LineWidth',2);
w=2*pi*fr/Fs;
Z=exp(-j*w(:)); %creating z^-1
Z2=Z.^2; %creating z^-2
Hsectv=zeros(length(w),size(Amd,2));
for k=1:size(Amd,2), % compute the transfer function of the second-order sections
    Hsectv(:,k) = (Bmd(1,k)+Bmd(2,k)*Z)./(Amd(1,k)+Amd(2,k)*Z+Amd(3,k)*Z2);
end;
semilogx(fr,db(Hsectv),'Linewidth',0.5);  % plot the transfer functions of the second-order sections
HFIR=freqz(FIRd,1,w); % plot the FIR part
semilogx(fr,db(HFIR),'k','Linewidth',1.5); 
axis([20 20000 -60 0]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
puttext('(a)');

start=434; % compute the mean dB error from 20 Hz to Fs/2
e1=mean(abs(db(h(start:end))-db(hp1(start:end))));
fprintf('\nMean abs. dB error for the PFE based conversion, 20Hz to 22kHz: %g dB\n',e1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% delayed filter - least squares procedure

NLS=NB*2; % the length of the impulse response used in the LS fit - must be larger than the number of free parameters NB
[Bmd2,Amd2,FIRd2]=tf2delparf_ls(B,A,NLS); % converting to the delayed parallel form by the LS fit
parfrespd2=delparfilt(Bmd2,Amd2,FIRd2,imp); % computing the parallel filter impulse reponse

% plotting the results
subplot(2,1,2);
[fr,h]=tfplot(filtresp); % plot the original filter response
semilogx(fr,db(h),'b','LineWidth',2);
hold on;
[fr,hp2]=tfplot(parfrespd2); % plot the parallel filter response
semilogx(fr,db(hp2),'r--','LineWidth',2);
w=2*pi*fr/Fs;
Z=exp(-j*w(:)); %creating z^-1
Z2=Z.^2; %creating z^-2
Hsectv=zeros(length(w),size(Amd2,2));
for k=1:size(Amd2,2), % compute the transfer function of the second-order sections
    Hsectv(:,k) = (Bmd2(1,k)+Bmd2(2,k)*Z)./(Amd2(1,k)+Amd2(2,k)*Z+Amd2(3,k)*Z2);
end;
semilogx(fr,db(Hsectv),'Linewidth',0.5); % plot the transfer function of the second-order sections
HFIR=freqz(FIRd2,1,w);
semilogx(fr,db(HFIR),'k','Linewidth',1.5); % plot the transfer function of the FIR part
axis([20 20000 -60 0]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
puttext('(b)');

start=434; % compute the mean dB error from 20 Hz to Fs/2
e2=mean(abs(db(h(start:end))-db(hp2(start:end))));
fprintf('Mean abs. dB error for the LS based conversion, 20Hz to 22kHz: %g dB\n \n',e2);