# length 27 (-13,...,13) linear phase FIR filter

param n := 14;
param pi := 4*atan(1);
param eps := 1.0e-4;
param sinc_w {k in -2*n+2..2*n-2} := 
    if k != 0 then (sin(2*pi*k*0.9)/(pi*k)) else (0.8);
param sinc_m {k in -2*n+2..2*n-2} := 
    if k != 0 then ( sin(2*pi*k*0.1)/(pi*k)+ sin(2*pi*k*0.8)/(pi*k)) else (0.8);
param sinc_t {k in -2*n+2..2*n-2} := 
    if k != 0 then (sin(2*pi*k*0.2)/(pi*k)) else (0.4);

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

minimize power_bnd: rho;

subject to passband: 
    ((hw[0]+hm[0]+ht[0]-1)^2 + 2*sum {k in 1..n-1} (hw[k]+hm[k]+ht[k])^2) 
    <= (eps);

subject to wooferband: 
    ( sum {k in -n+1..n-1, kk in -n+1..n-1} 
	           hw[abs(k)] * hw[abs(kk)] * sinc_w[k+kk]
	) /rho
    <= (0.6)*rho;

subject to midrangeband: 
    ( sum {k in -n+1..n-1, kk in -n+1..n-1} 
	           hm[abs(k)]*hm[abs(kk)] * sinc_m[k+kk]
	) /rho
    <= (0.4)*rho;

subject to tweeterband: 
    ( sum {k in -n+1..n-1, kk in -n+1..n-1} 
	           ht[abs(k)]*ht[abs(kk)] * sinc_t[k+kk]
	) /rho
    <= (0.6)*rho;

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

solve;

printf "%3d \n", card({-n+1..n-1}) > hw;
printf "%3d \n", card({-n+1..n-1}) > hm;
printf "%3d \n", card({-n+1..n-1}) > ht;
printf "%3d \n", card({-n+1..n-1}) > hwmt;
printf {k in -n+1..n-1}: "%10.6f \n", hw[abs(k)] > hw;
printf {k in -n+1..n-1}: "%10.6f \n", hm[abs(k)] > hm;
printf {k in -n+1..n-1}: "%10.6f \n", ht[abs(k)] > ht;
printf {k in -n+1..n-1}: "%10.6f \n", hw[abs(k)]+hm[abs(k)]+ht[abs(k)] > hwmt;

printf {nu in 0..1 by 1/1000}: "%7.4f %10.3e \n", 
    nu, 20*log10(abs(
    hw[0] + 2* sum {k in 1..n-1} (hw[k]*cos(-2*pi*k*nu))
    )) > 3aw.out;

printf {nu in 0..1 by 1/1000}: "%7.4f %10.3e \n", 
    nu, 20*log10(abs(
    hm[0] + 2* sum {k in 1..n-1} (hm[k]*cos(-2*pi*k*nu))
    )) > 3am.out;

printf {nu in 0..1 by 1/1000}: "%7.4f %10.3e \n", 
    nu, 20*log10(abs(
    ht[0] + 2* sum {k in 1..n-1} (ht[k]*cos(-2*pi*k*nu))
    )) > 3at.out;