function [g, err_inf_norm] = ssge3_gains(gc,param)
% ssge3_gains.m
%
% Optimize the command gains for third-octave symmetric graphic equalizer
% as presented in "Sparse Graphic Equalizer Design" by M. Antonelli, 
% J. Liski, and V. Valimaki in IEEE SPL (2022).
% 
% Input parameters:
% gc = command gains in dB, size 31x1
% param = struct containing fixed design parameters for GEQ
% 
% Output:
% g  = optimized band filter gains in dB, size 31x1
% err_inf_norm = maximum absolute error at band filter center frequencies 
%
% Uses interp1qr.m, ODOMP.m, min0norm4PolySet.m, interactionMatrix_SSGE3.m 
%
% Created by Mario Antonelli, Moricone, Roma, Italia, 28 March 2022
% for the
% Acoustics Laboratory
% Dept. of Signal Processing and Acoustics
% Aalto University

t = interp1qr(param.wg',gc,param.wc');

%% Sparse solution computation by using Orthogonal Matching Pursuit) 
B = param.B;
W = param.W;
thr = 0.2;
[gopt, feasible, err] = ODOMP(numel(gc), B, t, thr, W);

%% If not feasible, Sparse solution computation by using Linear Programming
if not(feasible)
    lambda = 1000;
    gopt = min0norm4PolySet(B, t, lambda, W);
end

%% Second interaction matrix computation by using Optimal Sparse Solution
[B2] = interactionMatrix_SSGE3(gopt, param);

%% Least Square Optimization of the Gains by using the Sparse Interaction Matrix
% warning('off', 'all')
inds = abs(gopt)>0.0;
A = B2(:,inds);
x = (W*A)\(W*t);
g = zeros(size(B2,2),1);
g(inds) = x;
% g =(W*B2)\(W*t);
% warning('on', 'all')
e = B2*g-t;
% g = gopt;
% e = B*g-t;
err_inf_norm = max(abs(e(1:2:61)));