Header imageAller au contenu principal

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):

image

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.

image

Commentaires