Source file `owl_signal.ml`

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
# 1 "src/owl/signal/owl_signal.ml"
(*
 * OWL - OCaml Scientific Computing
 * Copyright (c) 2016-2022 Liang Wang <liang@ocaml.xyz>
 *)

(** 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)