Discrete Fourier transform (DFT) is a widely-used approach of frequency analysis. DFT is defined by
$$z_k = \frac{1}{N} \sum^{N-1}_{t=0} x_t \exp\left( - \mathrm{i} \frac{2 \pi k t}{N} \right)$$where $N$ is the number of samples, $x_t \in \mathbb{C}$ ($t=0,\dots,N-1$) is a series of time-domain samples, and $z_k \in \mathbb{C}$ ($k=0,\dots,N-1$) is a series of frequency-domain points. The index $k$ is corresponds to frequency $f_k = f_s k / N$ where $f_s$ is a sample rate, and absolute value $|z_k|$ is amplitude of $f_k$.
Fast Fourier transform (FFT) is a set of algorithms for efficient computation of DFT. FFTW3 is a major FFT library faster than similar libraries. We show a simple example of fftw3-ocaml, a binding to FFTW3.
In [1]:
#require "fftw3" ;;
#use "archimedes_iocaml.ml" ;;
open Bigarray ;;
module FFT = Fftw3.D ;;
Out[1]:
In [2]:
let pi = 3.14159265358979 ;;
let n = 1000 ;; (* # of samples *)
let xdata = Array1.create FFT.complex C_layout n ;; (* source memory for time-domain data *)
let ydata = Array1.create FFT.complex C_layout n ;; (* destination memory for frequency-domain data *)
let plan = FFT.Array1.dft FFT.Forward xdata ydata ;; (* FFTW3 plan *)
Out[2]:
Out[2]:
Out[2]:
Out[2]:
Out[2]:
Out[2]:
This example shows the FFT of the hamming window given by
$$w(x) = 0.54 - 0.46 \cos(2 \pi x).$$We append zeros at the beginning and the end of the signal for simulation of masking by the window function.
In [3]:
let min = n / 10 in
let max = n - min in
for i = 0 to n - 1 do
if min <= i && i <= max then begin
let x = float (i - min) /. float (max - min) in
let y = 0.54 -. 0.46 *. cos (2.0 *. pi *. x) in (* the hamming window between `min` and `max` *)
xdata.{i} <- Complex.({ re = y; im = 0.0; })
end else begin
xdata.{i} <- Complex.zero
end
done
Out[3]:
In [4]:
let vp = A.init ~w:600. ~h:300. ["iocaml"] in
A.Axes.box vp ;
A.Viewport.title vp "Time domain" ;
A.Viewport.ylabel vp "Amplitude" ;
A.Viewport.yrange vp 0.0 1.0 ;
A.set_color vp A.Color.red ;
A.Array.y vp ~style:`Lines (Array.init n (fun i -> xdata.{i}.Complex.re)) ;
A.close vp
Out[4]:
In [5]:
FFT.exec plan ;; (* Execute FFT *)
Out[5]:
In [6]:
let spectrum = Array.init (n / 2) (fun i -> 10. *. log10 (Complex.norm2 ydata.{i})) in
let vp = A.init ~w:600. ~h:300. ["iocaml"] in
A.Axes.box vp ;
A.Viewport.title vp "Frequency domain" ;
A.Viewport.ylabel vp "Magnitude (dB)" ;
A.set_color vp A.Color.red ;
A.Array.y vp ~style:`Lines spectrum ;
A.close vp
Out[6]: