param n := 4;
param N := 3000;
param pi := 4*atan(1);

set NU_P := {0.0..0.25 by 1/N} union {0.75..1.0 by 1/N};
set NU_S := {0.25..0.75 by 1/N};
set NU   := NU_P union NU_S;

var rho >= 0;
var h_real {0..n-1};
var h_imag {1..n-1};

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

minimize power_bnd: rho;

subject to bnd_def: sqrt(sum {nu in NU_S} H[nu]^2) <= sqrt(N)*rho;

subject to DC_gain: H[0] = 1;

let h_real[0] := 1;

solve;

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

display h_real, h_imag;