# FIR filter design experiment using the simplex method

A starter about symetrical (linear phase) finite impulse response digital filter design using the simplex method. Much slower than Parks/McClellan/Remez exchange, but packaged in GLPK. 75 lines of OCaml code to design an FIR filter, what less could you ask for ?

The paper to have in hands is the "Meteor" one by Steiglitz. In short, from a vector dot matrix multiplication, we get the frequency response. The matrix beeing a cosine table, the vector beeing the filter coefficients.

The simplex method beeing a relatively simple optimization procedure for linear problems.

### Building

You need ocaml-glpk, ocaml-gnuplot. Using findlib, I build it with:

ocamlfind opt -o simplex_fir -package glpk,gnuplot -linkpkg simplex_fir0.ml -ccopt \"-L/opt/local/lib\"

(* A first experiment with the simplex method for FIR filter design Philippe Strauss, philou at philou.ch, 2010 *) open Glpk module P = Gnuplot.Array type trigf = Cosine | Sine let foi = float_of_int let pi = 4. *. (atan 1.) let trig fct f harm = let om = 2. *. pi *. f *. (foi harm) in match fct with | Cosine -> cos om | Sine -> sin om (* _T for transposed, about the vector *) let mul_mat_vect_T a x = let lena_l = Array.length a in let lena_c = Array.length (a.(0)) in let lenx = Array.length x in assert (lena_c = lenx) ; let mul = Array.make lena_l 0. in for i = 0 to lena_l - 1 do let a_i = a.(i) in (* hoisting *) let accum = ref 0. in for j = 0 to lenx - 1 do accum := !accum +. x.(j) *. a_i.(j) done ; mul.(i) <- !accum ; done ; mul let init_trig_table fct nf ncol = let tbl = Array.make_matrix nf ncol 0. in let fstep = 0.5 /. (foi nf) in for i = 0 to nf - 1 do for j = 0 to ncol - 1 do tbl.(i).(j) <- trig fct (fstep *. (foi i)) j done ; done ; tbl let () = let ncoeffs = 151 in let ncoeffs2 = (ncoeffs - 1) / 2 in let nf = 1000 in let a = init_trig_table Cosine nf ncoeffs2 in let fstep = 0.5 /. (foi nf) in let target_mag_response f = if f < 0.15 then (0.99, 1.01) else if f < 0.18 then (0.0, 1.01) else (-0.00001, 0.00001) in let pbounds0 = Array.make nf (-.infinity, infinity) in let pbounds = Array.mapi (fun i x -> target_mag_response (fstep *. (foi i))) pbounds0 in let zcoeffs = Array.make ncoeffs2 1. in let xbounds = Array.make ncoeffs2 (-.1., 1.) in let lp = make_problem Maximize zcoeffs a pbounds xbounds in scale_problem lp ; use_presolver lp true ; simplex lp ; let coeffs = get_col_primals lp in for j = 0 to Array.length coeffs - 1 do Printf.printf "%4d = %.5f\n" j coeffs.(j) ; done ; let mag = mul_mat_vect_T a coeffs in let mag_db = Array.map (fun x -> 10. *. (if x = 0.0 then (-20.) else log10 (x ** 2.))) mag in let freq = Array.init nf (fun i -> foi i *. 0.5 /. foi nf) in let g = P.init ~xsize:800. ~ysize:400. Gnuplot.X in P.title g "Amplitude response" ; P.xlabel g "Freq" ; P.env g ~xlog:true 0.001 (0.5) (-120.) 10.; P.ylabel g "Amplitude [dB]" ; P.pen g 1 ; P.xy g freq mag_db ; P.close g

Here's a screenshot of the gnuplot invocation, showing the frequency response when no truncation of coefficients occurs (float32 is much enough accuracy here):

And for comparison, two FIR filters designed using the classical Remez exchange algorithm, one of 371 taps in red, the blue one beeing of 401 taps.

## Comments