FFT module revamp + updating dune-lang (#680)

* FFT module revamp.

- Added several new functionnalities to the FFT module by changing the dependency from FFTPACK to POCKETFTT.
	- new optionnal parameters for the API (nthreads, norm, ...)
	- new functions for cosine and sine transforms (dct, dst, ...)
- Switched from dune 2.0 to dune 3.16 (this was required as I ran throught issues with linking while using 2.0)

* Adding abstract types and documentation.

- Added the ttrig_transform type to specify the DCT and DST types
- Added the tnorm type to specify the normalization option for the FFTs.

* Adding 2 FFT example usage.

- One "complex" usage: computing a PSD spectrogram from input data
- One "simple" usage: computing the maximum frequency peak in time series data using rfft.

* Adding tests, fixing minor issues.

- Fixed an issue where the DCT/DST parameters weren't passed in the right order
- Fixed an issue where the norm factor of the DCT wasn't computed correctly (incorrect delta)
- Generated a unit_fft file that tests various parameters of FFT. Values are generated using scipy.fft module for the FFT, DCT and DST functions
- Changed pocketfft::detail:: namespace usage to just pocketfft:: in the C++ code for more readability

* Remove pocketfft submodule from .gitmodules

* Convert pocketfft submodule to regular directory

* Ubuntu 4.x build fix. MacOS fix attempt.

* Moved the definition outside extern C.
- Kept the declaration inside the extern C but moved away the declaration, as an attempt to fix MacOS builds.

* owl_core preprocessor for cpp compiling.

* Adapted the fft2 and ifft2 prototypes to match new FFT module.

- Added norm and nthreads parameters.

* Fixing documentation strings.

* Little code change for more similarities with scipy's implementation.

Author: gabyfle

.gitmodules

@@ -1,3 +1,3 @@
-(lang dune 2.0)
+(lang dune 3.16)
 
 (name owl)

dune-project

@@ -1,5 +1,5 @@
 (executables
- (names 
+ (names
   backprop
   checkpoint
   cifar10_vgg
@@ -19,12 +19,23 @@
   lazy_mnist
   linear_algebra
   lstm
+  max_freq
   mnist_cnn
   mnist_lenet
   newton_method
+  specgram
   squeezenet
   test_log
   tfidf
-  vgg16
- )
- (libraries owl))
+  vgg16)
+ (libraries owl)
+ (flags ; in order to make the examples compile correctly even with the warnings.
+  (:standard
+   -warn-error
+   -unused-value-declaration
+   -warn-error
+   -unused-var-strict
+   -warn-error
+   -unused-var
+   -warn-error
+   -unused-field)))

examples/dune

@@ -0,0 +1,21 @@
+(** This example shows how to compute the maximum peak frequency in time series data using the FFT module. *)
+
+module G = Owl.Dense.Ndarray.Generic
+module FFT = Owl.Fft.Generic
+
+let max_freq signal sampling_rate =
+  (* Apply FFT *)
+  let fft_result = FFT.rfft ~otyp:Bigarray.Complex32 signal in
+  (* Get magnitude spectrum *)
+  let magnitudes = G.abs fft_result in
+  (* Find peak frequency *)
+  let max_idx = ref 0 in
+  let max_val = ref (G.get magnitudes [|0|]) in
+  for i = 0 to G.numel magnitudes - 1 do
+    let curr_val = G.get magnitudes [|i|] in
+    if curr_val > !max_val then (
+      max_val := curr_val ;
+      max_idx := i )
+  done ;
+  (* Convert index to frequency *)
+  float_of_int !max_idx *. sampling_rate /. float_of_int (G.numel signal)

examples/max_freq.ml

@@ -0,0 +1,62 @@
+(** This example, extracted from the SoundML library, shows how to compute a specgram using the FFT module *)
+
+module G = Owl.Dense.Ndarray.Generic
+
+(* helper to compute the fft frequencies *)
+let fftfreq (n : int) (d : float) =
+  let nslice = ((n - 1) / 2) + 1 in
+  let fhalf =
+    G.linspace Bigarray.Float32 0. (float_of_int nslice) nslice
+  in
+  let shalf =
+    G.linspace Bigarray.Float32 (-.float_of_int nslice) (-1.) nslice
+  in
+  let v = G.concatenate ~axis:0 [|fhalf; shalf|] in
+  Owl.Arr.(1. /. (d *. float_of_int n) $* v)
+
+(* Computes a one-sided PSD spectrogram with no padding and no detrend *)
+let specgram (nfft : int) (fs : int) ?(noverlap : int = 0) (x : (float, Bigarray.float32_elt) G.t) =
+  let window = Owl.Signal.hann in (* we're using hann window *)
+  assert (noverlap < nfft) ;
+  (* we're making copies of the data from x and y to then use in place padding
+     and operations *)
+  let x = G.copy x in
+  (* We're making sure the arrays are at least of size nfft *)
+  let xshp = G.shape x in
+  ( if Array.get xshp 0 < nfft then
+      let delta = nfft - Array.get xshp 0 in
+      (* we're doing this in place in hope to gain a little bit of speed *)
+      G.pad_ ~out:x ~v:0. [[0; delta - 1]; [0; 0]] x ) ;
+  let scale_by_freq = true  in
+  let pad_to = nfft in
+  let scaling_factor = 2. in
+  let window = window nfft |> G.cast_d2s in
+  let window =
+    G.reshape G.(window * ones Bigarray.float32 [|nfft|]) [|-1; 1|]
+  in
+  let res =
+    G.slide ~window:nfft ~step:(nfft - noverlap) x |> G.transpose
+  in
+  (* if we'd add a detrend, we'd need to do it before applying the window *)
+  let res = G.(res * window) in
+  (* here comes the rfft compute, if you're processing large audio data, you might want to set ~nthreads to
+     something that suits both your hardware and your needs. *)
+  let res = Owl.Fft.S.rfft res ~axis:0 in
+  let freqs = fftfreq pad_to (1. /. float_of_int fs) in
+  let conj = G.conj res in
+  (* using in-place operations to avoid array copy *)
+  G.mul_ ~out:res conj res;
+  let slice = if nfft mod 2 = 0 then [[1; -1]; []] else [[1]; []] in
+  let gslice = G.get_slice slice res in
+  G.mul_scalar_ ~out:gslice gslice Complex.{re= scaling_factor; im= 0.} ;
+  G.set_slice slice res gslice ;
+  if scale_by_freq then (
+    let window = G.abs window in
+    G.div_scalar_ ~out:res res Complex.{re= float_of_int fs; im= 0.} ;
+    let n = G.sum' (G.pow_scalar window (float_of_int 2)) in
+    G.div_scalar_ ~out:res res Complex.{re= n; im= 0.} )
+  else (
+    let window = G.abs window in
+    let n = Float.pow (G.sum' window) 2. in
+    G.div_scalar_ ~out:res res Complex.{re= n; im= 0.} ) ;
+  (res, freqs)

examples/specgram.ml

@@ -13,7 +13,7 @@ type 'a node =
     mutable next : 'a node array
   ; (* children of the node *)
     mutable attr : 'a (* indicate the validity *)
-  }
+  } [@@warning "-69"]
 
 type order =
   | BFS