The Fourier View of Images

The Fourier transform re-expresses a signal as a combination of complex exponentials (equivalently, sines and cosines) instead of as a function of position. For a 1D continuous signal `f(x)`, the transform is ``` F(u) = ∫ f(x) e^{-i2πux} dx ``` Here `u` is a spatial frequency — cycles per unit length — and `F(u)` is a complex number encoding, for that frequency, an amplitude (how strong the oscillation is) and a phase (where the oscillation is offset). For a 2D discrete image `f[m,n]` of size `M × N`, the analogous Discrete Fourier Transform (DFT) is ``` F[u,v] = Σ_{m=0}^{M-1} Σ_{n=0}^{N-1} f[m,n] · exp(−i2π(um/M + vn/N)) ``` and the image is recoverable exactly via the inverse DFT ``` f[m,n] = (1/MN) Σ_{u=0}^{M-1} Σ_{v=0}^{N-1} F[u,v] · exp(+i2π(um/M + vn/N)) ``` Each basis function `exp(−i2π(um/M + vn/N))` is a 2D sinusoidal grating: `u` controls how fast the wave oscillates vertically, `v` controls horizontal oscillation, and the combination gives a wave of a particular orientation and spatial frequency. Low `(u,v)` near the origin correspond to slowly varying image content (broad shading, smooth gradients); high `(u,v)` far from the origin correspond to fine detail (edges, texture, noise). Every `F[u,v]` is a complex number, usually written in polar form as `|F[u,v]| e^{iφ[u,v]}`. The magnitude spectrum `|F[u,v]|` measures how much energy the image has at that frequency and orientation, independent of where in the image it occurs. The phase spectrum `φ[u,v]` encodes the spatial alignment of that frequency component — where its peaks and troughs fall. A classic and important demonstration is that swapping the phase of two unrelated images (keeping each one's own magnitude) produces an image that looks structurally like the one whose phase was used, not the one whose magnitude was used. Phase carries edges and shapes; magnitude carries texture-like statistics. This matters directly for Linear Filtering and Convolution: a filter that manipulates only magnitude (like an ideal low-pass) mostly changes texture and sharpness, while any perturbation to phase can shift or distort edges.

The convolution theorem

Let `f` be an image and `h` a filter kernel, and let `` denote discrete 2D convolution (as introduced in <a class="wiki-link" href="/wiki/cs-cv-linear-filtering-and-convolution" data-track="wiki_article.body.wiki_link">Linear Filtering and Convolution</a>): ``` (f h)[m,n] = Σ_j Σ_k f[j,k] · h[m−j, n−k] ``` The convolution theorem states: ``` DFT{ f * h } = DFT{f} · DFT{h} (pointwise complex multiplication) ``` and symmetrically, multiplication in space corresponds to convolution in frequency. This gives two equivalent ways to apply a filter: 1. Spatial domain: slide `h` over `f`, computing a weighted sum at every pixel — cost `O(N² k²)` for an `N×N` image and `k×k` kernel. 2. Frequency domain: compute `F = DFT(f)`, `H = DFT(h)`, multiply pointwise, then invert — cost dominated by the FFT, `O(N² log N)`. For small kernels (3×3, 5×5) the spatial approach wins. For large kernels — a big Gaussian blur, a motion-blur kernel spanning hundreds of pixels, or correlating against a large template — the frequency approach is dramatically faster, and this crossover is a standard engineering decision in image-processing libraries.

Why the theorem is true (proof sketch)

Substitute the convolution definition into the DFT and swap the order of summation: ``` DFT{f*h}[u,v] = Σ_m Σ_n ( Σ_j Σ_k f[j,k] h[m−j,n−k] ) e^{−i2π(um/M+vn/N)} ``` Re-index with `m' = m−j`, `n' = n−k` (a valid change of variables on a periodic/circular domain, which is exactly what the DFT assumes): ``` = Σ_j Σ_k f[j,k] e^{−i2π(uj/M+vk/N)} · Σ_{m'} Σ_{n'} h[m',n'] e^{−i2π(um'/M+vn'/N)} = F[u,v] · H[u,v] ``` The key algebraic fact enabling the split is that the complex exponential is multiplicative: `e^{−i2πu(j+m')/M} = e^{−i2πuj/M} · e^{−i2πum'/M}`. This is why sinusoids — and only sinusoids — are the eigenfunctions of convolution: convolving a pure sinusoid with any filter yields the same sinusoid, only rescaled and phase-shifted. That eigenfunction property is the entire reason the theorem holds, and it is why the frequency domain is the natural coordinate system for any linear shift-invariant (LSI) system, exactly the class of systems studied in Linear Filtering and Convolution.

Reading a filter's frequency response

Because convolving with `h` means multiplying by `H = DFT(h)`, the function `H[u,v]` is literally the filter's transfer function: at each frequency, it tells you the gain (how much that frequency is amplified or attenuated) and phase shift applied. This lets you classify filters immediately:

• Low-pass (`H` near 1 at low frequencies, near 0 at high frequencies): blurs, removes noise and fine texture. Gaussian and box blurs are low-pass.
• High-pass (`H` near 0 at low frequencies, near 1 at high): removes flat shading, keeps edges and noise. The Laplacian is high-pass.
• Band-pass: keeps a middle range of frequencies — used for texture and edge detectors tuned to particular scales, e.g., Gabor filters, difference-of-Gaussians.

Algorithm: frequency-domain filtering

``` function FREQUENCY_FILTER(image f, kernel h):

# 1. Pad kernel to image size (zero-pad), align origin H_full = zero_pad(h, size(f)) H_full = circular_shift(H_full, -center(h)) # so DFT phase is correct

# 2. Transform both to frequency domain F = FFT2D(f) H = FFT2D(H_full)

# 3. Pointwise complex multiplication G = F .* H # elementwise, complex

# 4. Inverse transform back to spatial domain g = Real( IFFT2D(G) ) return g ```

Correctness caveat — linear vs. circular convolution. The DFT is inherently periodic: it implicitly treats the image as if it tiled infinitely, so multiplication in frequency corresponds to circular convolution, which wraps content from one edge of the image onto the other. To reproduce ordinary ("linear") convolution — the kind used in spatial filtering, where the kernel does not wrap around — the image and kernel must be zero-padded to at least `M + k − 1` by `N + k − 1` before transforming. Skipping this step is one of the most common correctness bugs in frequency-domain implementations.

Algorithm: the Fast Fourier Transform (why `O(N log N)` is achievable)

The Cooley–Tukey FFT computes the 1D DFT of length `n` (n a power of two, for the simplest case) by recursively splitting the sum into even- and odd-indexed terms: ``` X[k] = Σ_{j=0}^{n-1} x[j] · W_n^{jk}, W_n = e^{−i2π/n} = Σ_{j even} x[j] W_n^{jk} + Σ_{j odd} x[j] W_n^{jk} = E[k] + W_n^k · O[k] for k = 0 .. n/2 - 1 X[k + n/2] = E[k] − W_n^k · O[k] (aliasing/symmetry trick) ``` where `E` and `O` are the length-`n/2` DFTs of the even- and odd-indexed samples. ``` function FFT(x): # length n = 2^m n = length(x) if n == 1: return x even = FFT(x[0], x[2], x[4], ...) odd = FFT(x[1], x[3], x[5], ...) X = array of size n for k = 0 to n/2 - 1: twiddle = exp(-2πi k / n) odd[k] X[k] = even[k] + twiddle X[k + n/2] = even[k] - twiddle return X ``` The 2D FFT is obtained by applying the 1D FFT along all rows, then along all columns (separability of the DFT kernel — exactly as convolution with separable filters is applied in Linear Filtering and Convolution).

Time complexity of direct convolution. For an `N × N` image and `k × k` kernel, each output pixel needs `k²` multiply-adds, giving `O(N²k²)`. For small fixed `k` (3, 5, 7) this is effectively linear in image size and is usually fastest in practice.

Time complexity of FFT-based filtering. Computing a 1D FFT of length `n` costs `O(n log n)`: the recursion `T(n) = 2T(n/2) + O(n)` (the `O(n)` being the combine step of the twiddle-factor multiplications) solves, by the Master Theorem, to `T(n) = O(n log n)`. Applying this along `N` rows and `N` columns of an `N × N` image costs `O(N² log N)` total (2D FFT via row–column separability). Multiplying two frequency-domain images pointwise is `O(N²)`, negligible next to the transforms. So end-to-end frequency-domain filtering is `O(N² log N)`, versus `O(N²k²)` spatially.

Crossover point. Frequency-domain filtering wins asymptotically once `k² ≫ log N` — in practice, once kernels exceed roughly 15–30 pixels on a side for typical megapixel images, the FFT route becomes faster, which is why libraries like OpenCV switch strategies (or let the user choose) based on kernel size.

Correctness of the convolution theorem rests on two facts proved above: (1) complex exponentials are eigenfunctions of any LSI system, and (2) the DFT's implicit periodicity means the "convolution" it realizes is circular, not linear — hence the zero-padding requirement to make frequency-domain filtering agree exactly with the direct method on a finite image. With correct zero-padding, the two methods are numerically identical up to floating-point rounding, not merely similar — this is a hard mathematical equality, not an approximation.

Correctness of the FFT algorithm follows by induction: the base case (`n=1`) is trivial (the DFT of a single sample is itself). The inductive step correctly halves the problem because of the periodicity of the twiddle factor `W_n^k`: since `W_n^{k+n/2} = -W_n^k`, the two halves `X[k]` and `X[k+n/2]` can be recovered from the same even/odd sub-transforms with only a sign flip, which is exactly the "butterfly" step in the pseudocode above. No information is lost or duplicated — every one of the `n` inputs contributes to every one of the `n` outputs exactly once, matching the direct `O(n²)` definition, but computed with `O(n log n)` work via `log n` levels of `n/2` butterflies each.

Honest tradeoffs. Frequency-domain filtering assumes the kernel is shift-invariant (the same everywhere in the image); it cannot directly express spatially varying filters (e.g., a blur that increases toward the image edges) without tiling or other tricks. It also requires care with boundary handling (padding introduces its own artifacts, e.g., wrap-around or edge discontinuities that appear as high-frequency ringing — a manifestation of the Digital Images and Sampling issue that a hard, discontinuous edge injects broadband frequency content). Finally, floating-point FFTs accumulate rounding error over `log N` stages, generally negligible for images but relevant in high-precision scientific applications.

Consider a tiny `4 × 4` grayscale image with a single vertical edge — dark on the left, bright on the right — a stand-in for a real photographic edge: ``` f = 10 10 200 200 10 10 200 200 10 10 200 200 10 10 200 200 ```

Step 1 — separability. Every row is identical, so the 2D DFT reduces to a 1D DFT of one row, repeated: the image has zero variation along the vertical axis, so all frequency energy lies on the `u = 0` line (no vertical frequency), and varies only with `v` (horizontal frequency).

Step 2 — 1D DFT of a row `[10, 10, 200, 200]` using `F[v] = Σ_n f[n] e^{-i2πvn/4}`:

Step 3 — interpretation. The DC term (`v=0`, `F=420`) is `4 ×` the row mean `(10+10+200+200)/4 = 105`, confirming that zero frequency captures average brightness. The large magnitude at `v = 1` and `v = 3` (`|F| = √(190²+190²) ≈ 268.7`) reflects the sharp transition — an edge injects energy across many frequencies, especially the fundamental. The zero at `v = 2` (the Nyquist-adjacent frequency for this length) happens because the pattern `10,10,200,200` has no component that alternates every single sample; it is exactly antisymmetric to that basis function's phase, a coincidence of this particular numeric edge width, not a general law.

Step 4 — filtering in frequency. Suppose we apply a simple low-pass filter that zeroes out all but the DC and lowest AC term (`v=1,3`), then inverse-transform. This is equivalent, by the convolution theorem, to convolving the row with a small averaging-like kernel in space. The result is a smoothed transition — the sharp jump from 10 to 200 becomes a gentler ramp — visually and numerically the same operation as applying a blur kernel directly in the spatial domain, verifying the theorem on a case where the input is not trivially symmetric across rows and columns (the image varies richly in one axis and not at all in the other, forcing the reader to track a genuine, non-degenerate 1D DFT rather than an obviously symmetric toy case).

The mathematics underlying this chapter predates computer vision by over a century, developed to solve heat conduction, but was reborn as a computational tool once digital images demanded fast, exact filtering.

JPEG and video compression transform image blocks with the Discrete Cosine Transform (a real-valued cousin of the DFT) precisely because natural images concentrate energy in low frequencies, letting high-frequency coefficients be quantized aggressively with little perceptual loss. Astronomical and medical imaging pipelines use FFT-based deconvolution to remove known blur kernels (telescope point-spread functions, MRI reconstruction from k-space data, which literally is Fourier-domain data). Computational photography uses frequency analysis to design and evaluate deblurring and denoising filters before ever running them on pixels. Texture synthesis and analysis in graphics and vision exploit the magnitude spectrum to characterize repetitive or oriented patterns (wood grain, fabric weave, brick walls). Antialiasing and mipmapping in graphics rendering, and the resampling/warping operations covered next in Image Resampling and Warping, are justified entirely by frequency-domain arguments about which frequencies a given sampling rate can faithfully represent — a direct descendant of the sampling theorem most cleanly proved using the machinery introduced here. Even large convolutional kernels in some numerical libraries and GPU frameworks are executed via FFT-based convolution when kernel sizes are large enough to cross the complexity threshold derived above.

Ask the tutor: Can you show, step by step, why swapping the phase spectra of two different images produces an image that looks like the one whose phase (not magnitude) was used?