% Figs67.m
%
% Plot the comparison between EQ proposed in Z. Chen, G. S. Geng, F. L. 
% Yin, and J. Hao, “A pre-distortion based design method for digital audio 
% graphic equalizer,” Digital Signal Process., vol. 25, pp. 296–302, Feb. 
% 2014. and APGE (octave versions).
%
% This was done for the article J. Liski, B. Bank, J. O. Smith, and V. 
% Valimaki, "Converting series biquad filters into delayed parallel form: 
% Application to graphic equalizers," IEEE Transactions on Signal 
% Processing, 2019.
%
% Uses interactionMatrix2.m and pareq.m
%
% Created by Juho Liski, Otaniemi, Espoo, Finland, 6 May 2019
%
% Aalto University, Dept. of Signal Processing and Acoustics

% Select (a) or (b) plot:
fig_mode = 1; % (a)
%fig_mode = 0; % (b)

% EQ by Chen et al.
wdk = 2*pi*[31.25 62.50 125 250 500 1000 2000 4000 8000 16000]/48000;
Qk = [1.4142 1.4142 1.4141 1.4139 1.4130 1.4092 1.3941 1.3350 1.1156 0.4738];
if fig_mode
    Gk = [1.2502 1.2208 1.2229 1.2229 1.2220 1.2182 1.2032 1.1476 0.9364 1.5759]; % (a)
else
    Gk = [1.2632 0.1947 1.0613 0.9837 1.0882 0.2209 1.0781 1.0543 -0.0336 1.3828]; % (b)
end

ChenNum = zeros(10,3);
ChenDen = zeros(10,3);
for k = 1:10
    ChenNum(k,:) = Gk(k)*[sin(wdk(k)) 0 -sin(wdk(k))];
    ChenDen(k,:) = [2*Qk(k)+sin(wdk(k)) -4*Qk(k)*cos(wdk(k)) 2*Qk(k)-sin(wdk(k))];
end

Nfreq = 2^11;
w2 = logspace(log10(9),log10(24000),Nfreq);
W2 = 2*pi*w2/48000; W2  = W2(:);
Zz  = exp(-1i*W2); % Creating z^-1
Zz2 = Zz.^2; % Creating z^-2

H_Chentot = zeros(length(W2),1);
H_Chenpar = zeros(length(W2),10);

for k = 1:10
    H_Chenpar(:,k) = (ChenNum(k,1)+ChenNum(k,2)*Zz+ChenNum(k,3)*Zz2)./(ChenDen(k,1)+ChenDen(k,2)*Zz+ChenDen(k,3)*Zz2);
    H_Chentot = H_Chentot + H_Chenpar(:,k);
end

%%% Design APGE
if fig_mode
    Gdb = [5 5 5 5 5 5 5 5 5 5]; % (a)
else
    Gdb = [3 -3 3 3 3 -3 3 3 -3 3]; % (b)
end

fs  = 44.1e3;  % Sample rate
fc1 = 16000./(2.^(9:-1:0)); % Exact log center frequencies for filters

[B1,A1,F,~,~,~] = ACGE2APGE(Gdb);

%%% Number of sections and length of FIR part
NSEC = length(B1);
NFIR = length(F);

%%% Initialization
w = logspace(log10(9),log10(22050),Nfreq);
W = 2*pi*w/fs; W  = W(:);
Z  = exp(-1i*W); % Creating z^-1
Z2 = Z.^2; % Creating z^-2
H1partot = zeros(length(W),1);
H1par = zeros(length(W),NSEC);

%%% Evaluate sections separately and when combined
for k = 1:NSEC
    H1par(:,k) = (B1(1,k)*Z+B1(2,k)*Z2)./(A1(1,k)+A1(2,k)*Z+A1(3,k)*Z2);
    H1partot = H1partot + H1par(:,k);
end

%%% Adding the FIR part
for k = 1:NFIR
    H1partot = H1partot + F(k); 
end

fontsize = 15;
linewidth = 2; 

figure(6);
%subplot(1.6,1,1)
semilogx(fc1,Gdb,'ro','linewidth',2,'markersize',8); hold on
semilogx(w2,db(H_Chentot),'k','linewidth',linewidth); hold off 
set(gca,'fontname','Times','fontsize',fontsize);
if fig_mode
    xlabel({'Frequency (Hz)';'(a)'});
else
    xlabel({'Frequency (Hz)';'(b)'});
end
ylabel('Magnitude (dB)')
set(gca,'XTick',[30 100 300 1000 3000 10000]);
set(gca,'YTick',-20:2:20)
axis([15 22050 -5 7])

% if fig_mode
%     print(gca,'-dpdf','Chen_a')
% else
%     print(gca,'-dpdf','Chen_b')
% end

figure(7);
%subplot(1.6,1,1)
semilogx(fc1,Gdb,'ro','linewidth',2,'markersize',8); hold on
semilogx(w,db(H1partot),'k','linewidth',linewidth); hold off 
set(gca,'fontname','Times','fontsize',fontsize);
if fig_mode
    xlabel({'Frequency (Hz)';'(a)'});
else
    xlabel({'Frequency (Hz)';'(b)'});
end
ylabel('Magnitude (dB)')
set(gca,'XTick',[30 100 300 1000 3000 10000]);
set(gca,'YTick',-20:2:20)
axis([15 22050 -5 7])

% if fig_mode
%     print(gca,'-dpdf','APGE_a')
% else
%     print(gca,'-dpdf','APGE_b')
% end
