%   Example file for the paper:
%
%   Balazs Bank and Julius O. Smith, III, "A delayed parallel filter structure
%   with an FIR part having improved numerical properties", 136th AES
%   Convention, Preprint No. 9084, Berlin, April 2014.
%
%   http://www.mit.bme.hu/~bank/delparf
%
%   In this example we design a direct-form IIR filter by the
%   Steiglitz-McBride method and convert it to a parallel second-order
%   form. Since NB>=NA, an FIR part of order NB-NA arises, which is in
%   parallel with the IIR sections. However, this leads to numerical
%   problems, since the responses of the sections can be much larger then the 
%   net transfer function (see the paper).
%
%   This can be avoided by an alternative expansion where there is no
%   overlap between the FIR and the IIR parts. In other words, the
%   responses of the second-order sections are delayed. This is shown in
%   the second figure.
%
%   C. Balazs Bank, 2014.


load fe167;
Fs=44100;
N=1000;
impresp=impresp(1:N);

NB=50;
NA=20;

[B,A]=stmcb(impresp,NB,NA);
imp=zeros(1,N);
imp(1)=1;
filtresp=filter(B,A,imp);

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

[Bm,Am,FIR]=tf2parf(B,A);
parfresp=parfilt(Bm,Am,FIR,imp);


figure(1);
[fr,h]=tfplot(impresp);
semilogx(fr,db(h),'b','LineWidth',1.5);
hold on;
[fr,h]=tfplot(parfresp);
semilogx(fr,db(h),'r','LineWidth',1.5);
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),':'); 
HFIR=freqz(FIR,1,w);
semilogx(fr,db(HFIR),'k--'); % FIR part
axis([20 20000 -60 60]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
title('Normal parallel sections');

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

[Bmd,Amd,FIRd]=tf2delparf(B,A);
parfrespd=delparfilt(Bmd,Amd,FIRd,imp);


figure(2);
[fr,h]=tfplot(impresp);
semilogx(fr,db(h),'b','LineWidth',1.5);
hold on;
[fr,h]=tfplot(parfrespd);
semilogx(fr,db(h),'r','LineWidth',1.5);
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),':'); 
HFIR=freqz(FIRd,1,w);
semilogx(fr,db(HFIR),'k--'); % FIR part
axis([20 20000 -60 60]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
title('Delayed parallel sections');



