%   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 fixed-pole parallel filter with logarithmic
%   pole positioning having an additional FIR part. Unfortunatelly this
%   can lead to numerical problems, since the responses of the sections can be much larger then the 
%   net transfer function (see the paper). This is shown in Fig. 1.
%
%   This can be avoided if there is no overlap between the FIR and the IIR parts.
%   In other words, the responses of the second-order sections are delayed. 
%   The parameters of the delayed parallel filter can be obtained by conversion
%   from the original parallel filter parameters with PARF2DELPARF. The responses
%   are shown in Figure 2.
%
%   Alternatively, we can also design the filter directly in the delayed
%   form, which is shown in Figure 3. In theory, this is equivallent to
%   conversion, but in practice this is the preferred way because the
%   problem is numerically better conditioned and the estimation is faster
%   (the dimension of the LS problem to be solved is smaller).
%
%   C. Balazs Bank, 2014.


Fs=44100;
load fe167; % speaker
N=length(impresp);

P=30;
NFIR=20;
fp=logspace(log10(20),log10(20000),P);
p=freqpoles(fp,Fs);

imp=zeros(1,N);
imp(1)=1;
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% normal parallel filter

[Bm,Am,FIR]=parfiltdes(impresp,p,NFIR); % normal parallel filter
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--');
axis([20 20000 -60 40]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
title('Normal parallel filter');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% delayed parallel filter by conversion

[Bmc,Amc,FIRc]=parf2delparf(Bm,Am,FIR); % converting to a delayed one
parfrespc=delparfilt(Bmc,Amc,FIRc,imp);

figure(2);
[fr,h]=tfplot(impresp);
semilogx(fr,db(h),'b','LineWidth',1.5);
hold on;
[fr,h]=tfplot(parfrespc);
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(Amc,2));
for k=1:size(Amc,2), %second-order sections
    Hsectv(:,k) = (Bmc(1,k)+Bmc(2,k)*Z)./(Amc(1,k)+Amc(2,k)*Z+Amc(3,k)*Z2);
end;
semilogx(fr,db(Hsectv),':'); 
HFIR=freqz(FIRc,1,w);
semilogx(fr,db(HFIR),'k--');
axis([20 20000 -60 40]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
title('Delayed parallel filter by conversion');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% delayed parallel filter by direct design

[Bmd,Amd,FIRd]=delparfiltdes(impresp,p,NFIR); % normal parallel filter
parfrespd=delparfilt(Bmd,Amd,FIRd,imp);

figure(3);
[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--');
axis([20 20000 -60 40]);
hold off;
SIZE=14;
set(gca,'FontName','Times','Fontsize',SIZE);
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
title('Delayed parallel filter by direct design');