c12-4 (Numerical Recipes in C)

12.4 FFT in Two or More Dimensions521}}}CITED REFERENCES AND FURTHER READING:Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall), §10–10.Sorensen, H.V., Jones, D.L., Heideman, M.T., and Burris, C.S. 1987, IEEE Transactions onAcoustics, Speech, and Signal Processing, vol.

ASSP-35, pp. 849–863.Hou, H.S. 1987, IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35,pp. 1455–1461 [see for additional references].Hockney, R.W. 1971, in Methods in Computational Physics, vol. 9 (New York: Academic Press).Temperton, C. 1980, Journal of Computational Physics, vol. 34, pp. 314–329.Clarke, R.J. 1985, Transform Coding of Images, (Reading, MA: Addison-Wesley).Gonzalez, R.C., and Wintz, P.

1987, Digital Image Processing, (Reading, MA: Addison-Wesley).Chen, W., Smith, C.H., and Fralick, S.C. 1977, IEEE Transactions on Communications, vol. COM25, pp. 1004–1009. [1]12.4 FFT in Two or More DimensionsGiven a complex function h(k1 , k2 ) defined over the two-dimensional grid0 ≤ k1 ≤ N1 − 1, 0 ≤ k2 ≤ N2 − 1, we can define its two-dimensional discreteFourier transform as a complex function H(n1 , n2 ), defined over the same grid,H(n1 , n2 ) ≡NX2 −1 N1 −1Xexp(2πik2 n2 /N2 ) exp(2πik1 n1 /N1 ) h(k1 , k2 )k2 =0 k1 =0(12.4.1)By pulling the “subscripts 2” exponential outside of the sum over k1 , or by reversingthe order of summation and pulling the “subscripts 1” outside of the sum over k2 ,we can see instantly that the two-dimensional FFT can be computed by taking onedimensional FFTs sequentially on each index of the original function.

Symbolically,H(n1 , n2 ) = FFT-on-index-1 (FFT-on-index-2 [h(k1 , k2)])= FFT-on-index-2 (FFT-on-index-1 [h(k1 , k2)])(12.4.2)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.

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).An alternative way of implementing this algorithm is to form an auxiliaryfunction by copying the even elements of fj into the first N/2 locations, and theodd elements into the next N/2 elements in reverse order. However, it is not easyto implement the alternative algorithm without a temporary storage array and weprefer the above in-place algorithm.Finally, we mention that there exist fast cosine transforms for small N that donot rely on an auxiliary function or use an FFT routine.

Instead, they carry out thetransform directly, often coded in hardware for fixed N of small dimension [1].522Chapter 12.Fast Fourier TransformH(n1 , . . . , nL ) ≡NXL −1kL =0···NX1 −1exp(2πikL nL /NL ) × · · ·k1 =0(12.4.3)× exp(2πik1 n1 /N1 ) h(k1 , . . . , kL )where n1 and k1 range from 0 to N1 − 1, .

. . , nL and kL range from 0 to NL − 1.How many calls to a one-dimensional FFT are in (12.4.3)? Quite a few! For eachvalue of k1 , k2 , . . . , kL−1 you FFT to transform the L index. Then for each value ofk1 , k2 , . . . , kL−2 and nL you FFT to transform the L − 1 index. And so on.

It isbest to rely on someone else having done the bookkeeping for once and for all.The inverse transforms of (12.4.1) or (12.4.3) are just what you would expectthem to be: Change the i’s in the exponentials to −i’s, and put an overallfactor of 1/(N1 × · · · × NL ) in front of the whole thing. Most other featuresof multidimensional FFTs are also analogous to features already discussed in theone-dimensional case:• Frequencies are arranged in wrap-around order in the transform, but nowfor each separate dimension.• The input data are also treated as if they were wrapped around. If they arediscontinuous across this periodic identification (in any dimension) thenthe spectrum will have some excess power at high frequencies becauseof the discontinuity. The fix, if you care, is to remove multidimensionallinear trends.• If you are doing spatial filtering and are worried about wrap-around effects,then you need to zero-pad all around the border of the multidimensionalarray.

However, be sure to notice how costly zero-padding is in multidimensional transforms. If you use too thick a zero-pad, you are going towaste a lot of storage, especially in 3 or more dimensions!• Aliasing occurs as always if sufficient bandwidth limiting does not existalong one or more of the dimensions of the transform.The routine fourn that we furnish herewith is a descendant of one written by N.M. Brenner. It requires as input (i) a scalar, telling the number of dimensions, e.g.,2; (ii) a vector, telling the length of the array in each dimension, e.g., (32,64). Notethat these lengths must all be powers of 2, and are the numbers of complex valuesin each direction; (iii) the usual scalar equal to ±1 indicating whether you want thetransform or its inverse; and, finally (iv) the array of data.A few words about the data array: fourn accesses it as a one-dimensionalarray of real numbers, that is, data[1..(2N1N2 .

. . NL )], of length equal to twiceSample 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).For this to be practical, of course, both N1 and N2 should be some efficient lengthfor an FFT, usually a power of 2.

Programming a two-dimensional FFT, using(12.4.2) with a one-dimensional FFT routine, is a bit clumsier than it seems at first.Because the one-dimensional routine requires that its input be in consecutive orderas a one-dimensional complex array, you find that you are endlessly copying thingsout of the multidimensional input array and then copying things back into it. Thisis not recommended technique. Rather, you should use a multidimensional FFTroutine, such as the one we give below.The generalization of (12.4.1) to more than two dimensions, say to Ldimensions, is evidently52312.4 FFT in Two or More Dimensionsdata [1]row 1f1 = 0row 2f1 =row N1 / 2f1 =1N1 ∆11⁄ 2 Nrow N1 / 2 + 1f1 = ±row N1 / 2 + 2f1 = −row N1f1 = −1−1N1 ∆112∆11⁄ 2 N1−1N1 ∆11N1 ∆1data [2N1N2 ]Figure 12.4.1.

Storage arrangement of frequencies in the output H(f1, f2 ) of a two-dimensional FFT.The input data is a two-dimensional N1 × N2 array h(t1 , t2 ) (stored by rows of complex numbers).The output is also stored by complex rows. Each row corresponds to a particular value of f1 , as shownin the figure. Within each row, the arrangement of frequencies f2 is exactly as shown in Figure 12.2.2.∆1 and ∆2 are the sampling intervals in the 1 and 2 directions, respectively.

The total number of (real)array elements is 2N1 N2 . The program fourn can also do more than two dimensions, and the storagearrangement generalizes in the obvious way.the product of the lengths of the L dimensions. It assumes that the array representsan L-dimensional complex array, with individual components ordered as follows:(i) each complex value occupies two sequential locations, real part followed byimaginary; (ii) the first subscript changes least rapidly as one goes through the array;the last subscript changes most rapidly (that is, “store by rows,” the C norm); (iii)subscripts range from 1 to their maximum values (N1 , N2 , .

. . , NL , respectively),rather than from 0 to N1 − 1, N2 − 1, . . . , NL − 1. Almost all failures to get fournto work result from improper understanding of the above ordering of the data array,so take care! (Figure 12.4.1 illustrates the format of the output array.)#include <math.h>#define SWAP(a,b) tempr=(a);(a)=(b);(b)=temprvoid fourn(float data[], unsigned long nn[], int ndim, int isign)Replaces data by its ndim-dimensional discrete Fourier transform, if isign is input as 1.nn[1..ndim] is an integer array containing the lengths of each dimension (number of complexvalues), which MUST all be powers of 2. data is a real array of length twice the product ofthese lengths, in which the data are stored as in a multidimensional complex array: real andimaginary parts of each element are in consecutive locations, and the rightmost index of thearray increases most rapidly as one proceeds along data.