% Plots Fig. 8 shown in 'High-Precision Parallel Graphic Equalizer' by J.
% Rämö, V. Välimäki, and B. Bank

% 20.12.2013, Jussi Rämö 

% load impulse responses of the analog, second-order, and fourth-order
% graphic equalizers
load filter_data;

Fs=44100;       % Sample rate
FSTART=20;      % Lowest band
FSTOP=20000;    % Highest band

% EQ center frequencies
eqfr = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 ...
    4000 5000 6300 8000 10000 12500 16000 20000];
BANDS=length(eqfr);
eqfr=[0 eqfr Fs/2]; % Add extra frequency points to prevent target interpolation errors

% Zigzag gains
eqgain = zeros(1,31);
eqgain(1:2:31) = 12;
eqgain(2:2:30) = -12;
eqgain = [eqgain(1) eqgain eqgain(end)];

fir = 1; % Parallel direct path b_o is in use


% creating target response via
% piecewise cubic hermite interpolation
targetfr=logspace(log10(10),log10(Fs/2),BANDS*10);
dBtarget=pchip(eqfr,eqgain,targetfr);
wtarget=2*pi*targetfr/Fs;

polefr = [20 22.4 25 28.2 31.5 35.5 40 44.7 50 56.2 63 70.8 80 89.1 100 112 125 141 160 178 200 224 250 282 315 355 400 447 ...
    500 562 630 708 800 891 1000 1120 1250 1410 1600 1780 2000 2240 2500 2820 3150 3550 4000 4470 5000 5620 6300 7080 8000 ...
    8910 10000 11200 12500 14100 16000 17800 20000];    % Pole frequencies
p=freqpoles([10 polefr]); % Compute poles
POLENUM=length(polefr); % Number of poles

% the design is done on double amount of points compared to pole freqs
designfr=logspace(log10(10),log10(Fs/2),POLENUM*2);
wdesign=2*pi*designfr/Fs;

% making the target minimumphase and resample at design points
minphtarget=minphasen(10.^(dBtarget/20),wtarget,2^14,wdesign);

% Frequency-dependent weghting
Wt = 1./abs(minphtarget).^2;

% designs the second-order sections [Bm,Am] and the FIR coefficient b_0
[Bm,Am,FIR]=parfiltdesfr(wdesign,minphtarget,p,fir,Wt);
parfresp=parfiltfresp(Bm,Am,FIR,wtarget);

% Impulse
len = 330;
imp = zeros(len-1,1);
imp = [1; imp];

% Impulse response of the parallel equalizer
impresp=parfilt(Bm,Am,FIR,imp); 

% Time axis
t = 1:len;
xt = t./Fs*1000;

% Plot figure
figure(1);
plot(xt,filter_data.analog.zigzag(772+36:772+36+len-1),'k','Linewidth',2);
hold on;
plot(xt,filter_data.second.zigzag(1:len)-2,'k','Linewidth',2)
plot(xt,filter_data.fourth.zigzag(1:len)-4,'k','Linewidth',2)
plot(xt,impresp-6,'k','Linewidth',2);

plot([0 6],[0 0],'k')
plot([0 6],[-2 -2],'k')
plot([0 6],[-4 -4],'k')
plot([0 6],[-6 -6],'k')

xlabel('Time (ms)','Fontname','Times','Fontsize',18);
set(gca,'YTick',[])
ylim([-7 1.4])
xlim([0 6])
set(gca,'XTick',0:6,'XTickLabel',0:6,'Fontname','Times','Fontsize',18)
text(3,0.5,'(a)','Fontname','Times','Fontsize',18)
text(3,-1.5,'(b)','Fontname','Times','Fontsize',18)
text(3,-3.5,'(c)','Fontname','Times','Fontsize',18)
text(3,-5.5,'(d)','Fontname','Times','Fontsize',18)
hold off;
