% Figs45.m
%
% Plot the zigzag gain configuration with 31-band ACGE and APGE and their
% difference in dB.
%
% 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 interactionMatrix2.m and 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
Gdb = -12*ones(1,31); Gdb(1:2:31) = -Gdb(:,1);  % Zig-zag (dB)

% Design APGE and ACGE
[B1,A1,F,nums,dens,G0] = ACGE2APGE(Gdb);

% Evaluate and plot cascade structure
Nfreq = 2^11; w = logspace(log10(9),log10(22050),Nfreq);   % for plotting Nfreq = 2^11
H1 = ones(Nfreq,31); H1tot = ones(Nfreq,1);
for k = 1:31
    H1(:,k) = freqz(nums(:,k),dens(:,k),w,fs);
    H1tot = H1(:,k) .* H1tot;
end
H1tot = G0*H1tot;

% Evaluate and plot parallel structure:

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

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

H1partot = zeros(length(W),1);
H1par = zeros(length(W),NSEC);

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

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

figure(41); %subplot(1.6,1,1)
for k = 1:31
    semilogx(w,db(H1(:,k)),'linewidth',1.5);hold on
end
semilogx(w,db(G0)*ones(length(w),1),'linewidth',1.5,'Color',[0 0 0]);
semilogx(w,db(H1tot),'k','linewidth',2.5);
set(gca,'FontName','Times','FontSize',15); %grid on
xlabel({'Frequency (Hz)';'(a)'}); ylabel('Magnitude (dB)')
set(gca,'XTick',[30 100 300 1000 3000 10000]); axis([15 22050 -30 30])
set(gca,'YTick',-60:12:60)
hold off

figure(42); %subplot(1.6,1,1)
for k = 1:31
    semilogx(w,db(H1par(:,k)),'linewidth',1.5);hold on
end
semilogx(w,db(F)*ones(length(w),1),'linewidth',1.5,'Color',[0 0 0]);
semilogx(w,db(H1partot),'k','linewidth',2.5);
set(gca,'FontName','Times','FontSize',15); %grid on
xlabel({'Frequency (Hz)';'(b)'}); ylabel('Magnitude (dB)')
set(gca,'XTick',[30 100 300 1000 3000 10000]); axis([15 22050 -30 30])
set(gca,'YTick',-60:12:60)
hold off

%print(figure(41),'-dpdf','ACGE_zigzag')
%print(figure(42),'-dpdf','APGE_zigzag')

figure(5); %subplot(1.6,1,1)
semilogx(w,db(H1tot)-db(H1partot),'k','linewidth',2);
set(gca,'FontName','Times','FontSize',15); %grid on
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)')
set(gca,'XTick',[30 100 300 1000 3000 10000]); xlim([15 22050]); ylim([-1.2 1.2]*10^(-9))

%print(figure(5),'-dpdf','difference')
