f12-3

504 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).integer arithmetic modulo some large prime N+1, and the Nth root of 1by the modulo arithmetic equivalent. Strictly speaking, these are not Fouriertransforms at all, but the properties are quite similar and computational speed can be far superior. On the other hand, their use is somewhat restricted to quantities like correlations and convolutions since the transform itself is not easily interpretableas a “frequency” spectrum. CITED REFERENCES AND FURTHER READING: Nussbaumer,H.J.1982, FastFourierTransformandConvolutionAlgorithms (NewYork:Springer- Verlag). Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork: Academic Press). Brigham, E.O. 1974, The Fast FourierTransform (Englewood Cliffs, NJ: Prentice-Hall). [1] Bloomfield, P. 1976, Fourier Analysis of Time Series – An Introduction (New York: Wiley). Van Loan, C. 1992, Computational Frameworks for the Fast Fourier Transform (Philadelphia: S.I.A.M.). Beauchamp, K.G. 1984, Applications of Walsh Functions and Related Functions (New York: Academic Press) [non-Fourier transforms]. Heideman, M.T., Johnson, D.H., and Burris, C.S. 1984, IEEE ASSP Magazine , pp. 14–21 (Oc- tober). 12.3 FFT of Real Functions, Sine and Cosine Transforms It happensfrequentlythat the data whose FFT is desired consist of real-valued samples fj,j =0...N −1. To use four1, we put these into a complex array with all imaginary parts set to zero. The resulting transform Fn,n =0...N −1 satisfies FN−n*= Fn. Since this complex-valued array has real values for F0 andFN/2, and (N/2)−1other independent values F1...F N/2−1, it has the same 2(N/2−1) + 2 = N“degrees of freedom”as the original, real data set. However, theuseofthefullcomplexFFTalgorithmforrealdataisinefficient,bothinexecution time and in storage required. You would think that there is a better way. There are twobetter ways. The first is “mass production”: Pack two separate real functionsinto the input array in such a way that their individualtransformscan be separated from the result. This is implemented in the program twofftbelow. This may remind you of a one-cent sale, at which you are coerced to purchase two of an item when you only need one. However, remember that for correlations and convolutions the Fourier transforms of two functions are involved, and this is a handy way to do them both at once. The second method is to pack the real input arraycleverly,withoutextra zeros, into a complexarray of half its length. One thenperforms a complex FFT on this shorter length; the trick is then to get the required answer out of the result. This is done in the program realftbelow. 12.3FFTofRealFunctions,SineandCosineTransforms 505Sample 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).Transformof TwoReal Functions Simultaneously First we show how to exploit the symmetry of the transform Fnto handle two real functions at once: Since the input data fjare real, the components of the discrete Fourier transform satisfy FN−n=(Fn)* ( 12.3.1 ) where the asterisk denotes complex conjugation. By the same token, the discrete Fouriertransformof a purelyimaginaryset of gj’s has the opposite symmetry. GN−n=−(Gn)* ( 12.3.2 ) Therefore we can take the discrete Fourier transform of two real functions each of length Nsimultaneously by packing the two data arrays as the real and imaginary parts,respectively,ofthecomplexinputarrayof four1. Thentheresultingtransform array can be unpackedinto two complex arrays with the aid of the two symmetries. Routine twofftworks out these ideas. SUBROUTINE twofft(data1,data2,fft1,fft2,n) INTEGER nREAL data1(n),data2(n) COMPLEX fft1(n),fft2(n) C USES four1 Given two real input arrays data1(1:n) anddata2(1:n) , this routine calls four1and returns two complex output arrays, fft1(1:n) andfft2(1:n) , each of complex length n (i.e., real length 2*n), which contain the discrete Fourier transforms of the respective data arrays. nMUST be an integer power of 2. INTEGER j,n2 COMPLEX h1,h2,c1,c2 c1=cmplx(0.5,0.0)c2=cmplx(0.0,-0.5)do 11j=1,n fft1(j)=cmplx(data1(j),data2(j)) Pack the two real arrays into one complex array. enddo 11 call four1(fft1,n,1) Transform the complex array. fft2(1)=cmplx(aimag(fft1(1)),0.0) fft1(1)=cmplx(real(fft1(1)),0.0)n2=n+2do 12j=2,n/2+1 h1=c1*(fft1(j)+conjg(fft1(n2-j))) U s es y m m e t r i e st os e p a r a t et h et w ot r a n s - forms. h2=c2*(fft1(j)-conjg(fft1(n2-j))) fft1(j)=h1 Ship them out in two complex arrays. fft1(n2-j)=conjg(h1) fft2(j)=h2fft2(n2-j)=conjg(h2) enddo 12 returnEND What about the reverse process? Suppose you have two complex transform arrays, each of which has the symmetry (12.3.1),so that you know that the inverses of both transforms are real functions. Can you invert both in a single FFT? This is even easier than the other direction. Use the fact that the FFT is linear and form the sum of the first transform plus itimes the second. Invert using four1with 506 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).isign =−1. The real and imaginary parts of the resulting complex array are the two desired real functions. FFTof Single Real Function To implement the second method, which allows us to perform the FFT of asinglereal function without redundancy, we split the data set in half, thereby forming two real arrays of half the size. We can apply the program above to thesetwo, but of course the result will not be the transform of the original data. It will be a schizophrenic combination of two transforms, each of which has half of the informationwe need. Fortunately,this schizophreniais treatable. It workslike this: The right way to split the original data is to take the even-numbered f jas one data set, and the odd-numbered fjas the other. The beauty of this is that we can take the original real array and treat it as a complex array hjof half the length. The first data set is the real part of this array, and the second is the imaginarypart, as prescribedfor twofft. No repackingis required. In otherwords hj=f2j+if2j+1,j =0,...,N / 2−1. We submit this to four1, and it will return a complex array Hn=Fe n+iFo n,n =0,...,N / 2−1with Fe n=N/2−1/summationdisplay k=0f2ke2πikn/ (N/2) Fo n=N/2−1/summationdisplay k=0f2k+1e2πikn/ (N/2)(12.3.3 ) Thediscussionofprogram twoffttellsyouhowtoseparatethetwotransforms Fe nandFo nout of Hn. How doyouworkthem intothe transform Fnofthe original data set fj? Simply glance back at equation (12.2.3): Fn=Fe n+e2πin/NFo n n=0,...,N −1( 12.3.4 ) Expressed directly in terms of the transform Hnof our real (masquerading as complex) data set, the result is Fn=1 2(Hn+HN/2−n*)−i 2(Hn−HN/2−n*)e2πin/Nn=0,...,N −1 (12.3.5 ) A few remarks: •Since FN−n*=Fnthere is no point in saving the entire spectrum. The positivefrequencyhalf is sufficientand canbe storedin thesame arrayas the original data. The operation can, in fact, be done in place. •Evenso,weneedvalues Hn,n =0,...,N / 2whereas four1returnsonly the values n=0,...,N / 2−1. Symmetryto the rescue, HN/2=H0. •Thevalues F0andFN/2arerealandindependent. Inordertoactuallyget the entire Fnin the original array space, it is convenient to return FN/2 as the imaginary part of F0. 12.3FFTofRealFunctions,SineandCosineTransforms 507Sample 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).•Despite its complicated form, the process above is invertible. First peel FN/2out of F0. Then construct Fe n=1 2(Fn+F* N/2−n) Fo n=1 2e−2πin/N(Fn−F* N/2−n)n=0,...,N / 2−1 (12.3.6 ) and use four1to find the inverse transform of Hn=F(1) n +iF(2) n. Surprisingly, the actual algebraic steps are virtually identical to those of the forward transform. Here is a representation of what we have said: SUBROUTINE realft(data,n,isign) INTEGER isign,nREAL data(n) C USES four1 Calculates the Fourier transform of a set of nreal-valued data points. Replaces this data (which is stored in array data(1:n) ) by the positive frequency half of its complex Fourier transform. The real-valued first and last components of the complex transform are returned as elements data(1) anddata(2) , respectively. nmust be a power of 2. This routine also calculates the inverse transform of a complex data array if it is the transform of realdata. (Result in this case must be multiplied by 2/n.) INTEGER i,i1,i2,i3,i4,n2p3 REAL c1,c2,h1i,h1r,h2i,h2r,wis,wrs DOUBLE PRECISION theta,wi,wpi,wpr, * wr,wtemp Doubleprecisionforthetrigonometric recurrences. theta=3.141592653589793d0/dble(n/2) Initialize the recurrence. c1=0.5if (isign.eq.1) then c2=-0.5 call four1(data,n/2,+1) The forward transform is here. else c2=0.5 Otherwise set up for an inverse transform. theta=-theta endifwpr=-2.0d0*sin(0.5d0*theta)**2wpi=sin(theta) wr=1.0d0+wpr wi=wpin2p3=n+3do 11i=2,n/4 Case i=1done separately below. i1=2*i-1 i2=i1+1i3=n2p3-i2 i4=i3+1 wrs=sngl(wr)wis=sngl(wi)h1r=c1*(data(i1)+data(i3)) The two separate transforms are separated out of data. h1i=c1*(data(i2)-data(i4)) h2r=-c2*(data(i2)+data(i4))h2i=c2*(data(i1)-data(i3)) data(i1)=h1r+wrs*h2r-wis*h2i Here they are recombined to form the true trans- form of the original real data. data(i2)=h1i+wrs*h2i+wis*h2r data(i3)=h1r-wrs*h2r+wis*h2idata(i4)=-h1i+wrs*h2i+wis*h2r wtemp=wr The recurrence. wr=wr*wpr-wi*wpi+wrwi=wi*wpr+wtemp*wpi+wi enddo 11 508 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).if (isign.eq.1) then h1r=data(1) data(1)=h1r+data(2) data(2)=h1r-data(2) Squeeze the first and last data together to get them all within the original array. else h1r=data(1) data(1)=c1*(h1r+data(2))data(2)=c1*(h1r-data(2))call four1(data,n/2,-1) Thisistheinversetransformforthecase isign=-1 . endif returnEND Fast Sine andCosine Transforms Amongtheirotheruses,theFouriertransformsoffunctionscanbeusedtosolve differential equations (see §19.4). The most common boundary conditions for the solutions are 1) they have the value zero at the boundaries, or 2) their derivatives are zero at the boundaries. In the first instance, the natural transform to use is thesinetransform, given by F k=N−1/summationdisplay j=1fjsin(πjk/N )sine transform (12.3.7 ) where fj,j =0,...,N −1is the data array, and f0≡0. AtfirstblushthisappearstobesimplytheimaginarypartofthediscreteFourier transform. However, the argument of the sine differs by a factor of two from thevalue that would make this so. The sine transformuses sines only as a completeset of functions in the interval from 0to2π, and, as we shall see, the cosine transform usescosinesonly . Bycontrast,thenormalFFTusesbothsinesandcosines,butonly half as many of each. (See Figure 12.3.1.) Theexpression(12.3.7)canbe“force-fit”intoaformthatallowsitscalculation viatheFFT. Theideais to extendthegivenfunctionrightwardpast its last tabulated value. We extend the data to twice their length in such a way as to make them an oddfunction about j=N, with f N=0, f2N−j≡−fj j=0,...,N −1( 12.3.8 ) Consider the FFT of this extended function: Fk=2N−1/summationdisplay j=0fje2πijk/ (2N)(12.3.9 ) The half of this sum from j=Ntoj=2N−1can be rewritten with the substitution j/prime=2N−j 2N−1/summationdisplay j=Nfje2πijk/ (2N)=N/summationdisplay j/prime=1f2N−j/primee2πi(2N−j/prime)k/(2N) =−N−1/summationdisplay j/prime=0fj/primee−2πij/primek/(2N)(12.3.10 ) 12.3FFTofRealFunctions,SineandCosineTransforms 509Sample 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).(a)+1 0 −1 +1 0 −1 +1 0 −1(b) (c) 0 2π54 21 3 1 23 4 5 1 2 3 4 5 Figure12.3.1. BasisfunctionsusedbytheFouriertransform(a),sinetransform(b),andcosinetransform (c), are plotted. The firstfive basis functions are shown in each case. (For the Fourier transform, the real and imaginary parts of the basis functions are both shown.) While some basis functions occur in more than one transform, the basis sets are distinct. For example, the sine transform functions labeled (1), (3), (5) are not present in the Fourier basis. Any of the three sets can expand any function in the intervalshown; however, the sine or cosine transform best expands functions matching the boundary conditionsof the respective basis functions, namely zero function values for sine, zero derivatives for cosine. so that Fk=N−1/summationdisplay j=0fj/bracketleftBig e2πijk/ (2N)−e−2πijk/ (2N)/bracketrightBig =2iN−1/summationdisplay j=0fjsin(πjk/N )(12.3.11 ) Thus,uptoafactor 2iwegetthesinetransformfromtheFFToftheextendedfunction. This method introduces a factor of two inef ficiency into the computation by extending the data. This inef ficiency shows up in the FFT output, which has zeros for the real part of every element of the transform. For a one-dimensionalproblem, the factor of two may be bearable, especially in view of the simplicity of the method. When we work with partial differential equations in two or three dimensions, though, the factor becomes four or eight, so efforts to eliminate the inefficiency are well rewarded. 510 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).From the originalreal data array fjwe will constructan auxiliaryarray yjand applyto it theroutine realft. Theoutputwill thenbeusedtoconstructthedesired transform. Forthesinetransformofdata fj,j =1,...,N −1,theauxiliaryarrayis y0=0 yj=s i n ( jπ/N )(fj+fN−j)+1 2(fj−fN−j) j=1,...,N −1(12.3.12 ) This array is of the same dimension as the original. Notice that the first term is symmetric about j=N/2and the second is antisymmetric. Consequently, when realftisappliedto yj,theresulthasrealparts Rkandimaginaryparts Ikgivenby Rk=N−1/summationdisplay j=0yjcos(2 πjk/N ) =N−1/summationdisplay j=1(fj+fN−j)s i n ( jπ/N )c o s ( 2 πjk/N ) =N−1/summationdisplay j=02fjsin(jπ/N )c o s ( 2 πjk/N ) =N−1/summationdisplay j=0fj/bracketleftbigg sin(2k+1 )jπ N−sin(2k−1)jπ N/bracketrightbigg =F2k+1−F2k−1 (12.3.13 ) Ik=N−1/summationdisplay j=0yjsin(2 πjk/N ) =N−1/summationdisplay j=1(fj−fN−j)1 2sin(2 πjk/N ) =N−1/summationdisplay j=0fjsin(2 πjk/N ) =F2k (12.3.14 ) Therefore Fkcan be determined as follows: F2k=Ik F2k+1=F2k−1+Rk k=0,..., (N/2−1) (12.3.15 ) The even terms of Fkare thus determined very directly. The odd terms require a recursion, the starting point of which follows from setting k=0in equation (12.3.15) and using F1=−F−1: F1=1 2R0 (12.3.16 ) The implementing program is 12.3FFTofRealFunctions,SineandCosineTransforms 511Sample 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 sinft(y,n) INTEGER n REAL y(n) C USES realft Calculates the sine transform of a set of nreal-valued data points stored in array y(1:n). The number nmust be a power of 2. On exit yis replaced by its transform. This program, without changes, also calculates the inverse sine transform, but in this case the output arrayshould be multiplied by 2/n. INTEGER j REAL sum,y1,y2 DOUBLE PRECISION theta,wi,wpi,wpr, * wr,wtemp Double precision in the trigonometric recurrences. theta=3.141592653589793d0/dble(n) Initialize the recurrence. wr=1.0d0 wi=0.0d0wpr=-2.0d0*sin(0.5d0*theta)**2 wpi=sin(theta) y(1)=0.0do 11j=1,n/2 wtemp=wr wr=wr*wpr-wi*wpi+wr Calculate the sine for the auxiliary array. wi=wi*wpr+wtemp*wpi+wi The cosine is needed to continue the recurrence. y1=wi*(y(j+1)+y(n-j+1)) Construct the auxiliary array. y2=0.5*(y(j+1)-y(n-j+1)) y(j+1)=y1+y2 Terms jand N−jare related. y(n-j+1)=y1-y2 enddo 11 call realft(y,n,+1) Transform the auxiliary array. sum=0.0y(1)=0.5*y(1) Initialize the sum used for odd terms below. y(2)=0.0 do 12j=1,n-1,2 sum=sum+y(j)y(j)=y(j+1) Even terms in the transform are determined directly. y(j+1)=sum Odd terms are determined by this running sum. enddo 12 returnEND The sine transform, curiously, is its own inverse. If you apply it twice, you get the original data, but multiplied by a factor of N/2. The other common boundary condition for differential equations is that the derivative of the function is zero at the boundary. In this case the natural transform is thecosinetransform. There are several possible ways of de fining the transform. Each can be thought of as resulting from a different way of extendinga given array tocreateanevenarrayofdoublethelength,and/orfromwhethertheextendedarray contains 2N−1,2N, or some other number of points. In practice, only two of the numerouspossibilities are useful so we will restrict ourselvesto just these two. Thefirst form of the cosine transform uses N+1data points: Fk=1 2[f0+(−1)kfN]+N−1/summationdisplay j=1fjcos(πjk/N )( 12.3.17 ) It results from extendingthe given array to an even array about j=N, with f2N−j=fj,j =0,...,N −1( 12.3.18 ) 512 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).Ifyousubstitutethisextendedarrayintoequation(12.3.9),andfollowstepsanalogous to those leading up to equation(12.3.11),you will find that the Fourier transformis justtwicethecosinetransform(12.3.17). Anotherwayofthinkingabouttheformula (12.3.17)istonoticethatitistheChebyshevGauss-Lobattoquadratureformula(see §4.5),oftenusedinClenshaw-Curtisadaptivequadrature(see §5.9,equation5.9.4). Onceagainthetransformcanbecomputedwithoutthefactoroftwoinef ficiency. In this case the auxiliary function is yj=1 2(fj+fN−j)−sin(jπ/N )(fj−fN−j) j=0,...,N −1(12.3.19 ) Instead of equation (12.3.15), realftnow gives F2k=Rk F2k+1=F2k−1+Ik k=0,..., (N/2−1) (12.3.20 ) The starting value for the recursion for odd kin this case is F1=1 2(f0−fN)+N−1/summationdisplay j=1fjcos(jπ/N )( 12.3.21 ) This sum does not appear naturally among the RkandIk, and so we accumulate it during the generation of the array yj. Once again this transform is its own inverse, and so the following routine works for both directions of the transformation. Note that although this form of the cosine transform has N+1input and output values, it passes an array only of length Ntorealft. SUBROUTINE cosft1(y,n) INTEGER n REAL y(n+1) C USES realft Calculates the cosine transform of a set y(1:n+1) of real-valued data points. The trans- formed data replace the original data in array y.nmust be a power of 2. This program, without changes, also calculates the inverse cosine transform, but in this case the output array should be multiplied by 2/n. INTEGER jREAL sum,y1,y2 DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp For trig. recurrences. theta=3.141592653589793d0/n Initialize the recurrence. wr=1.0d0 wi=0.0d0 wpr=-2.0d0*sin(0.5d0*theta)**2wpi=sin(theta)sum=0.5*(y(1)-y(n+1)) y(1)=0.5*(y(1)+y(n+1)) do 11j=1,n/2-1 j=n/2 unnecessary since y(n/2+1) unchanged. wtemp=wr wr=wr*wpr-wi*wpi+wr Carry out the recurrence. wi=wi*wpr+wtemp*wpi+wiy1=0.5*(y(j+1)+y(n-j+1)) Calculate the auxiliary function. y2=(y(j+1)-y(n-j+1)) y(j+1)=y1-wi*y2 The values for jand N−jare related. y(n-j+1)=y1+wi*y2sum=sum+wr*y2 Carry along this sum for later use in unfolding the transform. enddo 11 12.3FFTofRealFunctions,SineandCosineTransforms 513Sample 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).call realft(y,n,+1) Calculate the transform of the auxiliary function. y(n+1)=y(2) y(2)=sum sum is the value of F1in equation (12.3.21). do12j=4,n,2 sum=sum+y(j) Equation (12.3.20). y(j)=sum enddo 12 returnEND The second important form of the cosine transform is de fined by Fk=N−1/summationdisplay j=0fjcosπk(j+1 2) N(12.3.22 ) with inverse fj=2 NN−1/summationdisplay/prime k=0Fkcosπk(j+1 2) N(12.3.23 ) Here the prime on the summation symbol means that the term for k=0has a coefficient of1 2in front. This form arises by extending the given data, de fined for j=0,...,N −1,toj=N,..., 2N−1insuchawaythatitisevenaboutthepoint N−1 2and periodic. (It is therefore also even about j=−1 2.) The form (12.3.23) is related to Gauss-Chebyshev quadrature (see equation 4.5.19), to Chebyshev approximation( §5.8, equation 5.8.7),and Clenshaw-Curtis quadrature( §5.9). This formof the cosine transformis usefulwhen solving differentialequations on“staggered ”grids,wherethevariablesarecenteredmidwaybetweenmeshpoints. It is also the standardformin the field ofdata compressionandimage processing. The auxiliary functionused in this case is similar to equation (12.3.19): yj=1 2(fj+fN−j−1)+s i nπ(j+1 2) N(fj−fN−j−1) j=0,...,N −1 (12.3.24 ) Carryingoutthestepssimilartothoseusedtogetfrom(12.3.12)to(12.3.15),we find F2k=c o sπk NRk−sinπk NIk (12.3.25 ) F2k−1=s i nπk NRk+c o sπk NIk+F2k+1 (12.3.26 ) Note that equation (12.3.26) gives FN−1=1 2RN/2 (12.3.27 ) Thus the even components are found directly from (12.3.25), while the odd com- ponents are found by recursing (12.3.26)down from k=N/2−1, using (12.3.27) to start. Since the transform is not self-inverting, we have to reverse the above steps to find the inverse. Here is the routine: 514 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 cosft2(y,n,isign) INTEGER isign,n REAL y(n) C USES realft Calculates the “staggered” cosine transform of a set y(1:n)of real-valued data points. The transformed data replace the original data in array y.nmust be a power of 2. Set isignto+1for a transform, and to −1for an inverse transform. For an inverse transform, the output array should be multiplied by 2/n. INTEGER i REAL sum,sum1,y1,y2,ytemp DOUBLE PRECISION theta,wi,wi1,wpi,wpr,wr,wr1,wtemp,PI Double precision for the trigonometric recurrences. PARAMETER (PI=3.141592653589793d0) theta=0.5d0*PI/n Initialize the recurrences. wr=1.0d0wi=0.0d0 wr1=cos(theta) wi1=sin(theta)wpr=-2.0d0*wi1**2wpi=sin(2.d0*theta) if(isign.eq.1)then Forward transform. do 11i=1,n/2 y1=0.5*(y(i)+y(n-i+1)) Calculate the auxiliary function. y2=wi1*(y(i)-y(n-i+1)) y(i)=y1+y2y(n-i+1)=y1-y2wtemp=wr1 Carry out the recurrence. wr1=wr1*wpr-wi1*wpi+wr1 wi1=wi1*wpr+wtemp*wpi+wi1 enddo 11 call realft(y,n,1) Calculate the transform of the auxiliary function. do12i=3,n,2 Even terms. wtemp=wrwr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi y1=y(i)*wr-y(i+1)*wiy2=y(i+1)*wr+y(i)*wiy(i)=y1 y(i+1)=y2 enddo 12 sum=0.5*y(2) Initialize recurrence for odd terms with1 2RN/ 2. do13i=n,2,-2 Carry out recurrence for odd terms. sum1=sumsum=sum+y(i)y(i)=sum1 enddo 13 else if(isign.eq.-1)then Inverse transform. ytemp=y(n) do14i=n,4,-2 Form difference of odd terms. y(i)=y(i-2)-y(i) enddo 14 y(2)=2.0*ytemp do15i=3,n,2 Calculate Rkand Ik. wtemp=wrwr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi y1=y(i)*wr+y(i+1)*wiy2=y(i+1)*wr-y(i)*wiy(i)=y1 y(i+1)=y2 enddo 15 call realft(y,n,-1)do 16i=1,n/2 Invert auxiliary array. y1=y(i)+y(n-i+1) 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 thefirstN/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 fixedNof 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 de fine its two-dimensional discrete Fouriertransformas a complexfunction H(n1,n2),d efinedoverthe 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,