% Fig1.m
%
% Plot the amplitude response of an octave EQ based on Oliver-Jot design
% with two different command gain settings. Also, the amplitude responses
% of the individual filters are plotted.
% 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, 23 December 2016
%
% Aalto University, Dept. of Signal Processing and Acoustics

fs  = 44.1e3;  % Sample rate
fc = 16000./(2.^(9:-1:0)); % Exact log center frequencies for filters
wg = 2*pi*fc/fs;  % Command gain frequencies in radians
wc = 2*pi*fc/fs;  % Center frequencies in radians
gw = 0.5; % Bandwidth is defined as midpoint gain (sqrt(Gdb))
bw = 2*pi*[32.76 65.51 131.0 262.1 524.1 1048 2096 4193 7212 10730]/fs; % Bandwidths in radians

%%% Fig. 1(a)
Gdb1 = 12*ones(1,10)';  % Constant (dB)
leak = interactionMatrix(10.^(10/20)*ones(1,10),gw,wg,wc,bw); % Calculate interaction matrix using 10 dB as prototype gain
Goptdb = leak'\Gdb1;  % Solve optimal dB gains based on leakage
Gopt = 10.^(Goptdb/20);    % Convert to linear gain factors
Gwoptdb = gw*Goptdb;       % Gain at crossover bandwidth
Gwopt = 10.^(Gwoptdb/20);  % Convert to linear gain factors

nums = zeros(3,10);  % 3 num coefficients for each 10 filters
dens = 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 the filters
    nums(:,k) = num;
    dens(:,k) = den;
end

% Evaluate frequency response
Nfreq = 2^11;  % Number of frequency points
w = logspace(log10(9),log10(22050),Nfreq);  % Log frequency points
H1 = ones(Nfreq,10);   % Frequency response of individual filters
H1tot = ones(Nfreq,1); % Frequency response of the cascade EQ
for k = 1:10
    H1(:,k) = freqz(nums(:,k),dens(:,k),w,fs);
    H1tot = H1(:,k) .* H1tot;
end

%%% Fig. 1(b)
Gdb2 = [12 -12 -12 12 -12 -12 -12 12 -12 -12]';  % Oliver&Jot worst case (dB)
Goptdb = leak'\Gdb2;  % Solve optimal dB gains based on leakage using the same interaction matrix as above
Gopt = 10.^(Goptdb/20);    % Convert to linear gain factors
Gwoptdb = gw*Goptdb;       % Gain at crossover bandwidth
Gwopt = 10.^(Gwoptdb/20);  % Convert to linear gain factors

nums = zeros(3,10);  % 3 num coefficients for each 10 filters
dens = 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 the filters
    nums(:,k) = num;
    dens(:,k) = den;
end

% Evaluate frequency response
H2 = ones(Nfreq,10);   % Frequency response of individual filters
H2tot = ones(Nfreq,1); % Frequency response of the cascade EQ
for k = 1:10
    H2(:,k) = freqz(nums(:,k),dens(:,k),w,fs);
    H2tot = H2(:,k) .* H2tot;
end

% Plot responses for the Oliver&Jot design
color = [0.5 0.5 0.5]; % Color for individual filters
figure(1); clf; subplot(2,1,1)
for k = 1:10 % Individual filters
    semilogx(w,db(H1(:,k)),'linewidth',1,'Color',color,'Linestyle','-');hold on
end
semilogx(w,db(H1tot),'k','linewidth',3);    % Total response
plot(fc,Gdb1,'ko','linewidth',2)            % Command gains
set(gca,'fontname','Times','fontsize',15);
xlabel('(a)'); ylabel('Magnitude (dB)')
set(gca,'XTick',[30 100 300 1000 3000 10000],'YTick',[0 6 12])
axis([20 22050 0 14])

subplot(2,1,2)
for k = 1:10 % Individual filters
    semilogx(w,db(H2(:,k)),'linewidth',1,'Color',color,'Linestyle','-');hold on
end
semilogx(w,db(H2tot),'k','linewidth',3);    % Total response
plot(fc,Gdb2,'ko','linewidth',2)            % Command gains
set(gca,'fontname','Times','fontsize',15);
xlabel({'Frequency (Hz)';'(b)'});ylabel('Magnitude (dB)')
set(gca,'XTick',[10 30 100 300 1000 3000 10000],'YTick',[-12 0 12 24])
axis([20 22050 -19 27])

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

