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




% THIS CREATES FIG 1
NB=20; % order of the numerator
NA=20; % order of the denominator

% THIS CREATES FIG 2
% NB=25; % order of the numerator
% NA=20; % order of the denominator


% anechoic response of the Fostex FE-167E full-range loudspeaker
load fe167; 
Fs=44100; % sampling rate
N=1000;
impresp=10*impresp(1:N); % we use the first 1000 samples



[B,A]=stmcb(impresp,NB,NA); % filter design with the Steiglitz-McBride algorithm
imp=zeros(1,N);
imp(1)=1;
filtresp=filter(B,A,imp); % computing the impulse response of the direct-form filter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% non-delayed parallel filter

[Bm,Am,FIR]=tf2parf(B,A); % conversion to the traditional parallel form
parfresp=parfilt(Bm,Am,FIR,imp); % computing the impulse response of the parallel filter


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,h]=tfplot(parfresp); % plot the non-delayed parallel filter response
semilogx(fr,db(h),'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(Am,2));
for k=1:size(Am,2), %second-order sections
    Hsectv(:,k) = (Bm(1,k)+Bm(2,k)*Z)./(Am(1,k)+Am(2,k)*Z+Am(3,k)*Z2);
end;
semilogx(fr,db(Hsectv),'Linewidth',0.75);  % plot the transfer function of the second-order sections
HFIR=freqz(FIR,1,w); % plot the FIR part
semilogx(fr,db(HFIR),'k--','Linewidth',1.5); 
ax=axis; % set the axis of the plot
upperlim=max(max(db(Hsectv)));
upperlim=max(upperlim,max(db(HFIR)));
upperlim=10*ceil(upperlim/10);
axis([20 20000 -30 upperlim]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
puttext('(a)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% delayed parallel filter

[Bmd,Amd,FIRd]=tf2delparf(B,A); % conversion to the delayed parallel form
parfrespd=delparfilt(Bmd,Amd,FIRd,imp); % computing the delayed parallel filter response


figure(1);
subplot(2,1,2);
[fr,h]=tfplot(filtresp); % plot the original filter response
semilogx(fr,db(h),'b','LineWidth',2);
hold on;
[fr,h]=tfplot(parfrespd);  % plot the delayed parallel filter response
semilogx(fr,db(h),'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(Am,2));
for k=1:size(Amd,2), %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.75); % plot the transfer function 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 -30 upperlim]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
puttext('(b)');