(* Digital filters *)
(* Generating optimal FIR filters
This code is intended as introductory and hence the resulting
filters aren't quite optimal. More can be done -- which is left
as a homework.
*)
(*{{{ Background on digital filters *)
(*
Introduction to Digital Filters with Audio Applications
JULIUS O. SMITH III
Center for Computer Research in Music and Acoustics (CCRMA)
https://ccrma.stanford.edu/~jos/filters/filters.html
Also:
http://www.diku.dk/~jda/biosignal/lab07.pdf
*)
(*}}}*)
(*{{{ General, obviously correct FIR filter *)
(* filter b x
FIR filter the signal x[i] according to the filter with coefficients b[k]
*)
let filter : float array -> float array -> float array = fun b x ->
let m = Array.length b in
let y i =
if i < m-1 then x.(i) else
let sum = ref 0.0 in
for k = 0 to m-1 do
sum := !sum +. b.(k) *. x.(i-k)
done;
!sum
in
Array.init (Array.length x) y (* essentially, for-loop *)
;;
(*{{{ What to do at the border, for x[-i]? *)
(*
There are several options:
-- start filtering at the k'th sample (what we do now);
-- set them to zero;
-- set them to x.(0);
-- apply a smoothing window (e.g., Hamming window)
EXERCISE: try one of these other options for boundary processing
*)
(*}}}*)
(*{{{ Testing filter *)
(* impulses i j n
n-sample signal with an impulse at i and j (i and j may be the same)
*)
let impulses : int -> int -> int -> float array = fun i j n ->
let x = Array.make n 0.0 in
x.(i) <- 1.0;
x.(j) <- 1.0;
x
;;
let [|1.; 0.; 0.2; 0.; 0.|] =
filter [|0.5; 0.3; 0.2|] (impulses 0 0 5);;
(* The impulse response of a FIR filter is its coefficients *)
let [|0.; 0.; 0.5; 0.3; 0.2; 0.; 0.|] =
filter [|0.5; 0.3; 0.2|] (impulses 2 2 7);;
let [|1.; 0.; 0.7; 0.3; 0.2; 0.; 0.|] =
filter [|0.5; 0.3; 0.2|] (impulses 0 2 7);;
(*}}}*)
(*}}}*)
(*{{{ Specializing to the known filter order *)
(*
Specializing the filter when the length of 'b' is known statically.
The length of 'b' minus one is called the order of the filter.
*)
(* Start with the signature:
-- what is the type of array about which nothing is known statically?
-- what should be the type of an array with the statically known length?
*)
(*{{{ Answer *)
type spine_arr = float code array
type dyn_arr = float array code
;;
(*}}}*)
(*{{{ Attempt 1 *)
(*
let filter_spine : spine_arr -> (float array -> float array) code =
fun b -> .
.~(let m = Array.length b in
..)>.
;;
What is the problem?
*)
(*}}}*)
(*{{{ Attempt 2 *)
(*
let filter_spine : spine_arr -> (float array -> float array) code =
fun b -> .
.~(let m = Array.length b in
..)>.
;;
*)
(*}}}*)
(*{{{ Successful attempt *)
let filter_spine : spine_arr -> (float array -> float array) code =
fun b -> .
.~(let m = Array.length b in
.. in
for k = 0 to m-1 do
sum := .< .~(!sum) +. .~(b.(k)) *. x.(i-k) >.
done;
!sum)
in
Array.init (Array.length x) y (* essentially, for-loop *)
>.)>.
;;
(*
. .~(filter_spine [| b1; b2; b3 |])>.;;
MetaOCaml is typed
*)
let _ = . .~(filter_spine [| ..; ..; .. |])>.;;
(*
- : (float -> float -> float -> float array -> float array) code = .<
fun b1_1 ->
fun b2_2 ->
fun b3_3 ->
fun x_4 ->
let y_6 i_5 =
if i_5 < (3 - 1)
then x_4.(i_5)
else
((0.0 +. (b1_1 *. (x_4.(i_5 - 0)))) +. (b2_2 *. (x_4.(i_5 - 1))))
+. (b3_3 *. (x_4.(i_5 - 2)))
in
Array.init (Array.length x_4) y_6>.
*)
(*}}}*)
(*
Application: request the coefficients from the user (GUI)
and generate a filter.
The loop is inlined.
*)
(*}}}*)
(*{{{ Specialization to known coefficients *)
(* What if the coefficients are known statically?
This seems to be the usual case.
Again, start with the original filter, write the signature,
and then add the staging annotations where the compiler tells us to.
*)
(*{{{ First, change the signature *)
(*
let filter_staged : float array -> (float array -> float array) code =
fun b -> .
.~(let m = Array.length b in
.. in
for k = 0 to m-1 do
sum := .< .~(!sum) +. .~(b.(k)) *. x.(i-k) >.
done;
!sum)
in
Array.init (Array.length x) y (* essentially, for-loop *)
>.)>.
;;
*)
(*}}}*)
(*{{{ Then fix the typing problems *)
let filter_staged : float array -> (float array -> float array) code =
fun b -> .
.~(let m = Array.length b in
.. in
for k = 0 to m-1 do
sum := .< .~(!sum) +. b.(k) *. x.(i-k) >.
done;
!sum)
in
Array.init (Array.length x) y (* essentially, for-loop *)
>.)>.
;;
let _ = filter_staged [| 0.5; 0.3; 0.2 |];;
(* What is the problem? How to fix it? *)
(*
- : (float array -> float array) code = .<
fun x_7 ->
let y_9 i_8 =
if i_8 < (3 - 1)
then x_7.(i_8)
else
((0.0 +. (((* CSP b *).(0)) *. (x_7.(i_8 - 0)))) +.
(((* CSP b *).(1)) *. (x_7.(i_8 - 1))))
+. (((* CSP b *).(2)) *. (x_7.(i_8 - 2)))
in
Array.init (Array.length x_7) y_9>.
*)
(*}}}*)
(*{{{ Successful specialization *)
let filter_staged : float array -> (float array -> float array) code =
fun b -> .
.~(let m = Array.length b in
.. in
for k = 0 to m-1 do
sum := .< .~(!sum) +. .~(let bk = b.(k) in ..) *. x.(i-k) >.
done;
!sum)
in
Array.init (Array.length x) y (* essentially, for-loop *)
>.)>.
;;
let _ = filter_staged [| 0.5; 0.3; 0.2 |];;
(*
- : (float array -> float array) code = .<
fun x_10 ->
let y_12 i_11 =
if i_11 < (3 - 1)
then x_10.(i_11)
else
((0.0 +. (0.5 *. (x_10.(i_11 - 0)))) +. (0.3 *. (x_10.(i_11 - 1)))) +.
(0.2 *. (x_10.(i_11 - 2)))
in
Array.init (Array.length x_10) y_12>.
*)
(* Or, we may lift the let statement up ... *)
let filter_staged : float array -> (float array -> float array) code =
fun b -> .
.~(let m = Array.length b in
.. in
for k = 0 to m-1 do
let bk = b.(k) in
sum := .< .~(!sum) +. bk *. x.(i-k) >.
done;
!sum)
in
Array.init (Array.length x) y (* essentially, for-loop *)
>.)>.
;;
let _ = filter_staged [| 0.5; 0.3; 0.2 |];;
(* Lesson: let-insertion.
Let-insertion is useful on many occasions. Sometimes we can do
it automatically (generating code with let-insertion).
That is beyond the scope of the tutorial though...
It requires effects.
*)
(*}}}*)
(*{{{ How not to unroll the loop *)
(* The two earlier specializations, filter_spine and filter_staged,
completely unrolled the loop in the |y| functions, which computed
the inner product (convolution) of the coefficients with the window
into the signal. The unrolling removes the overhead of the iteration,
at the expense of the code size. If there are many coefficients,
the code bloat may turn out detrimental. We now attempt to design
a `smarter' filter_spine, which fully unrolls the iteration if the
number of coefficients is less than some threshold.
*)
(*
let filter_spine_adaptive ?(threshold=3) :
spine_arr -> (float array -> float array) code = fun b ->
let m = Array.length b in
if m <= threshold then
filter_spine b
else
.
let y i =
if i < m-1 then x.(i) else
let sum = ref 0.0 in
for k = 0 to m-1 do
sum := !sum +. b.(k) *. x.(i-k)
done;
!sum
in
Array.init (Array.length x) y>.
;;
Characters 328-333:
sum := !sum +. b.(k) *. x.(i-k)
^^^^^
Error: This expression has type float code
but an expression was expected of type float
*)
(* Another sort of an eta-expansion;
cf. the transition from spower to spowern
*)
let dyn_arr_of_spine_arr : spine_arr -> dyn_arr = fun xa ->
. .<.~x :: .~acc>.) xa .<[]>.)>.
;;
let _ = dyn_arr_of_spine_arr [| .<1.>.; .<2.>.; .<3.>. |];;
(*
- : dyn_arr = ..
*)
let filter_spine_adaptive ?(threshold=3) :
spine_arr -> (float array -> float array) code = fun b ->
let m = Array.length b in
if m <= threshold then
filter_spine b
else
.
let y i =
if i < m-1 then x.(i) else
let sum = ref 0.0 in
for k = 0 to m-1 do
sum := !sum +. bdyn.(k) *. x.(i-k)
done;
!sum
in
Array.init (Array.length x) y>.
;;
(*
val filter_spine_adaptive :
?threshold:int -> spine_arr -> (float array -> float array) code =
*)
let _ = .
.~(filter_spine_adaptive [| ..; ..; .. |])>.;;
(* as before *)
let _ = .
.~(filter_spine_adaptive [| ..; ..; ..; .. |])>.;;
(*
- : (float -> float -> float -> float -> float array -> float array) code =
.<
fun b1_28 ->
fun b2_29 ->
fun b3_30 ->
fun b4_31 ->
let bdyn_32 = Array.of_list [b1_28; b2_29; b3_30; b4_31] in
fun x_33 ->
let y_37 i_34 =
if i_34 < (4 - 1)
then x_33.(i_34)
else
(let sum_35 = Pervasives.ref 0.0 in
for k_36 = 0 to 4 - 1 do
sum_35 :=
((! sum_35) +. ((bdyn_32.(k_36)) *. (x_33.(i_34 - k_36))))
done;
! sum_35)
in
Array.init (Array.length x_33) y_37>.
*)
let _ = .
.~(filter_spine_adaptive ~threshold:2 [| ..; ..; .. |])>.;;
(*}}}*)
(*{{{ Echo filter and its challenges *)
(*
Echo filter:
y_i = 1 /(1+alpha) * x_i + alpha/(1+alpha) * x[i-P]
delay (corresponding to P) should be about 0.2 secs, alpha is about 0.95
*)
let alpha = 0.95 in let alpha1 = 1. +. alpha in
filter_staged [| 1. /. alpha1; 0.; 0.;0.;0.;0.;0.; alpha /. alpha1 |];;
(* Anyone sees the problem?
For echo to be audible, the delay should be at least 0.2 secs, which
translates for many -- thousands -- zero samples.
Can we count on an optimizing compiler to eliminate multiplications
by zero? (cf. William Kahan) So, we need domain knowledge.
In general, a lot of effective optimizations depend on domain knowledge
and so cannot be realistically implemented in a general-purpose
optimizing compiler.
*)
(*}}}*)
(*}}}*)
(* Problems to think:
-- filtering inplace
-- IIR (recursive) filters
-- filter design (echo below is the trivial example)
*)
(*{{{ A bit of fun *)
let sampling_rate = 8000. (* Default sampling rate of /dev/dsp: 8Khz *)
let la = 440. (* A of the first octave *)
let pi = 4. *. atan2 1. 1.;;
let a12 =
let duration = 3.0 in (* seconds *)
let coeff1 = 2. *. pi *. la /. sampling_rate in
let coeff2 = 2. *. pi *. (2.4 *. la) /. sampling_rate in
let n = int_of_float (duration *. sampling_rate) in
Array.init n (fun i ->
sin ( coeff1 *. float_of_int i) +. 0.7 *. cos ( coeff2 *. float_of_int i))
(* Frequency modulation *)
let siren =
let duration = 3.0 in (* seconds *)
let coeff1 = 2. *. pi *. 2.5 /. sampling_rate in
let coeff2 = 2. *. pi *. (2.4 *. la) /. sampling_rate in
let n = int_of_float (duration *. sampling_rate) in
Array.init n (fun i ->
let fr = coeff2 *.
(1. +. 0.005 *. sin ( coeff1 *. float_of_int i)) in
sin ( fr *. float_of_int i))
;;
let play : float array -> unit = fun xs ->
let attenuation = 120. in
let maxx = Array.fold_left (fun a x -> max (abs_float x) a) 0.0 xs in
let _ = Printf.printf "maxx is %g\n" maxx in
let cout = open_out "/dev/dsp" in
Array.iter (fun x ->
output_byte cout ( truncate (x *. attenuation /. maxx) + 128 )) xs;
close_out cout
;;
(*
Echo filter
y_i = 1 /(1+alpha) x_i + alpha/(1+alpha) * x[i-P]
delay should be about 0.2 secs, alpha is about 0.95
*)
let echo : float -> float -> float array = fun alpha delay ->
let p = truncate (delay *. sampling_rate) in
let b = Array.make p 0.0 in
b.(0) <- 1. /. (1. +. alpha);
b.(p-1) <- alpha /. (1. +. alpha);
b
;;
(*
play a12;;
play (filter [|1.1; 0.5; -0.6|] a12);; (* low-pass *)
play (filter [|-1.0; 0.5; 0.5|] a12);; (* high-pass *)
play siren;;
play (filter [|1.1; 0.5; -0.6|] siren);; (* low-pass *)
play (filter [|-0.9; 0.5; 0.4|] siren);; (* high-pass *)
play (filter (echo 0.95 0.2) siren);;
*)
(*}}}*)