c12-5 (Numerical Recipes in C)

52512.5 Fourier Transforms of Real Data in Two and Three DimensionsCITED REFERENCES AND FURTHER READING:Nussbaumer, H.J. 1982, Fast Fourier Transform and Convolution Algorithms (New York: SpringerVerlag).Two-dimensional FFTs are particularly important in the field of image processing. An image is usually represented as a two-dimensional array of pixel intensities,real (and usually positive) numbers.

One commonly desires to filter high, or low,frequency spatial components from an image; or to convolve or deconvolve theimage with some instrumental point spread function. Use of the FFT is by far themost efficient technique.In three dimensions, a common use of the FFT is to solve Poisson’s equationfor a potential (e.g., electromagnetic or gravitational) on a three-dimensional latticethat represents the discretization of three-dimensional space. Here the source terms(mass or charge distribution) and the desired potentials are also real.

In two andthree dimensions, with large arrays, memory is often at a premium. It is thereforeimportant to perform the FFTs, insofar as possible, on the data “in place.” Wewant a routine with functionality similar to the multidimensional FFT routine fourn(§12.4), but which operates on real, not complex, input data. We give such aroutine in this section. The development is analogous to that of §12.3 leading tothe one-dimensional routine realft.

(You might wish to review that material atthis point, particularly equation 12.3.5.)It is convenient to think of the independent variables n1 , . . . , nL in equation(12.4.3) as representing an L-dimensional vector ~n in wave-number space, withvalues on the lattice of integers. The transform H(n1 , . . . , nL ) is then denotedH(~n).It is easy to see that the transform H(~n) is periodic in each of its L dimensions.~2, P~ 3, . .

. denote the vectors (N1 , 0, 0, . . .), (0, N2 , 0, . . .),Specifically, if P~1 , P(0, 0, N3 , . . .), and so forth, thenH(~n ± P~j ) = H(~n)j = 1, . . . , L(12.5.1)Equation (12.5.1) holds for any input data, real or complex. When the data is real,we have the additional symmetryH(−~n) = H(~n)*(12.5.2)Equations (12.5.1) and (12.5.2) imply that the full transform can be trivially obtainedfrom the subset of lattice values ~n that have0 ≤ n1 ≤ N1 − 10 ≤ n2 ≤ N2 − 1(12.5.3)···0 ≤ nL ≤NL2Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).12.5 Fourier Transforms of Real Data in Twoand Three Dimensions526Chapter 12.Fast Fourier TransformRe(SPEC[i1][i2][i3]) = data[i1][i2][2*i3-1]Im(SPEC[i1][i2][i3]) = data[i1][i2][2*i3](12.5.4)The remaining “plane” of values, SPEC[1..nn1][1..nn2][nn3/2+1], is returnedin the two-dimensional float array speq[1..nn1][1..2*nn2],with the correspondenceRe(SPEC[i1][i2][nn3/2+1]) = speq[i1][2*i2-1]Im(SPEC[i1][i2][nn3/2+1]) = speq[i1][2*i2](12.5.5)Note that speq contains frequency components whose third component f3 is atthe Nyquist critical frequency ±fc .

In some applications these values will in factbe ignored or set to zero, since they are intrinsically aliased between positive andnegative frequencies.With this much introduction, the implementing procedure, called rlft3, issomething of an anticlimax. Look in the innermost loop in the procedure, and youwill see equation (12.3.5) implemented on the last transform index. The case ofi3=1 is coded separately, to account for the fact that speq is to be filled instead ofSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).In fact, this set of values is overcomplete, because there are additional symmetryrelations among the transform values that have nL = 0 and nL = NL /2.

Howeverthese symmetries are complicated and their use becomes extremely confusing.Therefore, we will compute our FFT on the lattice subset of equation (12.5.3),even though this requires a small amount of extra storage for the answer, i.e., thetransform is not quite “in place.” (Although an in-place transform is in fact possible,we have found it virtually impossible to explain to any user how to unscramble itsoutput, i.e., where to find the real and imaginary components of the transform atsome particular frequency!)We will implement the multidimensional real Fourier transform for the threedimensional case L = 3, with the input data stored as a real, three-dimensional arraydata[1..nn1][1..nn2][1..nn3].

This scheme will allow two-dimensional datato be processed with effectively no loss of efficiency simply by choosing nn1 = 1.(Note that it must be the first dimension that is set to 1.) The output spectrum comesback packaged, logically at least, as a complex, three-dimensional array that we cancall SPEC[1..nn1][1..nn2][1..nn3/2+1] (cf. equation 12.5.3).

In the first twoof its three dimensions, the respective frequency values f1 or f2 are stored in wraparound order, that is with zero frequency in the first index value, the smallest positivefrequency in the second index value, the smallest negative frequency in the lastindex value, and so on (cf. the discussion leading up to routines four1 and fourn).The third of the three dimensions returns only the positive half of the frequencyspectrum. Figure 12.5.1 shows the logical storage scheme. The returned portion ofthe complex output spectrum is shown as the unshaded part of the lower figure.The physical, as opposed to logical, packaging of the output spectrum is necessarily a bit different from the logical packaging, because C does not have a convenient,portable mechanism for equivalencing real and complex arrays.

The subscript rangeSPEC[1..nn1][1..nn2][1..nn3/2] is returned in the input array data[1..nn1][1..nn2][1..nn3], with the correspondence52712.5 Fourier Transforms of Real Data in Two and Three Dimensionsnn2, 1nn2, nn3Input data arrayOutput spectrumarrays (complex)1, nn3f3 = 0f 3 = – fcf2 = fc1,1f2 = 0returned in spec[1..nn1][1..2*nn2]1,1nn2,nn3/2returned indata[1..nn1][1..nn2][1..nn3]nn2,11,nn3/2f2 = − fcFigure 12.5.1.Input and output data arrangement for rlft3. All arrays shown are presumedto have a first (leftmost) dimension of range [1..nn1], coming out of the page. The input dataarray is a real, three-dimensional array data[1..nn1][1..nn2][1..nn3].

(For two-dimensionaldata, one sets nn1 = 1.) The output data can be viewed as a single complex array with dimensions[1..nn1][1..nn2][1..nn3/2+1] (cf. equation 12.5.3), corresponding to the frequency componentsf1 and f2 being stored in wrap-around order, but only positive f3 values being stored (others beingobtainable by symmetry).

The output data is actually returned mostly in the input array data, but partlystored in the real array speq[1..nn1][1..2*nn2]. See text for details.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.