7 Fourier Transform
Prerequisites: the discrete Fourier transform and complex numbers. The chapter uses standard concepts (spectrum, frequency bins, convolution theorem) without introduction.
La Linea bridges Fastmath’s transform API to ComplexTensors. Fastmath’s (:complex :fftr) transformer outputs interleaved double[] with [re₀ im₀ re₁ im₁ ...] — exactly the memory layout of ComplexTensor. So the bridge is zero-copy.
(ns lalinea-book.fourier-transform
(:require
;; La Linea (https://github.com/scicloj/lalinea):
[scicloj.lalinea.tensor :as t]
[scicloj.lalinea.elementwise :as el]
;; FFT bridge — Fastmath transforms <-> ComplexTensor:
[scicloj.lalinea.transform :as ft]
;; Dataset manipulation (https://scicloj.github.io/tablecloth/):
[tablecloth.api :as tc]
;; Interactive Plotly charts (https://scicloj.github.io/tableplot/):
[scicloj.tableplot.v1.plotly :as plotly]
;; Visualization annotations (https://scicloj.github.io/kindly-noted/):
[scicloj.kindly.v4.kind :as kind]
[clojure.math :as math]))Forward FFT
A real signal produces a complex spectrum.
(ft/forward [1.0 0.0 -1.0 0.0])#la/C [:float64 [4] [0.000 + 0.000 i 2.000 + 0.000 i 0.000 + 0.000 i 2.000 + 0.000 i]]The DC component (\(k=0\)) is \(\sum x_n = 0\).
Round-trip
FFT followed by inverse FFT recovers the original signal.
(let [signal [1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0]
spectrum (ft/forward signal)
recovered (ft/inverse-real spectrum)]
(el/reduce-max (el/abs (el/- recovered signal))))0.0Parseval’s theorem
Energy in time domain equals energy in frequency domain (up to a scale factor of \(N\)):
\[\sum_n |x_n|^2 = \frac{1}{N} \sum_k |\hat{x}_k|^2\]
(let [signal [1.0 2.0 3.0 4.0]
n (count signal)
spectrum (ft/forward signal)
time-energy (el/sum (el/* signal signal))
magnitudes (el/abs spectrum)
freq-energy (/ (el/sum (el/* magnitudes magnitudes)) n)]
(< (abs (- time-energy freq-energy)) 1e-10))trueLinearity
\(\mathcal{F}(\alpha x + \beta y) = \alpha \mathcal{F}(x) + \beta \mathcal{F}(y)\)
(let [x [1.0 2.0 3.0 4.0]
y [5.0 6.0 7.0 8.0]
alpha 2.0
beta -1.5
combined (el/+ (el/* alpha x) (el/* beta y))
lhs (ft/forward combined)
rhs (el/+ (el/scale (ft/forward x) alpha)
(el/scale (ft/forward y) beta))]
(and (< (el/reduce-max (el/abs (el/- (el/re lhs) (el/re rhs)))) 1e-10)
(< (el/reduce-max (el/abs (el/- (el/im lhs) (el/im rhs)))) 1e-10)))trueConvolution theorem
Pointwise multiplication in frequency domain corresponds to circular convolution in time domain:
\[\mathcal{F}(x * y) = \mathcal{F}(x) \cdot \mathcal{F}(y)\]
(let [x [1.0 2.0 0.0 0.0]
y [1.0 0.0 1.0 0.0]
Fx (ft/forward x)
Fy (ft/forward y)
product-spectrum (el/* Fx Fy)
conv-result (ft/inverse-real product-spectrum)
n (count x)
manual-conv (let [out (t/make-container :float64 n)]
(dotimes [k n]
(let [s (loop [j 0 acc 0.0]
(if (>= j n) acc
(recur (inc j)
(+ acc (* (double (x j))
(double (y (mod (- k j) n))))))))]
(t/set-value! out k s)))
out)]
(< (el/reduce-max (el/abs (el/- conv-result manual-conv))) 1e-10))trueVisualising a spectrum
A signal made of two sine waves — 3 Hz and 7 Hz:
(def N-vis 64)(def signal-composed
(t/->real-tensor
(t/materialize
(t/make-reader :float64 N-vis
(let [ti (/ (double idx) N-vis)]
(+ (math/sin (* 2 math/PI 3 ti))
(* 0.5 (math/sin (* 2 math/PI 7 ti)))))))))The time-domain waveform:
(-> (tc/dataset {:t (t/make-reader :float64 N-vis (/ (double idx) N-vis))
:amplitude signal-composed})
(plotly/base {:=x :t :=y :amplitude})
(plotly/layer-line)
plotly/plot)The magnitude spectrum reveals two peaks at frequencies 3 and 7:
(let [spectrum (ft/forward signal-composed)
mags (el/abs spectrum)]
(-> (tc/dataset {:frequency (range (/ N-vis 2))
:magnitude (take (/ N-vis 2) mags)})
(plotly/base {:=x :frequency :=y :magnitude})
(plotly/layer-bar)
plotly/plot))The two largest magnitude bins are at frequencies 3 and 7:
(let [spectrum (ft/forward signal-composed)
mags (el/abs spectrum)
half-n (/ N-vis 2)
peak-idx (el/argsort > (t/select mags (range half-n)))]
(= [3 7] (sort (take 2 peak-idx))))trueKnown transform pairs
The DFT of a constant signal \(x_n = c\) is \(\hat{x}_0 = Nc\) with all other bins zero.
(let [spectrum (ft/forward [3.0 3.0 3.0 3.0])]
{:dc (el/re (spectrum 0))
:others [(el/abs (spectrum 1)) (el/abs (spectrum 2)) (el/abs (spectrum 3))]}){:dc 12.0, :others [0.0 0.0 0.0]}The DFT of \(x = [1, -1, 1, -1]\) has energy only at Nyquist (\(k = N/2\)).
(let [spectrum (ft/forward [1.0 -1.0 1.0 -1.0])]
{:dc (double (el/abs (spectrum 0)))
:nyquist (double (el/re (spectrum 2)))}){:dc 0.0, :nyquist 4.0}Complex-to-complex FFT
When the input is already complex, use forward-complex.
(let [signal (t/complex-tensor [1.0 0.0] [0.0 1.0])
spectrum (ft/forward-complex signal)
recovered (ft/inverse spectrum)]
(and (< (el/reduce-max (el/abs (el/- (el/re recovered) (el/re signal)))) 1e-10)
(< (el/reduce-max (el/abs (el/- (el/im recovered) (el/im signal)))) 1e-10)))trueOther transforms
DCT, DST, and DHT are also available, returning real tensors.
(let [signal [1.0 2.0 3.0 4.0]
dct (ft/dct-forward signal)
recovered (ft/dct-inverse dct)]
(< (el/reduce-max (el/abs (el/- recovered signal))) 1e-10))true