% Fig89.m
%
% Plot the comparison between PGE and APGE (third-octave versions). PGE was
% proposed in J. Ramo, V. Valimaki, and B. Bank, “High-precision parallel 
% graphic equalizer,” IEEE/ACM Trans. Audio Speech Lang. Process., vol. 22,
% no. 12, pp. 1894–1904, Dec. 2014 and its optimized version in B. Bank, J.
% A. Belloch, and V. Valimaki, “Efficient design of a parallel graphic 
% equalizer,” J. Audio Eng. Soc., vol. 65, no. 10, pp. 817–825, Oct. 2017.
%
% 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 freqpoles.m, insertsamples.m, interactionMatrix2.m, minphasen.m,
% pareq.m, parfiltdesfr_precomp.m, and parfiltfresp.m.
%
% Created by Juho Liski, Otaniemi, Espoo, Finland, 6 May 2019
% Modified by Juho Liski and Vesa Valimaki, Espoo, Finland, 28 October 2024
% Aalto University, Dept. of Signal Processing and Acoustics

fs  = 44.1e3; % Sample rate

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

if fig_mode == 0
    Gdb = -12*ones(1,31); Gdb(1:2:31) = -Gdb(:,1);  % Zig-zag (dB)
elseif fig_mode == 1
    Gdb = [12 -12 12 12 -12 12 -12 -12 -12 12 -12 12 12 12 -12 12 -12 -12 ...
        12 12 12 -12 -12 12 12 -12 12 -12 -12 12 12];  % Modified Oliver&Jot worst case (dB)
else
    Gdb = [12 12 11 9 6 3 1 0 0 6 6 12 6 6 -12 12 -12 -12 -12 -12 0 0 0 0 ...
        -3 -6 -9 -12 0 0 0]; % Non-extreme case
end

% PGE (by Ramo et al. and optimized by Bank et al.)
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);

eqgain = Gdb;

polefr=insertsamples(eqfr,1,'log'); % insert one frequency point between each eq frequency
pp=freqpoles([10 polefr]); % Compute poles - we add one pole at 10 Hz

NFIR=0; %FIR length
pnum=length(pp); %number of poles
ppnum=2*floor(pnum/2); %the even part of pnum

% DESIGN FREQUENCIES - DOUBLE COMPARED TO POLE FREQUENCIES (four times compared to EQ fr.)
% we use this frequency grid for the LS design - this assures double number
% of specification points compared to free parameters
designfr=insertsamples(polefr,1,'log'); % insert one frequency point between each pole frequency
designfr=[10 16 18 designfr fs/2]; % add some low freq, plus fs/2
wdesign=2*pi*designfr/fs;

% precompute the modeling matrix + denominators
[Am,MR]=parfiltdesfr_precomp(wdesign,pp,NFIR);
L=size(MR,2);

% TARGET FREQUANCIES - TEN TIMES MORE COMPARED TO DESIGN FREQUENCIES
% we have even higher resolution for computing the minimum-phase response
% and for smooth plotting in the figures
TUPS=10; % how many times the target resolution is higher than the design resolution
targetfr=insertsamples(designfr,TUPS-1,'log'); 
wtarget=2*pi*targetfr/fs;

clear hannbasis;
clear hannbasis_design;
for N=1:BANDS
    if N==BANDS
        x = (log10(targetfr)-log10(eqfr(N))) / (log10(eqfr(N))-log10(eqfr(N-1)));
    elseif N==1
        x = (log10(targetfr)-log10(eqfr(N))) / (log10(eqfr(N+1))-log10(eqfr(N)));
    else
        tmp1 = (log10(targetfr)-log10(eqfr(N))) / (log10(eqfr(N+1))-log10(eqfr(N)));
        tmp1 = max(0,tmp1); % positive side
        tmp2 = (log10(targetfr)-log10(eqfr(N))) / (log10(eqfr(N))-log10(eqfr(N-1)));
        tmp2 = min(0,tmp2); % negative side
        x = tmp1 + tmp2; % sum the two
    end
    xm = max(x,-1);
    xm = min(xm,1);    

    if N==1
        xm = max(xm,0);
    end
    if N==BANDS
        xm = min(xm,0);
    end

    hannbasis(N,:)=(cos(pi*xm)+1)/2;
    hannbasis_design(N,:)=hannbasis(N,1:TUPS:end); % downsample according to wdesign frequencies
end

clear phasebasis;
clear phasebasis_design;
for N=1:BANDS
    minph=minphasen(10.^(hannbasis(N,:)/20),wtarget,2^16);
    % phasebasis(N,:)=phase(minph); % NOTE 28 Oct. 2024: phase.m is no longer available
	phasebasis(N,:)=unwrap(angle(minph)); % NOTE 28 Oct. 2024: phase was replaced with unwrap(angle(...))
    phasebasis_design(N,:)=phasebasis(N,1:TUPS:end); % downsample according to wdesign frequencies
end

dBdesign=eqgain*hannbasis_design;
phdesign=eqgain*phasebasis_design;
hanndesign=10.^(dBdesign/20).*exp(1i*phdesign);
ampl=10.^(dBdesign(:)/20);
HR=[ampl.*cos(phdesign(:)); ampl.*sin(phdesign(:))];

% weighting
Wt=1./abs(hanndesign(:)); % making column vector
Wt=[Wt; Wt];
for k=1:L
      MRw(:,k)=MR(:,k).*Wt;
end
HRw=HR.*Wt;

% LS fit
A=MRw'*MRw; 
b=MRw'*HRw;
par=A\b;

%constructing the Bm matrices
for k=1:ppnum/2
    Bm(:,k)=par(2*k-1:2*k);
end

FIR=0;
% constructing the FIR part
if NFIR>0
    FIR=par(pnum+1:pnum+NFIR);
end

% Evaluate frequency response
Nfreq = 2^11; w = logspace(log10(9),log10(22050),Nfreq);
Htot = parfiltfresp(Bm,Am,FIR,pi*w/(fs/2));  % Evaluate frequency response

fontsize = 15;
linewidth = 2; 

% Plot the frequency response
figure(8);
%subplot(1.6,1,1)
semilogx(eqfr,Gdb,'ro','linewidth',2,'MarkerSize',8); hold on;
semilogx(w,db(Htot),'k','LineWidth',linewidth); hold off
set(gca,'fontname','Times','fontsize',fontsize);
if fig_mode == 0
    xlabel({'Frequency (Hz)';'(a)'});
elseif fig_mode == 1
    xlabel({'Frequency (Hz)';'(b)'});
else
    xlabel({'Frequency (Hz)';'(c)'});
end
ylabel('Magnitude (dB)')
set(gca,'YTick',-20:4:20)
set(gca,'XTick',[30 100 300 1000 3000 10000]);
axis([15 22050 -14 14])

% if fig_mode == 0
%     print(gca,'-dpdf','PGE_a')
% elseif fig_mode == 1
%     print(gca,'-dpdf','PGE_b')
% else
%     print(gca,'-dpdf','PGE_c')
% end


% APGE
fc = [1000./(2.^((1/3)*(17:-1:0))) 1000.*(2.^((1/3)*(1:13)))];
[B1,A1,F,~,~,~] = ACGE2APGE(Gdb);

%%% Initialization
W = 2*pi*w/fs; W  = W(:);
Z  = exp(-1i*W); % Creating z^-1
Z2 = Z.^2; % Creating z^-2
NFIR = length(F);

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

%%% Initialization
Hpartot = zeros(length(W),1);
Hpar = zeros(length(W),NSEC);

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

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

% Plot responses for the proposed optimized design
figure(9);
%subplot(1.6,1,1)
semilogx(fc,Gdb,'ro','linewidth',2,'markersize',8); hold on
semilogx(w,db(Hpartot),'k','linewidth',linewidth); hold off 
set(gca,'fontname','Times','fontsize',fontsize);
if fig_mode == 0
    xlabel({'Frequency (Hz)';'(a)'});
elseif fig_mode == 1
    xlabel({'Frequency (Hz)';'(b)'});
else
    xlabel({'Frequency (Hz)';'(c)'});
end
ylabel('Magnitude (dB)')
set(gca,'XTick',[30 100 300 1000 3000 10000]);
set(gca,'YTick',-20:4:20)
axis([15 22050 -14 14])

% if fig_mode == 0
%     print(gca,'-dpdf','APGE_a')
% elseif fig_mode == 1
%     print(gca,'-dpdf','APGE_b')
% else
%     print(gca,'-dpdf','APGE_c')
% end

