The Discrete Fourier Transform (DFT)

Mix.install([
  {:nx_signal, "~> 0.2"},
  {:vega_lite, "~> 0.1"},
  {:kino_vega_lite, "~> 0.1"}
])

What is the Discrete Fourier Transform (DFT)?

This livebook will show you how to use the Discrete Fourier Transform (DFT) to analyze a signal.

Suppose we have a periodic signal which we want to analyze.

We will run a Fast Fourrier Transform, which is a fast algorithm to compute the DFT. It transforms a time-domain function into the frequency domain.

Building the signal

Let’s build a known signal that we will decompose and analyze later on.

The signal will be the sum of two sinusoidal signals, one at 5Hz, and one at 20Hz with the corresponding amplitudes (1, 0.5).

$f(t) = \sin(2\pi5t) + \frac{1}{2} \sin(2\pi20t)$

Suppose we can sample at fs=50Hz (meaning 50 samples per second) and our aquisition time is duration = 1s.

We build a time series of t equally spaced points with the given duration interval with Nx.linspace.

For each value of this serie (the discrete time $t$), we will synthesize the signal $f(t)$ through the module below.

defmodule Signal do
  import Nx.Defn
  import Nx.Constants, only: [pi: 0]

  defn source(t) do
    f1 = 5
    f2 = 20
    Nx.sin(2 * pi() * f1 * t) + 1/2 * Nx.sin(2 * pi() * f2 * t)
  end

  def sample(opts) do
    fs = opts[:fs]
    duration = opts[:duration]
    t = Nx.linspace(0, duration, n: trunc(duration * fs), endpoint: false, type: {:f, 32})
    source(t)
  end
end

We sample our signal at fs=50Hz during 1s:

fs = 50; duration= 1

sample = Signal.sample(fs: fs, duration: 1)

Analyzing the signal with the DFT

Because our signal contains many periods of the underlying function, the DFT results will contain some noise. This noise can stem both from the fact that we’re likely cutting of the signal in the middle of a period and from the fact that we have a specific frequency resolution which ends up grouping our individual components into frequency bins. The latter isn’t really a problem as we have chosen fs to be fast enough.

> The number at the index $i$ of the DFT results gives an approximation of the amplitude and phase of the sampled signal at the frequency $i$.

In other words, doing Nx.fft(sample) returns a list of numbers indexed by the frequency.

dft = Nx.fft(sample)

We will limit our study points to the first half of dft because it is symmetrical on the upper half. The phase doesn’t really matter to us because we don’t wish to reconstruct the signal, nor find possible discontinuities, so we’ll use Nx.abs to obtain the absolute values at each point.

n = Nx.size(dft)

max_freq_index = div(n, 2)

amplitudes =  Nx.abs(dft)[0..max_freq_index]

# the frequency bins, "n" of them spaced by fs/n=1 unit:
frequencies = NxSignal.fft_frequencies(fs, fft_length: n)[0..max_freq_index]

data1 = %{
  frequencies: Nx.to_list(frequencies),
  amplitudes: Nx.to_list(amplitudes)
}

VegaLite.new(width: 700, height: 300)
|> VegaLite.data_from_values(data1)
|> VegaLite.mark(:bar, tooltip: true)
|> VegaLite.encode_field(:x, "frequencies",
  type: :quantitative,
  title: "frequency (Hz)",
  scale: [domain: [0, 50]]
)
|> VegaLite.encode_field(:y, "amplitudes",
  type: :quantitative,
  title: "amplitutde",
  scale: [domain: [0, 30]]
)

We see the peaks at 5Hz, 20Hz with their amplitudes (the second is half the first).

This is indeed our synthesized signal 🎉

We can confirm this visual inspection with a peek into our data. We use Nx.top_k function.

{values, indices} = Nx.top_k(amplitudes, k: 5)

{
  values, 
  frequencies[indices]
}

Visualizing the original signal and the Inverse Discrete Fourier Transform

Let’s visualize our incoming signal over 400ms. This correspond to 2 periods of our 5Hz component and 8 periods of our 20Hz component.

We compute 200 points to have a smooth curve, thus every (400/200=) 2ms.

We also add the reconstructed signal via the Inverse Discrete Fourier Transform available as Nx.ifft.

This gives us 50 values spaced by 1000ms / 50 = 20ms.

Below, we display them as a bar chart under the line representing the ideal signal.

#----------- REAL SIGNAL
# compute 200 points of the "real" signal during 2/5=400ms (twice the main period)

t = Nx.linspace(0, 0.4, n: trunc(0.4 * 500))
sample = Signal.source(t)

#----------- RECONSTRUCTED IFFT
yr = Nx.ifft(dft) |> Nx.real()
fs = 50
tr = Nx.linspace(0, 1, n: 1 * fs, endpoint: false)

idx = Nx.less_equal(tr, 0.4)
xr = Nx.select(idx, tr, :nan)
yr = Nx.select(idx, yr, :nan)
#----------------


data = %{
  x: Nx.to_list(t),
  y: Nx.to_list(sample)
}

data_r = %{
  yr: Nx.to_list(yr),
  xr: Nx.to_list(xr)
}

VegaLite.new(width: 600, height: 300)
|> VegaLite.layers([
  VegaLite.new()
  |> VegaLite.data_from_values(data)
  |> VegaLite.mark(:line, tooltip: true)
  |> VegaLite.encode_field(:x, "x", type: :quantitative, title: "time (ms)", scale: [domain: [0, 0.4]])
  |> VegaLite.encode_field(:y, "y", type: :quantitative, title: "signal")
  |> VegaLite.encode_field(:order, "x"),
  VegaLite.new()
  |> VegaLite.data_from_values(data_r)
  |> VegaLite.mark(:bar, tooltip: true)
  |> VegaLite.encode_field(:x, "xr", type: :quantitative, scale: [domain: [0, 0.4]])
  |> VegaLite.encode_field(:y, "yr", type: :quantitative, title: "reconstructed")
  |> VegaLite.encode_field(:order, "xr")
])

We see that during 400ms, we have 2 periods of a longer period signal, and 8 of a shorter and smaller perturbation period signal.