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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% designing the Butterworth filter in the series form

Fs=44100; % sampling frequency
N=8; % filter order
fc=100; % cutoff frequency

[z,p,gain]=butter(N,fc/(Fs/2),'high');
[SOS,G]=zp2sos(z,p,gain); % series version

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

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% zpk conversion

z2=z(1:end-1); % removing one zero - now all are -1, so this is easy
% otherwise, we have to be careful to remove a real zero or a complex conjugate pole pair
% 

[r,p]=zp2rp(z2,p,gain); % converting the pole-zero form to residue form 
[Bm,Am]=rpk2parf(r,p,0); % recombining the residues to second-order sections


parfresp=parfilt(Bm,Am,0,imp); % the response of the parallel filter
parfresp=filter([1 -z(end)],1,parfresp); % putting back effect of the zero 

% plotting the response of the series implementation
[fr,h]=tfplot(serresp);
semilogx(fr,db(h),'b','LineWidth',2);
hold on;
% plotting the response of the parallel implementation
[fr,h]=tfplot(parfresp);
semilogx(fr,db(h),'r--','LineWidth',2);

% plotting the responses of the second-order sections + the one zero
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);
    % the effect of the zero that was factored out
    Hsectv(:,k) = (1 - Z*z(end)).* Hsectv(:,k) ; 
end;
semilogx(fr,db(Hsectv),'Linewidth',0.75); 
hold off;

axis([10 20000 -150 50]);
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');