The Fourier transform

Every function can be represented as a (maybe infinite) sum of sine and cosine waves with different amplitudes and phases, which is called the frequency spectrum (or frequency response) of the function. This property was discovered by Jean-Baptiste-Joseph Fourier who used it to describe conduction of heat in solid bodies. If a function is represented in that way, i.e., by describing the frequencies and amplitudes, it is called to be depicted in frequency domain, whereas if a function is defined by values at certain positions (which is a representation most people are used to) it is said to be represented in spatial domain.

The representation in frequency domain does not contain additional information, in fact it is just another way of looking at a function, i.e., a dual representation, but some properties of functions are easier understood in frequency domain (``which frequencies are present in the function with which amplitudes'', etc.), others are more obvious in spatial domain.

The Fourier transform(FT) now is a means to switch between spatial and frequency domain (although there exist others, e.g., the Hartley transform [6,18] which is said to be computationally more efficient with real data [8] and often preferred for use in volume visualization algorithms [32,59]). However, as the FT is the traditional, most often facilitated tool to derive spectra we will restrict our discussion to the Fourier transform. Apart from computer graphics it is also used in plasma physics, semiconductor physics, microwave acoustics, seismography and many more fields [7].

To switch from spatial domain to the frequency domain and back, two versions of the Fourier transform are needed, they only differ by a minus sign. Principally, it is not important where this minus sign is put as long as it is kept consistent. The two versions of the Fourier transform are defined by

$$\int_{-\infty}^{\infty} \mathrm{F}(\omega)e^{2\pi j\omega x} \textrm{d}\omega$$ f(x) = (3.7)

$$\mathrm{F}(\omega) \int_{-\infty}^{\infty} \mathrm{f}(x)e^{-2\pi j\omega x} \textrm{d}x$$ = (3.8)

We will use the usual notational conventions here, which are, lowercase letters will denote functions in spatial domain, uppercase letters corresponding functions in frequency domain, and the symbol $\rightleftharpoons$ will be used to assign transform pairs, i.e.,

$$\mathrm{f}(x) \rightleftharpoons \mathrm{F}(\omega)$$ (3.9)

Whenever there is no uppercase letter available for a function name (or it would be confusing, as in Greek letters) the operator $\mathcal{F}$ will be used, so that $\mathcal{F}(\delta)$ denotes the Fourier transform of the impulse signal.

Extending Eqs. 3.7 and 3.8 to higher dimensions is straightforward, e.g., for two dimensional functions ,images for instance, we get

$$\int\!\!\!\int\mathrm{F}(u,v)e^{ 2\pi j(ux + vy)}\,\textrm{d}u\,\textrm{d}v$$ f(x,y) = (3.10)

$$\int\!\!\!\int\mathrm{f}(x,y)e^{-2\pi j(ux + vy)}\,\textrm{d}x\,\textrm{d}y$$ F(u,v) = (3.11)

where u and v are now depicting frequencies. For three variables we would simply have to add another exponent term to the complex exponent and integrate over the additional third dimension. This then can be done similarly for arbitrary dimensions.

From Eqs. 3.7 and 3.8 can be seen that the Fourier transform generates complex output for real data input. The real and complex part of the output specify the magnitude and phase of each frequency component, i.e., sine wave. The magnitude (or amplitude) of $F(\omega)$ is defined by

$$\vert\mathrm{F}(\omega)\vert = \sqrt{[Re(\mathrm{F}(\omega))]^2 + [Im(\mathrm{F}(\omega))]^2}$$ (3.12)

This is often referred to as the frequency response and should not be mistaken with the Fourier spectrum $\mathrm{F}(\omega)$ itself. The phase spectrum is defined by

$$\phi(u) = \tan^{-1}\left( \frac{Im(\mathrm{F}(\omega))}{Re(\mathrm{F}(\omega))} \right )$$ (3.13)

Often, the Fourier transform is plotted by magnitude versus frequency only and thus ignoring the phase angle. Although this may seem problematic it has become conventional, because the bulk of the information of a function is embedded in its frequency content as given by the magnitude spectrum $\vert\mathrm{F}(\omega)\vert$ [64].

The computation of Fourier transforms in practice can often be facilitated by using already well-known transforms of functions (some will be introduced in Sec. 3.3) and applying useful theorems, for instance the addition theorem.

Theorem 1 (Addition Theorem)   The frequency response of the sum of two functions is the sum of the frequency responses of both functions, symbolically,

$$\mathrm{f_1}(x) + \mathrm{f_2}(x) \rightleftharpoons \mathrm{F_1}(\omega) + \mathrm{F_2}(\omega)$$ (3.14)

So far, we only discussed Fourier transforms of continuous functions. However, in practice we often have to deal with discrete functions, that is, a set of samples. Also for this case a Fourier transform pair exists, which is called the discrete Fourier transform (DFT):

$$\sum_{\omega =0}^{N-1} \mathrm{F}[\omega]e^{ 2\pi j\omega x/N}$$ f[x] = (3.15)

$$\mathrm{F}[\omega] \frac{1}{N}\sum_{x=0}^{N-1} \mathrm{f}[x]e^{-2\pi j\omega x/N}$$ = (3.16)

with $0 \leq \omega,x \leq N-1$, and N being the number of samples. This definition sometimes varies slightly by having the coefficient $\frac{1}{N}$ in front of the expression for f[x] and not for $\mathrm{F}[\omega]$, or by having $\frac{1}{\sqrt{N}}$ in front of both expressions. All of these definitions are correct as long as, again, they are kept consistently.

A naive implementation of Eqs. 3.15 and 3.16 results in an algorithm with O(N²) complexity. However, there exists an efficient algorithm which computes the DFT with $O(N \log N)$ complexity, This algorithm is known as the Fast Fourier transform (FFT). There exist a lot of implementations of this algorithm and some of them are available as free software [14,15].

An interesting property of the Fourier transform is its separability. When we extend the DFT (Eqs. 3.15 and 3.16) to two dimension we can rewrite it as

$$\sum_{v=0}^{M-1}\left( \sum_{u=0}^{N-1} \mathrm{F}[u,v]e^{ 2\pi ju x/N} \right)e^{ 2\pi jvy/M}$$ f[x,y] = (3.17)

$$\frac{1}{M}\sum_{y=0}^{M-1}\left( \frac{1}{N}\sum_{x=0}^{N-1} \mathrm{f}[x,y]e^{-2\pi ju x/N} \right)e^{-2\pi jvy/M}$$ F[u,v] = (3.18)

with $0 \leq u,x \leq N-1$ and $0 \leq v,y \leq M-1$. This means that we can compute a multi-dimensional Fourier transform by performing one-dimensional transforms in each dimension consecutively.