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

% 20.12.2013, Jussi Rämö & Balazs Bank 

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

% Set the equalizer sliders in dB
eqgain=zeros(1,BANDS);
eqgain(1)=12;
eqgain(2)=12;
eqgain(3)=11;
eqgain(4)=9;
eqgain(5)=6;
eqgain(6)=3;
eqgain(7)=1;
eqgain(10:14)=6;
eqgain(12)=12;
eqgain(15:20)=-12;
eqgain(16)=12;
eqgain(25)=-3;
eqgain(26)=-6;
eqgain(27)=-9;
eqgain(28)=-12;

eqgain=[eqgain(1) eqgain eqgain(end)];

% 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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SAME NUMBER OF POLES AS BANDS with no weighting


polefr=eqfr(2:end-1);   % Pole frequencies
p=freqpoles(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);

% designs the second-order sections [Bm,Am] and the FIR coefficient b_0
[Bm,Am,FIR]=parfiltdesfr(wdesign,minphtarget,p,1); % No frequency-dependent weights
parfresp=parfiltfresp(Bm,Am,FIR,wtarget); % Frequency response of second-order parallel filters

% Plot figure
figure(1)
subplot(2,2,1);
semilogx(targetfr,dBtarget,'Color',[0.5 0.5 0.5],'LineWidth',4);
hold on;
semilogx(eqfr,eqgain,'ok','LineWidth',1.5);
semilogx(targetfr,db(parfresp),'k','LineWidth',2);
hold off;
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
axis([20 20000 -20 15]);
set(gca,'XTick',[100 1000 10000]);
set(gca,'XTickLabel',{'100','1k ','10k'},'Fontname','Times','Fontsize',18)
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SAME NUMBER OF POLES AS BANDS with weighting

[Bm,Am,FIR]=parfiltdesfr(wdesign,minphtarget,p,1,1./abs(minphtarget).^2); % With weighting
parfresp=parfiltfresp(Bm,Am,FIR,wtarget);

subplot(2,2,2);
semilogx(targetfr,dBtarget,'Color',[0.5 0.5 0.5],'LineWidth',4);
hold on;
semilogx(eqfr,eqgain,'ok','LineWidth',1.5);
semilogx(targetfr,db(parfresp),'k','LineWidth',2);
hold off;
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
axis([20 20000 -20 15]);
set(gca,'XTick',[100 1000 10000]);
set(gca,'XTickLabel',{'100','1k','10k'},'Fontname','Times','Fontsize',18)
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DOUBLE NUMBER OF POLES AS BANDS with no weighting

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
polefr=[10 polefr]; % Add pole frequency to 10 Hz
p = freqpoles(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);

[Bm,Am,FIR]=parfiltdesfr(wdesign,minphtarget,p,1); % no weights
parfresp=parfiltfresp(Bm,Am,FIR,wtarget);

subplot(2,2,3);
semilogx(targetfr,dBtarget,'Color',[0.5 0.5 0.5],'LineWidth',4);
hold on;
semilogx(eqfr,eqgain,'ok','LineWidth',1.5);
semilogx(targetfr,db(parfresp),'k','LineWidth',2);
hold off;
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
axis([20 20000 -20 15]);
set(gca,'XTick',[100 1000 10000]);
set(gca,'XTickLabel',{'100','1k','10k'},'Fontname','Times','Fontsize',18)
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DOUBLE NUMBER OF POLES AS BANDS with weighting

[Bm,Am,FIR]=parfiltdesfr(wdesign,minphtarget,p,1,1./abs(minphtarget).^2); % With weights
parfresp=parfiltfresp(Bm,Am,FIR,wtarget);

subplot(2,2,4);
semilogx(targetfr,dBtarget,'Color',[0.5 0.5 0.5],'LineWidth',4);
hold on;
semilogx(eqfr,eqgain,'ok','LineWidth',1.5);
semilogx(targetfr,db(parfresp),'k','LineWidth',2);
hold off;
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');
axis([20 20000 -20 15]);
set(gca,'XTick',[100 1000 10000]);
set(gca,'XTickLabel',{'100','1k','10k'},'Fontname','Times','Fontsize',18)
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')