function [numsAPGE,densAPGE,F,numsACGE,densACGE,G0] = ACGE2APGE(Gdb,plot_mode)
% ACGE2APGE.m
%
% Design an octave or third-octave cascade EQ and convert it to the delayed
% parallel form based on 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.
%
% Input parameters:
% Gdb  = command gains in dB, size 1x10 (octave) or 1x31 (third-octave)
% plot_mode = determine if the responses are plotted (1) or not (0)
% 
% Output:
% numsAPGE = numerator parts of the parallel EQ filters
% densAPGE = denominator parts of parallel EQ filters
% F = direct gain for APGE
% numsACGE = numerator parts of the cascade EQ filters
% densACGE = denominator parts of the cascade EQ filters
% G0 = direct gain for ACGE
%
% Uses interactionMatrix2.m & pareq.m
%
% Usage:
% [B1,A1,F,B2,A2,G0] = ACGE2APGE(12*(2*rand(1,31)-1),0)     % produces third-octave EQ with random gains, no figure
% [B1,A1,F,B2,A2,G0] = ACGE2APGE(12*sign(2*rand(1,10)-1),1) % produces octave EQ with random +/-12 dB gains with figure
%
% Created by Juho Liski, Otaniemi, Espoo, Finland, 6 May 2019
%
% Aalto University, Dept. of Signal Processing and Acoustics

% If plot_mode isn't specified, do not plot
if (nargin<2), plot_mode = 0; end

% Input length must be 10 or 31:
if length(Gdb) ~= 10 && length(Gdb) ~= 31
    error('Gdb length must be either 10 (octave) or 31 (third-octave)!')
end

% Determine mode, aka octave (0) or third-octave (1):
design_ver = rem(length(Gdb),10);

%%% Sampling frequency
fs = 44.1e3;

% More example command gain configurations:

%%% Third-octave
%Gdb = 12*ones(1,31);  % Constant (dB)
%Gdb = -12*ones(1,31); Gdb(1:2:31) = -Gdb(:,1); % Zigzag
%Gdb = [12 -12 12 12 -12 12 -12 -12 -12 12 -12 12 12 12 -12 12 -12 -12 12 12 12 -12 -12 12 12 -12 12 -12 -12 12 12]; % Extreme case modified from Oliver & Jot
%Gdb = [12 12 11 9 6 3 1 0 0 6 6 12 6 6 -12 12 -12 -12 -12 -12 0 0 0 0 -3 -6 -9 -12 0 0 0]; % Nonextreme
%Gdb = 12*(2*rand(1,31)-1);  % Random (dB)
%Gdb = 12*sign(2*rand(1,31)-1); % Random +/-12 dB

%%% Octave
%Gdb = 12*ones(1,10);  % Constant (dB)
%Gdb = -12*ones(1,10); Gdb(1:2:10) = -Gdb(:,1); % Zigzag
%Gdb = [12 -12 12 12 12 -12 -12 -12 12 12]; % Extreme case from Liski & Valimaki
%Gdb = 12*(2*rand(1,10)-1);  % Random (dB)
%Gdb = 12*sign(2*rand(1,10)-1); % Random +/-12 dB

%%%%%%%%%%%%%%%%%%%%
%%% Desing ACGE: %%%
%%%%%%%%%%%%%%%%%%%%

if design_ver
    fc1 = [1000./(2.^((1/3)*(17:-1:0))) 1000.*(2.^((1/3)*(1:13)))];
    c = 0.4; % Gain factor at bandwidth
else
    fc1 = 16000./(2.^(9:-1:0));
    c = 0.3; % Gain factor at bandwidth
end
fc2 = zeros(1,2*length(Gdb)-1); fc2(1:2:2*length(Gdb)-1) = fc1;
for k = 2:2:2*length(Gdb)-1
    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
wc = 2*pi*fc2/fs;  % Center frequencies in radians for iterative design with extra points

if design_ver
    bw = (nthroot(2,3)-1/nthroot(2,3))*wg;
    bw(26) = 0.9615*bw(26); bw(27) = 0.939*bw(27);
    bw(28) = 0.905*bw(28); bw(29) = 0.851*bw(29);
    bw(30) = 0.7627*bw(30); bw(31) = 0.593*bw(31);
else
    bw = 1.5*wg; 
    bw(8) = 0.93*bw(8); bw(9) = 0.78*bw(9); bw(10) = 0.76*wg(10);
end

leak = interactionMatrix2(10^(17/20)*ones(1,length(Gdb)),c,wg,wc,bw); % Estimate leakage b/w bands
Gdb2 = zeros(2*length(Gdb)-1,1);
Gdb2(1:2:2*length(Gdb)-1) = Gdb;
for k = 2:2:2*length(Gdb)-1
    Gdb2(k) = (Gdb2(k-1)+Gdb2(k+1))/2; % Interpolate target gains linearly b/w command gains
end

Goptdb = leak'\Gdb2;      % Solve first estimate of dB gains based on leakage
Gopt = 10.^(Goptdb/20);    % Convert to linear gain factors

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

% Design filters with optimized gains
nums = zeros(3,length(Gdb));  % 3 num coefficients for each 31 filters
dens = zeros(3,length(Gdb));  % 3 den coefficients for each 31 filters
G0 = 1;
for k = 1:length(Gdb)
    [num,den] = pareq(G2opt(k), G2wopt(k), wg(k), bw(k)); % Design filters
    G0 = G0*num(1);
    nums(:,k) = num/num(1);
    dens(:,k) = den;
end

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

if plot_mode
    %%% Plot cascade structure
    figure(1); subplot(2,1,1)
    semilogx(fc1,Gdb,'ro','linewidth',2,'MarkerSize',8); hold on
    for k = 1:length(Gdb)
        semilogx(w,db(H1(:,k)),'linewidth',1.5);
    end
    semilogx(w,db(G0)*ones(length(w),1),'linewidth',1.5);
    semilogx(w,db(H1tot),'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]); axis([15 22050 -30 30])
    hold off
end

numsACGE = nums;
densACGE = dens;

%%%%%%%%%%%%%%%%%%%%%%%%
%%% Convert to APGE: %%%
%%%%%%%%%%%%%%%%%%%%%%%%

B1 = zeros(length(Gdb),2);
%%% Solve poles; initialize
p = zeros(2*length(Gdb),1);
%%% Solve poles using quadratic formula
for sec = 1:length(Gdb) % go through all sections (2 poles for each section)
    a = dens(:,sec);
    p(2*sec-1,1) = (-a(2)+sqrt(a(2)^2-4*a(3)))/2; 
    p(2*sec,1) = (-a(2)-sqrt(a(2)^2-4*a(3)))/2;
end

%%% Solve 1st-order residues
r_nums = ones(2*length(Gdb),length(Gdb));
r_dens = ones(2*length(Gdb),length(Gdb));
for k=1:length(Gdb)
    % Calculating "polyval" manually:
    for m=2:3
        r_nums(:,k) = p .* r_nums(:,k) + nums(m,k); % All p's evaluated at numerators
        r_dens(:,k) = p .* r_dens(:,k) + dens(m,k); % All p's evaluated at denominators
    end
    % Taking the pole cancellation in to account:
    r_dens(2*k-1,k) = p(2*k-1) - p(2*k);
    r_dens(2*k,k) = p(2*k) - p(2*k-1);
end
rNum = prod(r_nums,2);
rDen = prod(r_dens,2);
r1 = (G0*rNum)./rDen;

%%% Direct gain
F = G0; % Direct gains are obtained when z is set to inf

%%% Parallel sections combined into second-order
for k = 1:length(Gdb)
    B1(k,:) = [r1(2*k-1)+r1(2*k) -r1(2*k-1)*p(2*k)-r1(2*k)*p(2*k-1)];
end
A1 = dens';

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

%%% Evaluation 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(k,1)*Z+B1(k,2)*Z2)./(A1(k,1)+A1(k,2)*Z+A1(k,3)*Z2);
    H1partot = H1partot + H1par(:,k);
end

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

if plot_mode
    subplot(2,1,2)
    semilogx(fc1,Gdb,'ro','linewidth',2,'MarkerSize',8); hold on
    for k = 1:length(Gdb)
        semilogx(w,db(H1par(:,k)),'linewidth',1.5);
    end
    semilogx(w,db(F)*ones(length(w),1),'linewidth',1.5);
    semilogx(w,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]); axis([15 22050 -30 30])
    hold off
end

numsAPGE = B1';
densAPGE = A1';
