(p.187) A Discrete Fourier transforms
(p.187) A Discrete Fourier transforms
The discrete nature of computer technology means that the continuum formulation of Fourier transforms in Section 2.4 is best replaced with a digitized version for numerical calculations. In the simplest case, this takes the form
where the integers j and k = 0, 1, 2, …, N −1 and the f j represent a sampling of the function f(x) on a grid of N points that are spaced evenly along the x-axis, and
With the complementary complex exponentials in the summations, these clearly have the essential ingredients to be a Fourier pair. To see how they relate to the discussion in Chapter 2, we need to examine them more carefully.
The first thing to note is that both f j and F k are implicitly periodic, with a repeat unit of N pixels:
from eqn (A.2), since
(p.188) as illustrated in Fig. A.1; likewise,
from eqn (A.1).
Although the periodicity of f j means that an independent interval of N consecutive pixels can be chosen anywhere along the x-axis, the location assigned to the point j=0 is highly significant: as in the continuum case, it defines the origin. For example, the substitution of a δ-function at j=0,
into eqn (A.1) yields F k = 1 for all k, as would be expected from eqn (2.46) when f(x)=δ(x). By contrast, a displacement of the non-zero pixel to j ≠ 0 leads to the F k being complex: a uniform modulus, but an argument that increases linearly with k. The corresponding position of the origin in reciprocal space is easiest to verify by putting k=0 in eqn(A.1). It gives
which is the discrete analogue of eqn (2.49).
With the periodicity of eqn (A.3), and the origin at j = 0, an even function is one that satisfies
Hence, the symmetry relationships of eqn (2.48) become
where we have also made use of eqn (A.4). Similarly, the discrete form of eqn (2.48) is
This implies that, for a real function f j , F0 is real:
Thus the Nyquist component, F N/2 , must also be real.
(p.189) A consideration of Fig. A.1 shows that the location of pixel j corresponds to a value of x given by
where x max is the length of the repeat unit. While x max is chosen at our discretion, it has to be bigger than the (finite) extent of the continuous function, f(x), for the calculation to be meaningful. If a convolution is to be performed with g(x), say, then x max needs to be larger than the combined spread of f(x) ⊗ g(x) to ensure that aliasing artefacts, resulting from the inherent periodicity of the sum, are avoided. The equivalence between the integer k and the related continuous variable, q, can also be expressed in the manner of eqn (A.8):
where the reciprocity of Fourier length scales dictates that
The constant of proportionality is 2π if the continuum exponential is of the form exp(±iqx); it’s unity for exp(±i2πqx). Equation (A.10) requires qualification, in that q max must be at least twice the highest ‘measured’ frequency: due to eqn (A.4), the Fourier coefficients pertaining to
or |q| ≤ q max/2, define the F k spectrum completely. According to eqn (A.7), only the components with 0 ≤ k ≤ N/2 need be specified if the function f j is real.
The discrete Fourier transform (DFT) of eqns (A.1) and (A.2) can (algorithm. Whereas the direct computation of a DFT entails 𝒪 (N 2) operations, because each of the F k involves the sum of N products, the same result can be obtained with just 𝒪(N log2 N) calculations by using the clever factorization of an FFT. The associated gain in computational speed, of 𝒪(N/log2 N), is enormous even for modest N. The only thing to bear in mind is that, for the most common FFT algorithm (Cooley and Tukey, 1965 ), N has to be a power-of-two:
where M is a positive integer. This requirement can be fulfilled by ‘padding out’ the f j function, or the F k spectrum, with zeros in an appropriate fashion. A suitable procedure for the example of Fig. A.1 is illustrated in Fig. A.2: with the increment in x, Δx, unchanged between the j-pixels, q max = 1/Δx remains the same; the larger N, however, means that the F k sample the Fourier spectrum much more finely in q.
where the integers j and j ′ relate to x and y, k and k ′ correspond to q x and q y , and the digitization is on a grid of N ×N ′.