Multidimensional FFT

Every digital image is a grid of numbers. When you blur a photo, sharpen an edge, or compress a JPEG, something remarkable happens under the hood: the image is decomposed into frequencies. Low frequencies carry the broad shapes; high frequencies carry sharp edges and fine texture. That decomposition is the Discrete Fourier Transform (DFT).

For a one-dimensional signal of length $N$, the DFT of a sequence $x$ is

$$X[k] = \sum_{n=0}^{N-1} x[n]\, e^{-2\pi i \frac{kn}{N}}, \quad k = 0, 1, \dots, N-1.$$

Computing all $N$ outputs naively takes $O(N^2)$ multiplications. The Fast Fourier Transform (FFT), introduced by Cooley and Tukey in 1965, reduces this to $O(N \log N)$ by recursively splitting even and odd indices — one of the most celebrated algorithms of the 20th century.

But images are two-dimensional, and medical scans or climate models are three-dimensional or more. How does the FFT extend? The key insight is separability: the multidimensional DFT can be computed by applying the 1D FFT repeatedly along each axis in turn. No new mathematics is needed — just the same elegant divide-and-conquer, swept across every dimension.

Draw on the left grid — click or drag to toggle pixels on and off. The right panel shows the 2D FFT magnitude spectrum (log scale, center-shifted so DC sits in the middle). The transform is recomputed live, axis by axis, every time you change a pixel.

<!-- {{c_layout_comment}} -->
<div class="demo-wrap">
  <div class="panel">
    <div class="panel-label">{{label_input}}</div>
    <canvas id="inputCanvas" width="128" height="128" title="{{canvas_title}}"></canvas>
    <div class="controls">
      <button id="btnClear" type="button">{{btn_clear}}</button>
      <button id="btnHStripe" type="button">{{btn_hstripe}}</button>
      <button id="btnVStripe" type="button">{{btn_vstripe}}</button>
      <button id="btnChecker" type="button">{{btn_checker}}</button>
    </div>
  </div>
  <div class="panel">
    <div class="panel-label">{{label_spectrum}}</div>
    <canvas id="specCanvas" width="128" height="128" aria-label="{{spectrum_aria}}"></canvas>
    <div class="info" id="info">{{info_default}}</div>
  </div>
</div>

/* {{c_base_styles}} */
* { box-sizing: border-box; margin: 0; padding: 0; }
body { font-family: system-ui, sans-serif; color: #222; background: #f5f7fa; }
.demo-wrap {
  display: flex; gap: 16px; flex-wrap: wrap;
  justify-content: center; padding: 12px;
}
.panel {
  display: flex; flex-direction: column; align-items: center; gap: 8px;
}
.panel-label {
  font-size: .82rem; font-weight: 600; color: #5a7088; letter-spacing: .03em;
  text-transform: uppercase;
}
canvas {
  border: 1.5px solid #cdd9e3; border-radius: 6px;
  cursor: crosshair; image-rendering: pixelated;
  width: 160px; height: 160px;
}
/* {{c_spectrum_style}} */
#specCanvas { cursor: default; }
.controls {
  display: flex; gap: 6px; flex-wrap: wrap; justify-content: center;
}
button {
  font: 600 12px system-ui; padding: .35rem .7rem;
  border: 1px solid #1d3557; background: #1d3557;
  color: #fff; border-radius: 6px; cursor: pointer;
}
button:hover { background: #16294a; }
.info {
  font-size: .8rem; color: #5a7088; text-align: center;
  min-height: 1.2em; max-width: 200px; line-height: 1.4;
}

// Code not found

Try drawing a horizontal stripe: energy concentrates along the vertical axis of the spectrum (vertical frequencies). Draw a vertical stripe and energy shifts to the horizontal axis. A checkerboard pattern pushes energy to the four corners — the highest frequencies in both directions. This is separability made visible: the 2D spectrum is the outer product of two 1D spectra.

Naive 2D DFT. An $N \times N$ image has $N^2$ pixels and $N^2$ frequencies. Computing each output frequency directly requires summing over all $N^2$ input pixels: $O(N^4)$ operations in total. For a modest $1024 \times 1024$ image that is roughly $10^{12}$ multiplications — far too slow for interactive use.

The row-column (separable) algorithm. Because the 2D DFT factors as

$$X[k_1, k_2] = \sum_{n_1=0}^{N-1} \left( \sum_{n_2=0}^{N-1} x[n_1, n_2]\, e^{-2\pi i \frac{k_2 n_2}{N}} \right) e^{-2\pi i \frac{k_1 n_1}{N}},$$

you can compute the inner sum (a 1D DFT over each row) first, then the outer sum (a 1D DFT over each column of the result). Each of the $N$ rows costs $O(N \log N)$ with an FFT; so does each of the $N$ columns. Total: $O(N \cdot N \log N + N \cdot N \log N) = O(N^2 \log N)$.

In $d$ dimensions. An $N^d$ volume requires $d \cdot N^{d-1}$ one-dimensional FFTs of length $N$, giving $O(d \cdot N^d \log N)$ — or just $O(N^d \log N)$ when $d$ is fixed. The saving over the naive $O(N^{2d})$ approach grows exponentially with dimension.

Status: solved. The row-column decomposition was known before Cooley and Tukey (1965), but the FFT made it practically fast. FFTW (Frigo & Johnson, 1997) later showed that the optimal algorithm depends on the machine's cache hierarchy, not just the asymptotic count — auto-tuning at install time to find the best plan for each specific $N$.

See also compression for how the related Discrete Cosine Transform (DCT) dominates image codecs, and pattern matching for how 2D convolution — equivalent to pointwise multiplication in frequency space — powers template search.

The multidimensional FFT is one of the most widely deployed algorithms on the planet:

• Image compression (JPEG/HEIC): the Discrete Cosine Transform — a close relative of the DFT — is applied to $8 \times 8$ pixel blocks. Low-frequency coefficients are kept; high-frequency ones are quantized away. The eye barely notices.
• MRI reconstruction: a magnetic resonance scanner directly samples a 3D Fourier space (k-space). An inverse 3D FFT turns those raw measurements into a volumetric image.
• Radio astronomy: interferometer arrays like the Very Large Array (VLA) and the Event Horizon Telescope measure the 2D Fourier transform of the sky. Imaging = inverse FFT.
• Convolution in deep learning: large convolutional kernels are sometimes applied in frequency space — multiply the FFTs of the image and the filter, then inverse-FFT. When the kernel is big enough, $O(N^2 \log N)$ beats the $O(N^2 k^2)$ direct convolution.
• Signal processing and radar: pulse compression, synthetic aperture radar (SAR), and sonar all exploit multidimensional FFTs to separate targets by range and Doppler shift.
• Crystallography: X-ray diffraction patterns are the squared magnitudes of the 3D Fourier transform of electron density. Recovering the structure from diffraction data is (roughly) an inverse FFT — though the missing phase information makes it a genuinely hard problem.

The common thread: whenever a physical process mixes contributions linearly across space or time, the Fourier basis diagonalizes it — turning convolution into multiplication and turning hard global problems into cheap local ones in frequency space.

The multidimensional FFT is a triumph of mathematical structure exploited by algorithm design. The 2D DFT looks daunting — $N^4$ naive operations for an $N \times N$ image — but separability reveals that it is nothing more than $2N$ independent 1D problems, each solvable by the Cooley-Tukey FFT in $O(N \log N)$. The saving is not a constant factor; it grows with $N$ and with the number of dimensions.

The same trick carries to three, four, and any number of dimensions without modification. Everything from your phone's camera pipeline to a hospital's MRI machine to a radio telescope array quietly runs on this one observation: the multidimensional Fourier transform separates along axes.

Draw a stripe, watch the spectrum snap to a perpendicular line, and you have seen the separability theorem with your own eyes. That perpendicular line is the whole story — repeated once per axis, at $O(N \log N)$ per pass, for as many dimensions as you need.