f12-4

12.4FFTinTwoorMoreDimensions 515Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).y2=(0.5/wi1)*(y(i)-y(n-i+1)) y(i)=0.5*(y1+y2) y(n-i+1)=0.5*(y1-y2) wtemp=wr1wr1=wr1*wpr-wi1*wpi+wr1 wi1=wi1*wpr+wtemp*wpi+wi1 enddo 16 endifreturn END An alternative way of implementing this algorithm is to form an auxiliary function by copying the even elements of fjinto the first N/2locations, and the odd elements into the next N/2elements in reverse order. However, it is not easy to implement the alternative algorithm without a temporary storage array and we prefer the above in-place algorithm. Finally, we mention that there exist fast cosine transforms for small Nthat do not rely on an auxiliary function or use an FFT routine. Instead, they carry out the transformdirectly, oftencoded in hardwarefor fixed Nof small dimension [1]. CITED REFERENCES AND FURTHER READING: Brigham, E.O. 1974, TheFastFourierTransform (EnglewoodCliffs, NJ: Prentice-Hall), §10–10. Sorensen, H.V., Jones, D.L., Heideman, M.T., and Burris, C.S. 1987, IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. ASSP-35, pp. 849–863. Hou,H.S.1987, IEEETransactionsonAcoustics,Speech,andSignalProcessing ,vol.ASSP-35, pp. 1455–1461 [see for additional references]. Hockney, R.W. 1971,in Methods inComputational Physics , vol. 9 (NewYork: 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., andWintz,P.1987, DigitalImageProcessing ,(Reading,MA:Addison-Wesley). Chen,W.,Smith,C.H.,andFralick,S.C.1977, IEEETransactionsonCommunications ,vol.COM- 25, pp. 1004–1009. [1] 12.4 FFT in Two or More Dimensions Given a complex function h(k1,k2)defined over the two-dimensional grid 0≤k1≤N1−1,0≤k2≤N2−1, we can define its two-dimensional discrete Fouriertransformas a complexfunction H(n1,n2), definedoverthe same grid, H(n1,n2)≡N2−1/summationdisplay k2=0N1−1/summationdisplay k1=0exp(2 πik 2n2/N 2)e x p ( 2 πik 1n1/N 1)h(k1,k2) (12.4.1 ) Bypullingthe“subscripts2”exponentialoutsideofthesumover k1,orbyreversing the order of summation and pulling the “subscripts 1” outside of the sum over k2, 516 Chapter12. FastFourierTransformSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).we can see instantly that the two-dimensionalFFT can be computedby taking one- dimensionalFFTssequentiallyoneachindexoftheoriginalfunction. 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 ) For this to be practical, of course, both N1andN2should be some efficient length for an FFT, usually a power of 2. Programming a two-dimensional FFT, using (12.4.2)with a one-dimensionalFFT 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-dimensionalcomplex array, you find that you are endlessly copying things out of the multidimensional input array and then copying things back into it. This is 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 L- dimensions, is evidently H(n 1,...,n L)≡NL−1/summationdisplay kL=0···N1−1/summationdisplay k1=0exp(2 πik LnL/N L)×··· ×exp(2 πik 1n1/N 1)h(k1,...,k L)(12.4.3 ) where n1andk1range from 0 to N1−1,...,nLandkLrange from 0 to NL−1. How many calls to a one-dimensional FFT are in (12.4.3)? Quite a few! For eachvalueof k 1,k2,...,k L−1youFFT to transformthe Lindex. Thenforeachvalueof k1,k2,...,k L−2andnLyou FFT to transform the L−1index. And so on. It is best to rely on someoneelse havingdone the bookkeepingfor once and for all. The inverse transforms of (12.4.1) or (12.4.3) are just what you would expect them to be: Change the i’s in the exponentials to −i’s, and put an overall factor of 1/(N1×···× NL)in front of the whole thing. Most other features of multidimensional FFTs are also analogous to features already discussed in the one-dimensional case: •Frequencies are arranged in wrap-aroundorder in the transform, but now for each separate dimension. •Theinputdataarealsotreatedasiftheywerewrappedaround. Iftheyare discontinuous across this periodic identification (in any dimension) then the spectrum will have some excess power at high frequencies because of the discontinuity. The fix, if you care, is to remove multidimensional linear trends. •Ifyouaredoingspatialfilteringandareworriedaboutwrap-aroundeffects, then you need to zero-pad all around the border of the multidimensional array. However, be sure to notice how costly zero-padding is in multidi-mensional transforms. If you use too thick a zero-pad, you are going to waste alotof storage, especially in 3 or more dimensions! •Aliasing occurs as always if sufficient bandwidth limiting does not exist along one or more of the dimensions of the transform. 12.4FFTinTwoorMoreDimensions 517Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica)......... ........column of 2 N1 real numbers........ ................ ................ ................ ................ ........... ...N2 columns f2 = 0f2 = 1 N2∆2f2 = N2/2 − 1 N2∆2f2 = ± 1 2∆2f2 = −N2/2 − 1 N2∆2f2 = − 1 N2∆2array element 2N1N2array element 1 col. 1 col. 2 col.N2 2col. +1N2 2col. +2N2 2col. N2 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×N2array h(t1,t2)(stored by columns of complex numbers). The output is also stored by complex columns. Each column corresponds to a particular value of f2,a s showninthe figure. Withineach column,thearrangement offrequencies f1isexactly asshowninFigure 12.2.2 ∆1and∆2are the sampling intervals in the 1 and 2 directions, respectively. The total number of (real) array elements is 2N1N2. The program fourncan also do more than two dimensions, and the storage arrangement generalizes in the obvious way. Theroutine fournthatwefurnishherewithisadescendantofonewrittenbyN. 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). Note that these lengths must allbe powers of 2, and are the numbers of complexvalues in each direction; (iii) the usual scalar equal to ±1indicating whether you want the transform or its inverse; and, finally (iv) the array of data. A few words about the data array: fournaccesses it as a one-dimensional array of real numbers, of length equal to twice the product of the lengths of the Ldimensions. It assumes that the array represents an L-dimensional complex array, in normal FORTRAN order. Normal FORTRAN order means: (i) each complex value occupies two sequential locations, real part followed by imaginary; (ii) thefirst subscript changes most rapidly as one goes through the array; the last subscriptchangesleastrapidly;(iii)subscriptsrangefrom1totheirmaximumvalues(N 1,N 2,...,N L,respectively),ratherthanfrom0to N1−1,N 2−1,..., N L−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.) 518 Chapter12. FastFourierTransformSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE fourn(data,nn,ndim,isign) INTEGER isign,ndim,nn(ndim) REAL data(*) 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 complex values), which MUST all be powers of 2. data is a real array of length twice the product of these lengths, in which the data are stored as in a multidimensional complex FORTRAN array. If isign is input as −1,data is replaced by its inverse transform times the product of the lengths of all dimensions. INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2, * ip3,k1,k2,n,nprev,nrem,ntot REAL tempi,temprDOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp Double precision for trigonometric re- currences. ntot=1 do 11idim=1,ndim Compute total number of complex values. ntot=ntot*nn(idim) enddo 11 nprev=1 do18idim=1,ndim Main loop over the dimensions. n=nn(idim) nrem=ntot/(n*nprev) ip1=2*nprevip2=ip1*n ip3=ip2*nrem i2rev=1do 14i2=1,ip2,ip1 This is the bit-reversal section of the routine. if(i2.lt.i2rev)then do13i1=i2,i2+ip1-2,2 do12i3=i1,ip3,ip2 i3rev=i2rev+i3-i2 tempr=data(i3) tempi=data(i3+1)data(i3)=data(i3rev)data(i3+1)=data(i3rev+1) data(i3rev)=tempr data(i3rev+1)=tempi enddo 12 enddo 13 endifibit=ip2/2 1 if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then i2rev=i2rev-ibit ibit=ibit/2 goto 1endif i2rev=i2rev+ibit enddo 14 ifp1=ip1 Here begins the Danielson-Lanczos section of the routine. 2 if(ifp1.lt.ip2)then ifp2=2*ifp1theta=isign*6.28318530717959d0/(ifp2/ip1) Initialize for the trig. recur- rence. wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0wi=0.d0 do 17i3=1,ifp1,ip1 do16i1=i3,i3+ip1-2,2 do15i2=i1,ip3,ifp2 k1=i2 Danielson-Lanczos formula: k2=k1+ifp1 tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)data(k2)=data(k1)-tempr data(k2+1)=data(k1+1)-tempi 12.5FourierTransformsofRealDatainTwoandThreeDimensions 519Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).data(k1)=data(k1)+tempr data(k1+1)=data(k1+1)+tempi enddo 15 enddo 16 wtemp=wr Trigonometric recurrence. wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi enddo 17 ifp1=ifp2 goto 2 endifnprev=n*nprev enddo 18 returnEND CITED REFERENCES AND FURTHER READING: Nussbaumer,H.J.1982, FastFourierTransformandConvolutionAlgorithms (NewYork:Springer- Verlag). 12.5 Fourier Transforms of Real Data in Two and Three Dimensions Two-dimensionalFFTs areparticularlyimportantinthe fieldofimageprocess- ing. Animageis usuallyrepresentedas atwo-dimensionalarrayofpixelintensities, real (and usually positive) numbers. One commonly desires to filter high, or low, frequency spatial components from an image; or to convolve or deconvolve the image 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 equation for a potential (e.g., electromagneticor gravitational) on a three-dimensionallatticethat represents the discretization of three-dimensionalspace. Here the source terms (mass or charge distribution) and the desired potentials are also real. In two and three dimensions, with large arrays, memory is often at a premium. It is therefore important to perform the FFTs, insofar as possible, on the data “in place.”We wantaroutinewithfunctionalitysimilartothemultidimensionalFFTroutine fourn (§12.4), but which operates on real, not complex, input data. We give such a routinein this section. The developmentis analogousto that of §12.3 leadingto the one-dimensional routine realft. (You might wish to review that material at this point, particularly equation 12.3.5.) It is convenient to think of the independent variables n 1,...,n Lin equation (12.4.3)asrepresentingan L-dimensionalvector /vectorninwave-numberspace,withvalues on the lattice of integers. The transform H(n1,...,n L)is then denoted H(/vectorn). Itiseasytoseethatthetransform H(/vectorn)isperiodicineachofits Ldimensions. Specifically, if /vectorP1,/vectorP2,/vectorP3,...denote the vectors (N1,0,0,... ),(0,N 2,0,... ), (0,0,N 3,... ), and so forth, then H(/vectorn±/vectorPj)=H(/vectorn) j=1,...,L (12.5.1 )