function [B] = interactionMatrix_SSGE3(g, param)
% interactionMatrix_SSGE3.m
% 
% Compute the interaction matrix of a cascade graphic equalizer containing
% the leak factors to account for the band interaction when assigning
% filter gains. All filters are Orfanidis peak/notch filters with
% adjustable bandwidth gain and Nyquist gain.
% 
% Input parameters:
% g  = dB gain at which the leakage is determined
% param = struct containing the initialized parameters
% 
% Output:
% B = N by M matrix showing how the magnitude responses of the
% band filters leak to the design frequencies. N is determined from the
% length of the array wc (number of design frequencies) whereas M is 
% determined from the length of wg (number of filters)
%
% Created by Mario Antonelli, Moricone, Roma, Italia, 28 March 2022
% for the
% Acoustics Laboratory
% Dept. of Signal Processing and Acoustics
% Aalto University

inds = find(abs(g)>0.0);
vgains = param.vgains;
K = param.K;
G1 = param.G1;
gw = param.c; 
wc = param.wc;
wg = param.wg;
Z = exp(-1i*wc); % Creating z^-1
Z2 = Z.^2;       % Creating z^-2
B = zeros(numel(param.wc), numel(param.wg));
for i = 1:numel(inds)
    G  = 10^(g(inds(i))/20);
    GB = 10^(gw*g(inds(i))/20);
    w0 = wg(inds(i));
    k = interp1qr(vgains,K(:,inds(i)),abs(g(inds(i))));
    g1 = interp1qr(vgains,G1(:,inds(i)),abs(g(inds(i))));
    if g(inds(i))<0
        g1 = -db(g1);
        g1 = 10^(g1/20);
    end
    Dw = k*w0;
    [b, a] = peq(1, G, GB, w0, Dw, g1);
    H = (b(1) + b(2)*Z + b(3)*Z2)./(a(1) + a(2)*Z + a(3)*Z2);
    B(:,inds(i)) = 20*log10(abs(H))/g(inds(i));
end
end
