Add freqz and freqz_test

src/owl/signal/owl_signal.ml

@@ -35,3 +35,39 @@ let blackman m =
   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 freq n =         (*n is the number of frequencies where freqz is to be calculated*)
+  let w = Ndarray.D.linspace 0. Owl_const.pi (n+1) in
+  Ndarray.D.get_slice [[0;n-1]] w
+
+let freqf n =         (*n is the number of frequencies where freqz is to be calculated (if whole is true)*)
+  let w = Ndarray.D.linspace 0. (2. *. Owl_const.pi) (2*n+1) in
+  Ndarray.D.get_slice [[0;(2*n-1)]] w
+
+
+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.Z.div y x
+  else	
+    let x = dtft a n in
+    let y = dtft b n in
+    Ndarray.Z.div y x

test/unit_signal.ml

@@ -42,6 +42,27 @@ module To_test = struct
     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) 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 ())
@@ -50,5 +71,7 @@ let hamming () = Alcotest.(check bool) "hamming" true (To_test.hamming ())
 
 let hann () = Alcotest.(check bool) "hann" true (To_test.hann ())
 
-let test_set = 
+let freqz () = Alcotest.(check bool) "freqz" true (To_test.freqz ())
+
+let test_set =
   [ "blackman", `Slow, blackman; "hamming", `Slow, hamming; "hann", `Slow, hann ]