f12-5

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 areparticularlyimportantinthefieldofimageprocess- 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 ) 520 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).Equation (12.5.1) holds for any input data, real or complex. When the data is real, we have the additional symmetry H(−/vectorn)=H(/vectorn)* ( 12.5.2 ) Equations(12.5.1)and(12.5.2)implythatthefulltransformcanbetriviallyobtained from the subset of lattice values /vectornthat have 0≤n1≤N1 2 0≤n2≤N2−1 ··· 0≤nL≤NL−1(12.5.3 ) In fact, this set of values is overcomplete, because there are additional symmetry relations among the transform values that have n1=0andn1=N1/2. However these 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., thetransformisnot quite“inplace.” (Althoughanin-placetransformisinfactpossible, we have found it virtually impossible to explain to any user how to unscramble its output, i.e., where to find the real and imaginary components of the transform at some particular frequency!) Figure 12.5.1 shows the storage scheme that we will use for the input data and the output transform. The figure is specialized to the case of two dimensions, L=2, but the generalization to higher dimensions is obvious. The input data is a two-dimensional real array of dimensions N 1(called nn1)b y N2(called nn2). Noticethatthe FORTRAN subscriptsnumberfrom1to nn1,andnotfrom0to N1−1. The output spectrum is in two complex arrays, one two-dimensional and the other one-dimensional. The two-dimensional one, spec, has dimensions nn1/2bynn2. This is exactly half the size of the input data array; but since it is complex, it is the same amount of storage. In fact, specwill share storage with (and overwrite) the input data array. As the figure shows, speccontains those spectral components whose first component of frequency, f1, ranges from zero to just short of the Nyquist frequency fc. The full rangeof positive andnegativesecond-componentof frequencies, f2,isstored,inwrap-aroundorder(see §12.2),withnegativefrequencies shifted by exactly one period to put them “above” the positive frequencies, as the figure indicates. The figure also indicates how the additional L−1(here, one-) dimensional array speqstores only that single value of n1that corresponds to the Nyquist frequency, but all values of n2, etc. With this much introduction, the implementing procedure, called rlft3,i s somethingofananticlimax. Theroutineiswrittenforthecaseof L=3dimensions, but (we will explain below)it can be used without modificationfor L=2also; and it is quite trivial to generalizeit to larger L. Look at the innermost(“do 13”) loop in theprocedure,andyouwillseeequation(12.3.5)implementedonthe firsttransform index. The case of i1=1is coded separately, to account for the fact that speqis to be filled instead of spec(which is here called datasince it shares storage with 12.5FourierTransformsofRealDatainTwoandThreeDimensions 521Sample 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).Input data array Output spectrum arraysnn1, 1 nn1, nn2 REAL data(nn1,nn2) 1,1 1,11, nn2 nn1/2,1 nn1/2,nn2 1,nn2COMPLEX speq(nn2) COMPLEX spec(nn1/2,nn2)f1 = fc f1 = 0 f2 = fcf2 = 01 nn2f2 = –fc f1 = –fc Figure 12.5.1. Input and output data arrangement for rlft3in the case of two-dimensional data. The inputdataarrayisareal,two-dimensionalarray. Theoutputdataarray specisacomplex,two-dimensional array whose (1,1)element contains the f1=f2=0spectral component; a complete set of f2values are stored in wrap-around order, while only positive f1values are stored (others being obtainable by symmetry). The output array speqcontains components with f1equal to the Nyquist frequency. 522 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).the input array). The three enclosing doloops (indices i2,i1, and i3, from inside to outside) could in fact be done in any order —their actions all commute. We chosetheordershownbecauseofthefollowingconsiderations: (i) i1shouldnotbe the inner loop, because if it is, then the recurrence relations on wrandwibecome burdensome. (ii)Onvirtual-memorymachines, i3shouldbetheouterloop,because (with FORTRAN order of array storage) this results in the array data, which might be very large, being accessed in block sequential order. Notethattheworkdonein rlft3is quite(logarithmically)small,comparedto the associated complexFFT, fourn. For this reason, we allow ourselves the clarity of using FORTRAN complex arithmetic even when (as in the multiplications by c1 andc2) there are a few unnecessary operations. The routine rlft3is based on an earlier routine by G.B. Rybicki. SUBROUTINE rlft3(data,speq,nn1,nn2,nn3,isign) INTEGER isign,nn1,nn2,nn3COMPLEX data(nn1/2,nn2,nn3),speq(nn2,nn3) C USES fourn Given a two- or three-dimensional real array data whose dimensions are nn1,nn2,nn3 (where nn3 is 1 for the case of a two-dimensional array), this routine returns (for isign=1 ) the complex fast Fourier transform as two complex arrays: On output, data contains the zero and positive frequency values of the first frequency component, while speq contains the Nyquist critical frequency values of the first frequency component. Second (and third)frequency components are stored for zero, positive, and negative frequencies, in standardwrap-around order. For isign=-1 , the inverse transform (times nn1*nn2*nn3/2 as a constant multiplicative factor) is performed, with output data (viewed as a real array) deriving from input data (viewed as complex) and speq . For inverse transforms on data not generated first by a forward transform, make sure the complex input data array satisfies property (12.5.2). The dimensions nn1,nn2,nn3 must always be integer powers of 2. INTEGER i1,i2,i3,j1,j2,j3,nn(3)DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtempCOMPLEX c1,c2,h1,h2,w Note that data is dimensioned as complex , its output format. c1=cmplx(0.5,0.0) c2=cmplx(0.0,-0.5*isign)theta=6.28318530717959d0/dble(isign*nn1) wpr=-2.0d0*sin(0.5d0*theta)**2 wpi=sin(theta)nn(1)=nn1/2nn(2)=nn2 nn(3)=nn3 if(isign.eq.1)then Case of forward transform. call fourn(data,nn,3,isign) Here is where most all of the compute time is spent. do 12i3=1,nn3 Extend data periodically into speq . do11i2=1,nn2 speq(i2,i3)=data(1,i2,i3) enddo 11 enddo 12 endifdo 15i3=1,nn3 j3=1 Zero frequency is its own reflection, otherwise locate cor- responding negative frequency in wrap-around order. if (i3.ne.1) j3=nn3-i3+2 wr=1.0d0 Initialize trigonometric recurrence. wi=0.0d0 do14i1=1,nn1/4+1 j1=nn1/2-i1+2do 13i2=1,nn2 j2=1 if (i2.ne.1) j2=nn2-i2+2 if(i1.eq.1)then Equation (12.3.5). h1=c1*(data(1,i2,i3)+conjg(speq(j2,j3))) h2=c2*(data(1,i2,i3)-conjg(speq(j2,j3))) 12.5FourierTransformsofRealDatainTwoandThreeDimensions 523Sample 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).Figure 12.5.2. (a) Atwo-dimensional image with intensities either purely black or purely white. (b)The sameimage,afterithasbeenlow-pass filteredusing rlft3. Regionswith fine-scalefeaturesbecomegray. data(1,i2,i3)=h1+h2 speq(j2,j3)=conjg(h1-h2) else h1=c1*(data(i1,i2,i3)+conjg(data(j1,j2,j3)))h2=c2*(data(i1,i2,i3)-conjg(data(j1,j2,j3))) data(i1,i2,i3)=h1+w*h2 data(j1,j2,j3)=conjg(h1-w*h2) endif enddo 13 wtemp=wr Do the recurrence. wr=wr*wpr-wi*wpi+wrwi=wi*wpr+wtemp*wpi+wiw=cmplx(sngl(wr),sngl(wi)) enddo 14 enddo 15 if(isign.eq.-1)then Case of reverse transform. call fourn(data,nn,3,isign) endifreturnEND We now give some fragments from notional calling programs, to clarify the use of rlft3for two- and three-dimensional data. Note that the routine does not actuallydistinguishbetweentwoand threedimensions;twois treatedlike three,but with the third dimension having length 1. Since the third dimension is the outerloop, almost no inef ficiency is introduced. ThefirstprogramfragmentFFTsatwo-dimensionaldataarray,allowsforsome processing on it, e.g., filtering, and then takes the inverse transform. Figure 12.5.2 shows an example of the use of this kind of code: A sharp image becomes blurry when its high-frequency spatial components are suppressed by the factor (here)max (1−6f 2/f2 c,0). The second programexample illustrates a three-dimensional transform, where the three dimensions have different lengths. The third program exampleisanexampleofconvolution,asitmightoccurinaprogramtocomputethepotential generated by a three-dimensional distribution of sources. 524 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).PROGRAM exmpl1 This fragment shows how one might filter a 256 by 256 digital image. INTEGER N1,N2,N3 PARAMETER (N1=256,N2=256,N3=1) Note that the third component must be set to 1. C USES rlft3 REAL data(N1,N2) COMPLEX spec(N1/2,N2),speq(N2)EQUIVALENCE (data,spec) C ... Here the image would be loaded into data . call rlft3(data,speq,N1,N2,N3,1) C ... Here the arrays spec andspeq would be multiplied by a suit- able filter function (of frequency). call rlft3(data,speq,N1,N2,N3,-1) C ... Here the filtered image would be unloaded from data . END PROGRAM exmpl2 This fragment shows how one might FFT a real three-dimensional array of size 32 by 64 by 16. INTEGER N1,N2,N3PARAMETER (N1=32,N2=64,N3=16) C USES rlft3 REAL data(N1,N2,N3)COMPLEX spec(N1/2,N2,N3),speq(N2,N3) EQUIVALENCE (data,spec) C ... Here load data . call rlft3(data,speq,N1,N2,N3,1) C ... Here unload spec and speq . END PROGRAM exmpl3 This fragment shows how one might convolve two real, three-dimensional arrays of size 32 by 32 by 32, replacing the first array by the result. INTEGER NPARAMETER (N=32) C USES rlft3 INTEGER jREAL fac,data1(N,N,N),data2(N,N,N)COMPLEX spec1(N/2,N,N),speq1(N,N),spec2(N/2,N,N),speq2(N,N), * zpec1(N*N*N/2),zpeq1(N*N),zpec2(N*N*N/2),zpeq2(N*N) EQUIVALENCE (data1,spec1,zpec1), (data2,spec2,zpec2), * (speq1,zpeq1), (speq2,zpeq2) C ... call rlft3(data1,speq1,N,N,N,1) FFT both input arrays. call rlft3(data2,speq2,N,N,N,1)fac=2./(N*N*N) Factor needed to get normalized inverse. do 11j=1,N*N*N/2 The sole purpose of the zpec sa n d zpeq si st om a k e this a single do-loop instead of three-nested ones. zpec1(j)=fac*zpec1(j)*zpec2(j) enddo 11 do12j=1,N*N zpeq1(j)=fac*zpeq1(j)*zpeq2(j) enddo 12 call rlft3(data1,speq1,N,N,N,-1) Inverse FFT the product of the two FFTs. C ... END Toextend rlft3tofourdimensions,yousimplyaddanadditional(outer)nested doloopin i4,analogoustothepresent i3. (Modifyingtheroutinetodoan arbitrary numberofdimensions,as in fourn,is agoodprogrammingexerciseforthereader.) CITED REFERENCES AND FURTHER READING: Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall). Swartztrauber, P. N. 1986, Mathematics of Computation , vol. 47, pp. 323–346. 12.6ExternalStorageorMemory-LocalFFTs 525Sample 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).12.6 External Storage or Memory-Local FFTs Sometime in your life, you might have to compute the Fourier transform of a really largedata set, larger than the size of your computer ’s physical memory. In such a case, the data will be stored on some external medium, such as magnetic or optical tape or disk.Needed is an algorithm that makes some manageable number of sequential passes throughthe external data, processing it on the fly and outputting intermediate results to other external media, which can be read on subsequent passes. In fact, an algorithm of just this description was developed by Singleton [1]very soon after the discovery of the FFT. The algorithm requires four sequential storage devices, eachcapable of holding half of the input data. The first half of the input data is initially on one device, the second half on another. Singleton ’s algorithm is based on the observation that it is possible to bit-reverse 2 M values by the following sequence of operations: On the first pass, values are read alternately from the two input devices, and written to a single output device (until it holds half the data),and then to the other output device. On the second pass, the output devices become inputdevices, and vice versa. Now, we copy twovalues from the first device, then twovalues from the second, writing them (as before) first tofill one output device, then to fill a second. Subsequent passes read 4, 8, etc., input values at a time. After completion of pass M−1, the data are in bit-reverse order. Singleton ’s next observation is that it is possible to alternate the passes of essentially this bit-reversal technique with passes that implement one stage of the Danielson-Lanczoscombination formula (12.2.3). The scheme, roughly, is this: One starts as before with halfthe input data on one device, half on another. In the first pass, one complex value is read from each input device. Two combinations are formed, and one is written to each of twooutput devices. After this “computing ”pass, the devices are rewound, and a “permutation ” pass is performed, where groups of values are read from the first input device and alternately written to the first and second output devices; when the first input device is exhausted, the secondissimilarlyprocessed. ThissequenceofcomputingandpermutationpassesisrepeatedM−K−1times, where 2 Kis the size of internal buffer available to the program. The secondphaseofthecomputationconsistsofa final Kcomputationpasses. Whatdistinguishes the second phase from the first is that, now, the permutations are local enough to do in place during the computation. There are thus no separate permutation passes in the second phase.In all, there are 2M−K−2passes through the data. Here is an implementation of Singleton ’s algorithm, based on [1]: SUBROUTINE fourfs(iunit,nn,ndim,isign) INTEGER ndim,nn(ndim),isign,iunit(4),KBFPARAMETER (KBF=128) C USES fourew One- or multi-dimensional Fourier transform of a large data set stored on external media.On input, ndim is the number of dimensions, and nn(1:ndim) contains the lengths of each dimension (number of complex values), which must be powers of two. iunit(1:4) contains the unit numbers of 4 sequential files, each large enough to hold half of the data.The four units must be opened for FORTRAN unformatted access. The input data must be inFORTRAN normal order, with its first half stored on unit iunit(1) , its second half on iunit(2) , in unformatted form, with KBF real numbers per record. isign should be set to 1 for the Fourier transform, to −1for its inverse. On output, values in the array iunit may have been permuted; the first half of the result is stored on iunit(3) , the second half on iunit(4) . N.B.: For ndim >1, the output is stored by rows, i.e., notinFORTRAN normal order; in other words, the output is the transpose of that which would have been produced by routine fourn . INTEGER j,j12,jk,k,kk,n,mm,kc,kd,ks,kr,nr,ns,nv,jx, * mate(4),na,nb,nc,nd REAL tempr,tempi,afa(KBF),afb(KBF),afc(KBF) DOUBLE PRECISION wr,wi,wpr,wpi,wtemp,thetaSAVE mate DATA mate /2,1,4,3/