option loqo_options $loqo_options" convex";
# length 27 (-13,...,13) linear phase FIR filter

param n := 14;
param N := 1000;
param pi := 4*atan(1);
param eps := 1.0e-4;

set NU := {0..1 by 1/N};

param Sw2 {nu in NU} := (if 0.2<nu && nu<0.8 then 1 else 0);
param Sm2 {nu in NU} := 
    (if nu<0.1 || (0.4<nu && nu<0.6) || nu>0.9 then 1 else 0);
param St2 {nu in NU} := (if nu<0.3 || nu>0.7 then 1 else 0);

var rho >= 0;
var hw {0..n-1};
var hm {0..n-1};
var ht {0..n-1};

var Hw {nu in NU} =
    hw[0] + 2* sum {k in 1..n-1} (hw[k]*cos(-2*pi*k*nu));

var Hm {nu in NU} =
    hm[0] + 2* sum {k in 1..n-1} (hm[k]*cos(-2*pi*k*nu));

var Ht {nu in NU} =
    ht[0] + 2* sum {k in 1..n-1} (ht[k]*cos(-2*pi*k*nu));

minimize power_bnd: rho;

subject to passband: 
    sqrt(sum {nu in NU} (Hw[nu]+Hm[nu]+Ht[nu]-1)^2) 
    <= sqrt(N)*sqrt(eps);

subject to wooferband: 
    sqrt(sum {nu in NU} Hw[nu]^2 * Sw2[nu]) 
    <= sqrt(0.6*N)*rho;

subject to midrangeband: 
    sqrt(sum {nu in NU} Hm[nu]^2 * Sm2[nu]) 
    <= sqrt(0.4*N)*rho;

subject to tweeterband: 
    sqrt(sum {nu in NU} Ht[nu]^2 * St2[nu]) 
    <= sqrt(0.6*N)*rho;

let hw[0] := 2;
let hm[0] := 2;
let ht[0] := 2;

solve;

printf {nu in NU}: "%7.4f %10.3e \n", 
    nu, 20*log10(abs(Hw[nu])) > 3w.out;

printf {nu in NU}: "%7.4f %10.3e \n", 
    nu, 20*log10(abs(Hm[nu])) > 3m.out;

printf {nu in NU}: "%7.4f %10.3e \n", 
    nu, 20*log10(abs(Ht[nu])) > 3t.out;