Signal

src/owl/dune

@@ -1,3 +1,4 @@
+
 (copy_files# ext/*)
 
 (copy_files# core/*)
@@ -42,6 +43,8 @@
 
 (copy_files# working/*)
 
+(copy_files# signal/*)
+
 (library
  (name owl)
  (public_name owl)

src/owl/owl.ml

@@ -41,6 +41,7 @@ module Graph = Owl_graph
 module Nlp = Owl_nlp
 module Log = Owl_log
 module Computation = Owl_computation
+module Signal = Owl_signal
 
 (* backend modules *)
 

src/owl/signal/owl_signal.ml

@@ -0,0 +1,64 @@
+(*
+ * OWL - OCaml Scientific and Engineering Computing
+ * Copyright (c) 2016-2020 Liang Wang <[email protected]>
+ *)
+
+(** Signal processing helpers *)
+
+open Owl_dense
+
+let hamming m =
+  let open Ndarray in
+  let range = (D.sequential [|m|] |> D.mul_scalar) (Owl_const.pi *. 2.0 /. (Int.to_float (m - 1))) in
+  let inter = D.cos range in
+  let mulv = D.mul_scalar inter (-0.46) in
+  D.add_scalar mulv 0.54
+
+let hann m =
+  let open Ndarray in
+  let range = (D.sequential [|m|] |> D.mul_scalar) (Owl_const.pi *. 2.0 /. (Int.to_float (m - 1))) in
+  let inter = D.cos range in
+  let mulv = D.mul_scalar inter (-0.5) in
+  D.add_scalar mulv 0.5
+
+let blackman m =
+  let open Ndarray in
+  let pi = Owl_const.pi in
+  let tpi=2.0 *. pi in
+  let fpi=4.0 *. pi in
+  let range1 = D.linspace 0. tpi m in
+  let range2 = D.linspace 0. fpi m in
+  let inter1 = D.cos range1 in
+  let inter2 = D.cos range2 in
+  let mulv1 = D.mul_scalar inter1 (-0.5) in
+  let mulv2 = D.mul_scalar inter2 (0.08) in
+  let mulvf = D.add mulv1 mulv2 in
+  let ans = D.add_scalar mulvf 0.42 in
+  ans
+
+let resize r x =                    (*zero pad the array to 2*x length*)
+  let open Ndarray in
+  let z = Array.make (2*x-(Array.length r)) 0. in
+  let y= Array.append r z in
+  D.of_array y [|2*x|] |> Generic.cast_d2z
+
+
+let dtft r x =                  (*dtft for upper unit circle (i.e whole if false)*)
+  let open Ndarray in
+  let a = resize r x in
+  let b = Owl_fft.D.fft a in
+  Z.get_slice [[0;x-1]] b
+
+let dtftw r x =              (*dtft for full circle (i.e whole is true)*)
+  let a = resize r x in
+  Owl_fft.D.fft a
+
+let freqz ?(n=512) ?(whole=false) b a = (*b represents numerator array while a represent denominator array*)
+  if whole then
+    let x = dtftw a n in
+    let y = dtftw b n in
+    (Ndarray.D.linspace (-1.0 *. Owl_const.pi) Owl_const.pi (n+1), Ndarray.Z.div y x)
+  else	
+    let x = dtft a n in
+    let y = dtft b n in
+    (Ndarray.D.linspace 0. Owl_const.pi (n+1), Ndarray.Z.div y x)

src/owl/signal/owl_signal.mli

@@ -0,0 +1,36 @@
+(*
+ * OWL - OCaml Scientific and Engineering Computing
+ * Copyright (c) 2021 Chandra Shekhar <[email protected]>, Kumar Appaiah <[email protected]>
+ *)
+
+(** Signal: Fundamental Signal Processing functions. *)
+
+open Owl_dense
+
+(** {Basic window functions} *)
+
+val blackman : int -> Ndarray.D.arr
+(** Blackman window is a taper formed by using the first three terms of a summation of cosines. It was designed to have close to the minimal leakage possible. 
+``blackman m`` returns a blackman window.
+*)
+
+val hamming : int -> Ndarray.D.arr
+(** Hamming window is a taper formed by using a raised cosine with non-zero endpoints, optimized to minimize the nearest side lobe. ``hamming m``
+returns a hamming window.
+*)
+
+
+val hann : int -> Ndarray.D.arr
+(** Hann window is a taper formed by using a raised cosine or sine-squared with ends that touch zero. ``hann m``
+returns a hann window.
+*)
+
+(** {Filter response function} *)
+
+val freqz : ?n:int -> ?whole:bool -> float array -> float array -> Ndarray.D.arr * Ndarray.Z.arr
+(** freqz computes the frequency response of a digital filter. 
+
+``freqz b a`` computes the frequency response of digital filter with numerator filter coeffecient given by ``b`` (float array) while the denominator filter coeffecient given by ``a`` (float array), and returns the frequencies and the frequency response respectively in real and complex ndarrays. Two optional parameters may be specified: ``n`` is an integer that determines the number of frequencies where the frequency response is to be evaluated, and ``whole`` is a boolean that decides whether the frequency response is two-sided or one-sided. Default values of ``n`` and ``whole`` are 512 and false.
+*)
+
+(*ends here*)

test/test_runner.ml

@@ -39,4 +39,5 @@ let () =
     ; "base: complex", Unit_base_complex.test_set
     ; "base: ndarray core", Unit_base_ndarray_core.test_set
     ; "base: linalg", Unit_linalg_solver.test_set
+    ; "base: signal", Unit_signal.test_set
     ; ("algodiff matrix", Unit_algodiff_matrix.[ Reverse.test; Forward.test ]) ]

test/unit_signal.ml

@@ -0,0 +1,81 @@
+(*Used unit_dense_ndarray test code as reference*)
+open Bigarray
+open Owl
+open Owl_signal
+
+module M = Owl_dense_ndarray_generic
+
+
+let ndarray = Alcotest.testable (fun _p (_x : (float, float64_elt) M.t) -> ()) M.equal
+
+let blackman_reference = M.zeros Float64 [|4|]
+let hamming_reference = M.zeros Float64 [|4|]
+let hann_reference = M.zeros Float64 [|4|]       
+
+let _ = 
+  M.set blackman_reference [|0|] 0.;
+  M.set blackman_reference [|1|] 0.63;
+  M.set blackman_reference [|2|] 0.63;
+  M.set blackman_reference [|3|] 0.;
+  M.set hamming_reference [|0|] 0.08;
+  M.set hamming_reference [|1|] 0.77;
+  M.set hamming_reference [|2|] 0.77;
+  M.set hamming_reference [|3|] 0.08;	
+  M.set hann_reference [|0|] 0.;
+  M.set hann_reference [|1|] 0.75;
+  M.set hann_reference [|2|] 0.75;
+  M.set hann_reference [|3|] 0.;
+
+(*a module with functions to test*)
+module To_test = struct
+  let blackman () = 
+    let b = blackman 4 in
+    let max_err = (Arr.map2 (fun x y -> x -. y) blackman_reference b |> Arr.abs |> Arr.max |> Arr.get) [|0|] in
+    max_err < 1e-5
+
+  let hamming () = 
+    let h = hamming 4 in
+    let max_err = (Arr.map2 (fun x y -> x -. y) hamming_reference h |> Arr.abs |> Arr.max |> Arr.get) [|0|] in
+    max_err < 1e-5
+
+  let hann () = 
+    let h = hann 4 in
+    let max_err = (Arr.map2 (fun x y -> x -. y) hann_reference h |> Arr.abs |> Arr.max |> Arr.get) [|0|] in
+    max_err < 1e-5
+
+  let freqz () =
+    let freqz_ref = M.zeros Complex64 [|4|] in
+    let freqz_num = M.zeros Float64 [|4|] in
+    let freqz_den = M.zeros Float64 [|4|] in
+    M.set freqz_num [|0|] (1. /. 6.);
+    M.set freqz_num [|1|] 0.5;
+    M.set freqz_num [|2|] 0.5;
+    M.set freqz_num [|3|] (1. /. 6.);
+    M.set freqz_den [|0|] 1.;
+    M.set freqz_den [|1|] 0.;
+    M.set freqz_den [|2|] (1. /. 3.);
+    M.set freqz_den [|3|] 0.;
+    M.set freqz_ref [|0|] {Complex.re = 1.0; im = 0.0}; (*ref values are taken from gnu Octave*)
+    M.set freqz_ref [|1|] {Complex.re = 0.6536; im = -0.7536};
+    M.set freqz_ref [|2|] {Complex.re = -0.5; im = -0.5};
+    M.set freqz_ref [|3|] {Complex.re = -0.0536; im = 0.0464};
+    let f = (freqz ~n:4 ~whole:false (freqz_num |> M.to_array) (freqz_den |> M.to_array))
+          |> (fun (a, b) -> b) in
+    let max_err = (Owl_dense_ndarray.Z.map2 (fun x a-> Complex.sub x a) f freqz_ref |> Owl_dense_ndarray.Z.abs |> Owl_dense_ndarray.Z.max |> Owl_dense_ndarray.Z.re |> Arr.get) [|0|] in
+    max_err < 1e-3
+
+end
+
+let blackman () = Alcotest.(check bool) "blackman" true (To_test.blackman ())
+
+let hamming () = Alcotest.(check bool) "hamming" true (To_test.hamming ())
+
+let hann () = Alcotest.(check bool) "hann" true (To_test.hann ())
+
+let freqz () = Alcotest.(check bool) "freqz" true (To_test.freqz ())
+
+let test_set =
+  [ "blackman", `Slow, blackman;
+    "hamming", `Slow, hamming;
+    "hann", `Slow, hann;
+    "freqz", `Slow, freqz; ]