% Fig3.m
%
% Plot the amplitude response of the new octave EQ design with two
% different command gain settings.
% This was done for the article V. Valimaki and J. Liski, "Accurate Cascade
% Graphic Equalizer," IEEE Signal Processing Letters, 2017
%
% Uses pareq.m & interactionMatrix.m
%
% Created by Juho Liski and Vesa Valimaki, Otaniemi, Espoo, Finland, 26 September 2016
% Modified by Juho Liski, Otaniemi, Espoo, Finland, 10 February 2017
%   - Minor bug fix
%
% Aalto University, Dept. of Signal Processing and Acoustics

fs  = 44.1e3;  % Sample rate
fc1 = 16000./(2.^(9:-1:0)); % Exact log center frequencies for filters
fc2 = zeros(1,19); fc2(1:2:19) = fc1; % Extra points between filters (Hz)
for k = 2:2:19
    fc2(k) = sqrt(fc2(k-1)*fc2(k+1));  % Extra points are at geometric mean frequencies
end
wg = 2*pi*fc1/fs;  % Command gain frequencies in radians
wc1 = 2*pi*fc1/fs;  % Center frequencies in radians for non-iterative design at command frequencies
wc2 = 2*pi*fc2/fs;  % Center frequencies in radians for iterative design with extra points
gw = 0.3; % Gain factor at bandwidth
bw = 2*pi*[46.88 93.75 187.5 375.0 750.0 1500 3000 5580 9360 12160]/fs; % Bandwidths in radians

load('interaction_matrices.mat') % Load two interaction matrices that have been calculated
                                 % in advance and stored:
                                 % - leak1 is 10x10 matrix used in non-iterative design
                                 % - leak2 is 10x19 matrix used in iterative design

%%% Fig. 3(a)
% The non-iterative design with 10 design points
Gdb1 = 12*ones(1,10)';     % Constant (dB)
Goptdb = leak1'\Gdb1;      % Solve optimal dB gains based on leakage
Gopt = 10.^(Goptdb/20);    % Convert to linear gain factors
Gwoptdb = gw*Goptdb;       % Gain at bandwidth wg
Gwopt = 10.^(Gwoptdb/20);  % Convert to linear gain factor

% Design filters with optimized gains
numsopt = zeros(3,10);  % 3 num coefficients for each 10 filters
densopt = zeros(3,10);  % 3 den coefficients for each 10 filters
for k = 1:10,
    [num,den] = pareq(Gopt(k), Gwopt(k), wg(k), bw(k)); % Design filters
    numsopt(:,k) = num;
    densopt(:,k) = den;
end
% Evaluate frequency responses
Nfreq = 2^11;  % Number of frequency points
w = logspace(log10(9),log10(22050),Nfreq);  % Log frequency points
H1opt = ones(Nfreq,10);   % Frequency response of individual filters
H1opttot = ones(Nfreq,1); % Frequency response of the cascade EQ
for k = 1:10
    H1opt(:,k) = freqz(numsopt(:,k),densopt(:,k),w,fs);
    H1opttot = H1opt(:,k) .* H1opttot;
end

% The iterative design with 19 design points
Gdb2 = zeros(19,1);
Gdb2(1:2:19) = Gdb1;
for k = 2:2:19
    Gdb2(k) = (Gdb2(k-1)+Gdb2(k+1))/2; % Interpolate target gains linearly b/w command gains
end
Goptdb = leak2'\Gdb2;      % Solve first estimate of dB gains based on leakage
Gopt = 10.^(Goptdb/20);    % Convert to linear gain factors

% Iterate once
leak22 = interactionMatrix(Gopt,gw,wg,wc2,bw); % Use previous gains
G2optdb = leak22'\Gdb2;     % Solve optimal dB gains based on leakage
G2opt = 10.^(G2optdb/20);   % Convert to linear gain factors
G2woptdb = gw*G2optdb;      % Gain at bandwidth wg
G2wopt = 10.^(G2woptdb/20); % Convert to linear gain factor

% Design filters with optimized gains
numsopt = zeros(3,10);  % 3 num coefficients for each 10 filters
densopt = zeros(3,10);  % 3 den coefficients for each 10 filters
for k = 1:10,
    [num,den] = pareq(G2opt(k), G2wopt(k), wg(k), bw(k)); % Design filters
    numsopt(:,k) = num;
    densopt(:,k) = den;
end
% Evaluate frequency responses
H2opt = ones(Nfreq,10);   % Frequency response of individual filters
H2opttot = ones(Nfreq,1); % Frequency response of the cascade EQ
for k = 1:10
    H2opt(:,k) = freqz(numsopt(:,k),densopt(:,k),w,fs);
    H2opttot = H2opt(:,k) .* H2opttot;
end

%%% Fig. 3(b)
% The non-iterative design with 10 design points
Gdb3 = [12 -12 -12 12 -12 -12 -12 12 -12 -12]';  % Oliver&Jot worst case (dB)
Goptdb = leak1'\Gdb3;     % Solve optimal dB gains based on leakage
Gopt = 10.^(Goptdb/20);   % Convert to linear gain factors
Gwoptdb = gw*Goptdb;      % Gain at bandwidth wg
Gwopt = 10.^(Gwoptdb/20); % Convert to linear gain factor

% Design filters with optimized gains
numsopt = zeros(3,10);  % 3 num coefficients for each 10 filters
densopt = zeros(3,10);  % 3 den coefficients for each 10 filters
for k = 1:10,
    [num,den] = pareq(Gopt(k), Gwopt(k), wg(k), bw(k)); % Design filters
    numsopt(:,k) = num;
    densopt(:,k) = den;
end
% Evaluate frequency responses
H3opt = ones(Nfreq,10);   % Frequency response of individual filters
H3opttot = ones(Nfreq,1); % Frequency response of the cascade EQ
for k = 1:10
    H3opt(:,k) = freqz(numsopt(:,k),densopt(:,k),w,fs);
    H3opttot = H3opt(:,k) .* H3opttot;
end

% The iterative design with 19 design points
Gdb4 = zeros(19,1);
Gdb4(1:2:19) = Gdb3;
for k = 2:2:19
    Gdb4(k) = (Gdb4(k-1)+Gdb4(k+1))/2; % Interpolate target gains linearly b/w command gains
end
Goptdb = leak2'\Gdb4;   % Solve first estimate of dB gains based on leakage
Gopt = 10.^(Goptdb/20); % Convert to linear gain factors
Gwoptdb = gw*Goptdb;    % Gain at bandwidth wg

% Iterate once
leak22 = interactionMatrix(Gopt,gw,wg,wc2,bw); % Use previous gains
G2optdb = leak22'\Gdb4;     % Solve optimal dB gains based on leakage
G2opt = 10.^(G2optdb/20);   % Convert to linear gain factors
G2woptdb = gw*G2optdb;      % Gain at bandwidth wg
G2wopt = 10.^(G2woptdb/20); % Convert to linear gain factor

% Design filters with optimized gains
numsopt = zeros(3,10);  % 3 num coefficients for each 10 filters
densopt = zeros(3,10);  % 3 den coefficients for each 10 filters
for k = 1:10,
    [num,den] = pareq(G2opt(k), G2wopt(k), wg(k), bw(k)); % Design filters
    numsopt(:,k) = num;
    densopt(:,k) = den;
end
% Evaluate frequency responses
H4opt = ones(Nfreq,10);   % Frequency response of individual filters
H4opttot = ones(Nfreq,1); % Frequency response of the cascade EQ
for k = 1:10
    H4opt(:,k) = freqz(numsopt(:,k),densopt(:,k),w,fs);
    H4opttot = H4opt(:,k) .* H4opttot;
end

% Plot responses for the proposed optimized design
figure(3); clf; subplot(2,1,1)
semilogx(w,db(H1opttot),'k','linewidth',1.5);  % Total response of non-iterative design
hold on
semilogx(w,db(H2opttot),'k','linewidth',3);    % Total response of iterative design
plot(fc1,Gdb1,'ko','linewidth',2)              % Command gains
semilogx([fc1(1) fc1(end)],13*ones(2,1),'k--') % Tolerance lines
semilogx([fc1(1) fc1(end)],11*ones(2,1),'k--')
set(gca,'fontname','Times','fontsize',16);
xlabel('(a)'); ylabel('Magnitude (dB)')
set(gca,'XTick',[10 30 100 300 1000 3000 10000],'YTick',-20:1:20)
axis([20 22050 9.5 13.5])

subplot(2,1,2)
semilogx(w,db(H3opttot),'k','linewidth',1.5); % Total response of non-iterative design
hold on
semilogx(w,db(H4opttot),'k','linewidth',3);   % Total response of iterative design
plot(fc1,Gdb3,'ko','linewidth',2)             % Command gains
semilogx([fc1(2) fc1(3)],-13*ones(2,1),'k--') % Tolerance lines
semilogx([fc1(2) fc1(3)],-11*ones(2,1),'k--')
semilogx([fc1(5) fc1(7)],-13*ones(2,1),'k--')
semilogx([fc1(5) fc1(7)],-11*ones(2,1),'k--')
semilogx([fc1(9) fc1(10)],-13*ones(2,1),'k--')
semilogx([fc1(9) fc1(10)],-11*ones(2,1),'k--')
set(gca,'fontname','Times','fontsize',16);
xlabel({'Frequency (Hz)';'(b)'}); ylabel('Magnitude (dB)')
set(gca,'XTick',[10 30 100 300 1000 3000 10000],'YTick',-20:1:20)
axis([20 22050 -13.5 -9.5])

%print('Proposed','-dpdf')

