% Fig2.m
%
% Plot the example of peaking EQ self-similarity. 
% This was done for the article V. Valimaki and J. Liski, "Accurate Cascade
% Graphic Equalizer," IEEE Signal Processing Letters, 2017
%
% Uses pareq.m
%
% Created by Juho Liski and Vesa Valimaki, Otaniemi, Espoo, Finland, 21 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 = 1000; % Center frequency for the filters
wc = 2*pi*fc/fs;  % Center frequency in radians
c = 0.3;  % Gain factor at bandwidth
bw = 0.21; % Bandwidth in radians

% Plotting parameters
fontsize = 16; 
style = {'k' 'k-.' 'k--' 'k'};
linewidth = [1 1 1 2.5];

% Plot self-similarity
figure(2);clf;
Nfreq = 2048;  % Number of frequency points
w = logspace(log10(0.2),log10(22050),Nfreq); % Log frequency points
for k = 3:7:24  % Step through the four plotted gains
    gain = 10.^(k/20);  % Interpret k as gain in dB
    gainw = 10.^(c*k/20);  % Gain of c*k dB at bandwidth
    [num,den] = pareq(gain, gainw, wc, bw); % Design filter
    [H,f] = freqz(num,den,w,fs); % Evaluate filter response at w
    semilogx(w,db(H)/k,style{(k+4)/7},'LineWidth',linewidth((k+4)/7));
    hold on
end
plot(500,0.3,'ks','LineWidth',2,'MarkerSize',15)  % Plot the three squares
plot(2000,0.3,'ks','LineWidth',2,'MarkerSize',15) % at own center frequency
plot(1000,1,'ks','LineWidth',2,'MarkerSize',15)   % and the neighboring ones
set(gca,'FontName','Times','FontSize',fontsize)
axis([100 10000 0.0 1.1])
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)')
set(gca,'XTick',[31.25 63.25 125 250 500 1000 2000 4000 8000 16000])
set(gca,'YTick',[0 0.3 0.7 1])

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