% Fig3.m
%
% Plot the example of peaking EQ self-similarity. 
% 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 pareq.m
%
% Created by Juho Liski, Otaniemi, Espoo, Finland, 6 May 2019
%
% Aalto University, Dept. of Signal Processing and Acoustics

fs  = 44.1e3;  % Sample rate
fc = 1000; % Center frequency for the filter
wc = 2*pi*fc/fs;  % Center frequency in radians
c = 0.4;  % Gain factor at bandwidth
bw = 0.0664; % Bandwidth in radians

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

% Plot self-similarity
figure(3);clf;
%subplot(1.6,1,1)
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(1000/(2^(1/3)),c,'ks','LineWidth',2,'MarkerSize',12)  % Plot the three squares
plot(1000*(2^(1/3)),c,'ks','LineWidth',2,'MarkerSize',12) % at own center frequency
plot(1000,1,'ks','LineWidth',2,'MarkerSize',12)   % and the neighboring ones
set(gca,'FontName','Times','FontSize',fontsize)
axis([250 4000 -0.01 1.1])
xlabel('Frequency (Hz)'); ylabel('Normalized magnitude (dB)')
set(gca,'XTick',[250 315 397 500 630 794 1000 1260 1587 2000 2520 3175 4000])
set(gca,'YTick',[0 0.4 0.7 1])

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