41) Extra: Fast transforms

1. Discrete Signals

2. Approximation by Fourier Basis

3. Fast Fourier Transform

using LinearAlgebra
using Plots
default(linewidth=3, legendfontsize=12, xtickfontsize=12, ytickfontsize=12)
using Pkg; Pkg.add("LaTeXStrings")
using LaTeXStrings

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

1. Discrete signals

Consider a function \(f(x)\) defined on an integer grid \(x \in \mathbb Z\).

We define the Fourier modes

\[ \phi(x, \theta) = e^{i\theta x} \]

for continuous \(\theta\).

Question: Is it possible to distinguish \(e^{i\theta x}\) from \(e^{i(\theta+2\pi)x}\)?

phi(x, theta) = exp(1im * theta * x)
rphi(x, theta) = real(phi(x, theta)) # real part
x = -4:4
theta = 3.14
plot([x -> rphi(x, theta), x -> rphi(x, theta+2π), x -> rphi(x, theta-2π)], label = [L"Re(\phi(x,\theta))" L"Re(\phi(x,\theta+ 2 \pi))" L"Re(\phi(x,\theta - 2 \pi))"])
scatter!(x, rphi.(x, theta), color=:black, label = L"Re(\phi(x, \theta))")

We say \(\theta = \pm \pi\) is the Nyquist Frequency for the integer grid.

• It corresponds to “two points per wavelength”

2. Approximation by Fourier basis

Just like we can approximate functions using linear combinations of polynomials, we can approximate periodic functions using a linear combination of Fourier modes.

\[ f(x) \approx \sum_{k=1}^n \underbrace{\hat f(\theta_k)}_{\hat f_k} e^{i\theta_k x} . \]

This is reminiscent of linear algebra… We saw how linear combinations can be expressed as column-wise matrix-vector products:

(23)\[\begin{align} \Bigg[ f(x) \Bigg] &= \Bigg[ e^{i\theta_1 x} \Bigg| e^{i\theta_2 x} \Bigg| \dotsm \Bigg] \begin{bmatrix} \hat f_1 \\ \hat f_2 \\ \vdots \end{bmatrix} \\ &= \Bigg[ e^{i\theta_1 x} \Bigg] \hat f_1 + \Bigg[ e^{i\theta_2 x} \Bigg] \hat f_2 + \dotsb . \end{align}\]

Continuous \(\theta\): infinite domain

If we take \(\theta \in (-\pi, \pi]\) as a continuous quantity (instead of a discrete set of modes), the sum becomes and integral and we get equality (for “nice enough” \(f(x)\)),

\[ f(x) = \int_{-\pi}^\pi \hat f(\theta) e^{i\theta x} d\theta, \]

in which \(\hat f(\theta)\) is the Fourier transform (specifically, the discrete time transform) of \(f(x)\). This representation is valuable for analyzing convergence of multigrid methods, among other applications.

Computing \(\hat f(\theta)\)

If we select a finite number of points \(x\) and compute the square Vandermonde matrix

\[ \mathcal F = \Bigg[ e^{i\theta_1 x} \Bigg| e^{i\theta_2 x} \Bigg| \dotsm \Bigg] \]

then, knowing the vector \(f\), we could solve

\[ \mathcal F \hat f = f \]

for \(\hat f\). This would require \(O(n^3)\) where \(n\) is the number of points.

function vander_fourier(x, n=nothing)
    if isnothing(n)
        n = length(x)
    end
    theta = LinRange(-pi + 2pi/n, pi, n)
    F = exp.(1im * x * theta')
end

x = LinRange(-20, 20, 41)
F = vander_fourier(x)
plot(x, imag.(F[:, 19:22]))

\(\mathcal F\) as a matrix

x = LinRange(-2, 2, 5)
F = vander_fourier(x) / sqrt(5)
@show norm(F' * F - I)
@show norm(F * F' - I);

norm(F' * F - I) = 4.0067422874877247e-16
norm(F * F' - I) = 

4.594822717058422e-16

• Every \(\mathcal F\) (suitably normalized) is a unitary matrix

  • a unitary matrix is the complex-valued generalization of “orthogonal matrix”

  • \(\mathcal F^H \mathcal F = \mathcal F \mathcal F^H = I\)

• Typical notation is \(\mathcal F^*\) or \(\mathcal F^H\) representing “Hermitian transpose” or conjugate transpose

# Remember, the ' in Julia is the Hermitian adjoint
vander_fourier([-1, 0, 1])'

3×3 adjoint(::Matrix{ComplexF64}) with eltype ComplexF64:
  0.5-0.866025im     1.0+0.0im   0.5+0.866025im
  0.5+0.866025im     1.0-0.0im   0.5-0.866025im
 -1.0+1.22465e-16im  1.0-0.0im  -1.0-1.22465e-16im

What does this mean for cost?

Fitting a discrete signal in the Fourier basis requires solving

\[\mathcal F \hat y = y\]

3. Fast Fourier Transform

In this discrete context, the transform we need to evaluate is

\[ \hat f_k = \sum_\ell e^{-i\theta_k x_\ell} f_\ell \]

where \(f_\ell\) are samples \(f(x_\ell)\) at integers \(x_\ell = \ell\) and \(\theta_k\) are the frequencies \(2 \pi k/n\) (because the branch \(\theta \in (-\pi, \pi]\) is arbitrary).

(24)\[\begin{align} \hat f_k &= \sum_{\ell=0}^{n-1} e^{-2\pi i \frac{k \ell}{n}} f_\ell \\ &= \underbrace{\sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k (2\ell)}{n}} f_{2\ell}}_{\text{even}} + \underbrace{\sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k (2\ell+1)}{n}} f_{2\ell+1}}_{\text{odd}} \\ &= \underbrace{\sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k \ell}{n/2}} f_{2\ell}}_{\text{transform of even data}} + e^{-2\pi i \frac k n} \underbrace{\sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k \ell}{n/2}} f_{2\ell+1}}_{\text{transform of odd data}} \end{align}\]

Periodicity

When the original signal \(f\) is periodic, \(f_{\ell} = f_{(\ell + n) \bmod n}\), then

(25)\[\begin{align} \hat f_{k+n/2} &= \sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{(k+n/2) \ell}{n/2}} f_{2\ell} + e^{-2\pi i \frac{(k+n/2)}{n}} \sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{(k+n/2) \ell}{n/2}} f_{2\ell+1} \\ &= \sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k \ell}{n/2}} f_{2\ell} + e^{-2\pi i \frac{k}{n}} e^{-\pi i} \sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k\ell}{n/2}} f_{2\ell+1} \\ &= \sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k \ell}{n/2}} f_{2\ell} - e^{-2\pi i \frac{k}{n}} \sum_{\ell=0}^{n/2-1} e^{-2\pi i \frac{k \ell}{n/2}} f_{2\ell+1} \end{align}\]

where we have used that

\[ e^{-2\pi i \ell} = \Big( e^{-2\pi i} \Big)^\ell = 1^\ell = 1 .\]

We’ve reduced a Fourier transform of length \(n\) (at a cost of \(n^2\)) to two transforms of length \(n/2\) (at a cost of \(2 n^2/4 = n^2/2\)). Repeating this recursively yields a complexity of \(O(n\log n)\).

Visualization

using Pkg
Pkg.add("FastTransforms")
using FastTransforms

n = 64; m = 62
x = 0:n-1
f = exp.(2im * pi * m * x / n)
plot(x, real.(f), ylim=(-1.1, 1.1), marker=:circle)

   Resolving package versions...

   Installed FastTransforms_jll ────── v0.6.3+0

   Installed Bessels ───────────────── v0.2.8
   Installed IntelOpenMP_jll ───────── v2025.2.0+0
   Installed BandedMatrices ────────── v1.10.1
   Installed LazyArrays ────────────── v2.8.0
   Installed FastTransforms ────────── v0.17.0
   Installed DSP ───────────────────── v0.8.4
   Installed oneTBB_jll ────────────── v2022.0.0+1
   Installed MKL_jll ───────────────── v2025.2.0+0
   Installed RecurrenceRelationships ─ v0.2.0
   Installed ArrayLayouts ──────────── v1.12.0
   Installed GenericFFT ────────────── v0.1.6
   Installed ToeplitzMatrices ──────── v0.8.5

   Installed FFTW ──────────────────── v1.10.0

    Updating `~/.julia/environments/v1.10/Project.toml`
  [057dd010] + FastTransforms v0.17.0
    Updating `~/.julia/environments/v1.10/Manifest.toml`
  [4c555306] + ArrayLayouts v1.12.0
  [aae01518] + BandedMatrices v1.10.1
  [0e736298] + Bessels v0.2.8
  [717857b8] + DSP v0.8.4
  [7a1cc6ca] + FFTW v1.10.0
  [057dd010] + FastTransforms v0.17.0
  [a8297547] + GenericFFT v0.1.6
  [c8e1da08] + IterTools v1.10.0
  [5078a376] + LazyArrays v2.8.0
  [807425ed] + RecurrenceRelationships v0.2.0
  [c751599d] + ToeplitzMatrices v0.8.5
  [f5851436] + FFTW_jll v3.3.11+0
  [34b6f7d7] + FastTransforms_jll v0.6.3+0
  [1d5cc7b8] + IntelOpenMP_jll v2025.2.0+0
  [856f044c] + MKL_jll v2025.2.0+0
  [1317d2d5] + oneTBB_jll v2022.0.0+1
  [4af54fe1] + LazyArtifacts
  [781609d7] + GMP_jll v6.2.1+6
  [3a97d323] + MPFR_jll v4.2.0+1

    Building FastTransforms → `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/b41969ccec1379b33967c9b720a250d4687cfc2d/build.log`

Precompiling packages...

    262.5 ms  ✓ RecurrenceRelationships

    311.1 ms  ✓ IntelOpenMP_jll
    319.6 ms  ✓ oneTBB_jll
    328.3 ms  ✓ FastTransforms_jll

    227.9 ms  ✓ RecurrenceRelationships → RecurrenceRelationshipsLinearAlgebraExt

    250.9 ms  ✓ RecurrenceRelationships → RecurrenceRelationshipsFillArraysExt

   1191.8 ms  ✓ Bessels

   6782.7 ms  ✓ ArrayLayouts

    697.5 ms  ✓ ArrayLayouts → ArrayLayoutsSparseArraysExt

   7327.6 ms  ✓ MKL_jll

   1471.3 ms  ✓ LazyArrays

    780.6 ms  ✓ LazyArrays → LazyArraysStaticArraysExt

   2821.9 ms  ✓ FFTW

    969.2 ms  ✓ RecurrenceRelationships → RecurrenceRelationshipsLazyArraysExt

    419.2 ms  ✓ GenericFFT

    635.6 ms  ✓ Polynomials → PolynomialsFFTWExt

   1479.5 ms  ✓ DSP

   6647.9 ms  ✓ BandedMatrices

    900.6 ms  ✓ DSP → OffsetArraysExt

   1050.3 ms  ✓ ToeplitzMatrices

    830.7 ms  ✓ ArrayInterface → ArrayInterfaceBandedMatricesExt

    766.9 ms  ✓ ToeplitzMatrices → ToeplitzMatricesStatsBaseExt
    998.9 ms  ✓ BandedMatrices → BandedMatricesSparseArraysExt

   1233.5 ms  ✓ LazyArrays → LazyArraysBandedMatricesExt

   2584.5 ms  ✓ FastTransforms
  25 dependencies successfully precompiled in 19 seconds. 330 already precompiled.
  1 dependency had output during precompilation:
┌ MKL_jll
│   Downloading artifact: IntelOpenMP
│   Downloading artifact: oneTBB
└  

fhat = fft(f) # built-in Julia function int the FastTransforms package
scatter([abs.(fhat), abs.(fftshift(fhat))], marker=[:square :circle])

UndefVarError: `fft` not defined

Stacktrace:
 [1] top-level scope
   @ In[7]:1

Example: Transform a Gaussian bump \(e^{-x^2}\)

n = 64; w = 1
x = 0:n-1
f = exp.(-((x .- n/2)/w) .^ 2)
scatter(x, real.(f), ylim=(-1.1, 1.1))

fhat = fft(f)
scatter([abs.(fhat), abs.(fftshift(fhat))])

Compute derivatives using the Fourier transform

How do we differentiate this?

\[f(x) = \sum_k e^{i\theta_k x} \hat f_k\]

• Evidently, we need only compute

  \[f'(x) = \mathcal F^{-1} (i \theta_k \hat f_k)\]

Generalizations

• Non-power of 2

• Non-uniform grids

• Multiple dimensions

• Butterfly algorithms for integral operators

  \[ (\mathcal K f)(x) = \int_{Y} K(x,y) f(y) dy \]
  See Poulson et al., A parallel butterfly algorithm.

Applications

• Everywhere in signal processing

• ECMWF global climate and weather model

• Particle-Mesh Ewald for long-range forces in molecular dynamics

• Turbulence simulation

Parallel implications: Bisection bandwidth of the network

function vander_chebyshev(x, n=nothing)
    if isnothing(n)
        n = length(x) # Square by default
    end
    m = length(x)
    T = ones(m, n)
    if n > 1
        T[:, 2] = x
    end
    for k in 3:n
        #T[:, k] = x .* T[:, k-1]
        T[:, k] = 2 * x .* T[:,k-1] - T[:, k-2]
    end
    T
end

CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(0, π, n))

Derivative of a function on the Chebyshev (CosRange) points

function chebfft(y)
    # Adapted from Trefethen's Spectral Methods in Matlab
    n = length(y) - 1
    x = CosRange(-1, 1, n+1)
    Y = [y; reverse(y[2:n])]
    ii = 0:n-1
    U = real.(fft(Y))
    W = real.(ifft(1im * [ii; 0; 1-n:-1] .* U))
    w = zeros(n + 1)
    w[2:n] = -W[2:n] ./ sqrt.(1 .- x[2:n] .^ 2)
    w[1] = sum(ii .^ 2 .* U[ii .+ 1]) / n + .5 * n * U[n+1]
    w[n+1] = sum((-1) .^ (ii .+ 1) .* ii .^ 2 .* U[ii .+ 1]) / n + .5 * (-1) ^ (n+1) * n * U[n+1]
    w
end

chebfft([1,2,0,3])

x = CosRange(-1, 1, 40)
g(x) = cos(3x - 1)
dg(x) = -3 * sin(3x - 1) # for verification
y = g.(x)
dy = chebfft(y)
plot(x, dy, marker=:circle)
plot!(dg)