Jump to ContentJump to Main Navigation
Elementary Scattering TheoryFor X-ray and Neutron Users$

D.S. Sivia

Print publication date: 2011

Print ISBN-13: 9780199228676

Published to Oxford Scholarship Online: December 2013

DOI: 10.1093/acprof:oso/9780199228676.001.0001

Show Summary Details
Page of

PRINTED FROM OXFORD SCHOLARSHIP ONLINE (www.oxfordscholarship.com). (c) Copyright Oxford University Press, 2017. All Rights Reserved. Under the terms of the licence agreement, an individual user may print out a PDF of a single chapter of a monograph in OSO for personal use (for details see http://www.oxfordscholarship.com/page/privacy-policy). Subscriber: null; date: 26 February 2017

(p.187) A Discrete Fourier transforms

(p.187) A Discrete Fourier transforms

Elementary Scattering Theory
Oxford University Press

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

F k = j = 0 N 1 f j e i 2 π j k / N ,

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

f k = 1 N k = 0 N 1 F k e i 2 π j k / N .

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:

f j = f j + N

from eqn (A.2), since

e i2 π ( j + N ) k / N = e i2 π j k / N e i2 π k = e i2 π j k / N

                     A Discrete Fourier transforms

Fig. A.1 A continuous and periodic function sampled at N points per repeat unit.

(p.188) as illustrated in Fig. A.1; likewise,

F k = F k + N

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,

f j = { 1 if j = 0 , 0 otherwise,

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

F 0 = j = 0 N 1 f j ,

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

f j = f j = f N j .

Hence, the symmetry relationships of eqn (2.48) become

f j = { f N j f N j F k = { F N k , F N k ,

where we have also made use of eqn (A.4). Similarly, the discrete form of eqn (2.48) is

f j = f j * F N k = F k * .

This implies that, for a real function f j , F0 is real:

F 0 = F N = F 0 * ,

which is confirmed by eqn (A.5). Assuming that N is an even number, and putting k=N/2 in eqn (A.7),

F N / 2 = F N / 2 * .

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

j x j = x max N j ,

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):

k q k = q max N k ,

where the reciprocity of Fourier length scales dictates that

q max N x max .

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

| k | N / 2 ,

or |q| ≤ q max/2, define the F k spectrum completely. According to eqn (A.7), only the components with 0 ≤ kN/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:

N = 2 M ,

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.


                     A Discrete Fourier transforms

Fig. A.2 The function in Fig. A.1 zero-padded to make N a power-of-two.

The DFT of eqns (A.1) and (A.2), and its FFT implementation, can be generalized easily to many dimensions. Explicitly, the case of two takes the form

F k k = j = 0 N 1 j = 0 N 1 f j j e i 2 π ( j k / N + j k / N )


f j j = 1 N N k = 0 N 1 k = 0 N 1 F k k e i 2 π ( j k / N + j k / N )

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 .