%   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. 4.
%
%   http://www.mit.bme.hu/~bank/parconv/
%
%   C. Balázs Bank, 2018.


load sosfilt; % load the filter in the series second-order form

imp=zeros(1,10000); % create a unit pulse
imp(1)=1;

serresp=imp*G; % compute the impulse response of the series filter
for k=1:size(SOS,1),
    serresp=filter(SOS(k,1:3),SOS(k,4:6),serresp);
end;


[Bmd,Amd,FIRd]=sos2delparf_ls(SOS,G,10000); % perform the series -> parallel conversion
parfrespd=delparfilt(Bmd,Amd,FIRd,imp); % compute the parallel filter response

% display the results
[fr,h]=tfplot(serresp);
semilogx(fr,db(h),'b','LineWidth',1.5);
hold on;
[fr,hp]=tfplot(parfrespd);
semilogx(fr,db(hp)-20,'r','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('Series implementation',0.69,0.86,SIZE); 
puttext('Parallel implementation',0.68,0.525,SIZE); 

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