# This is the problem whose solution is shown in Figure 10
# of 
# "FIR Filter Design via Spectral Factorization and Convex Optimization"
# S.P. Wu, S. Boyd, and L. Vandenberghe

# In this variant we fix delta and minimize alpha
# The problem is formulated as an SOCP in the original
# complex variables w[0],...w[n-1].

param n := 12;
param N := 300;
param pi := 4*atan(1);
param delta := 0.11;
param d_over_lambda := 0.45;
param omega_0 := -2*pi*d_over_lambda*cos(0);
param omega_b := -2*pi*d_over_lambda*cos(pi/6);
param omega_s := -2*pi*d_over_lambda*cos(pi/4);
param omega_pi:= -2*pi*d_over_lambda*cos(pi);

set OMEGA_P := {omega_0..omega_b  by omega_pi/N};
set OMEGA_S := {omega_s..omega_pi by omega_pi/N};
set OMEGA_I := {omega_b..omega_s  by omega_pi/N};
set OMEGA := OMEGA_P union OMEGA_S union OMEGA_I;

var alpha >= 0;
var w_real {0..n-1};
var w_imag {0..n-1};

var H_real {o in OMEGA} = 
    sum {k in 0..n-1} ( w_real[k]*cos(k*o) + w_imag[k]*sin(k*o));
var H_imag {o in OMEGA} = 
    sum {k in 0..n-1} (-w_real[k]*sin(k*o) + w_imag[k]*cos(k*o));
var R {o in OMEGA} =
    H_real[o]^2 + H_imag[o]^2;

minimize stop_band_signal_bnd: delta;

subject to ripple_up_bnds {o in OMEGA_P}:            sqrt(R[o]) <= alpha;

subject to ripple_lo_bnds {o in OMEGA_P}: 1/alpha <= sqrt(R[o]);

subject to stop_bnd_def {o in OMEGA_S}: R[o] <= delta^2;

let alpha := 1;
#let {j in 0..n-1} w_real[j] := 1;
let w_real[0] := 1;
display {o in OMEGA} R[o];

solve;

display (delta);

printf {o in OMEGA}: "%7.4f %10.3e \n%7.4f %10.3e \n", 
    acos(-o/(2*pi*d_over_lambda)), 
    max(0,log10(sqrt( H_real[o]^2 + H_imag[o]^2)) + 3),
    -acos(-o/(2*pi*d_over_lambda)), 
    max(0,log10(sqrt( R[o])) + 3)
	> antenna.out;

printf "\n" >> antenna.out;
printf "%7.4f %10.3e \n", 0, (-0/20)+3 >> antenna.out;
printf "%7.4f %10.3e \n", 0, (-10/20)+3 >> antenna.out;
printf "%7.4f %10.3e \n", 0, (-20/20)+3 >> antenna.out;
printf "%7.4f %10.3e \n", 0, (-30/20)+3 >> antenna.out;