f12-2

498 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 discrete form of Parseval’s theorem is N−1/summationdisplay k=0|hk|2=1 NN−1/summationdisplay n=0|Hn|2(12.1.10 ) Therearealsodiscreteanalogstotheconvolutionandcorrelationtheorems(equations 12.0.9and 12.0.11),but we shall defer them to §13.1and §13.2, respectively. CITED REFERENCES AND FURTHER READING: Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall). Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork: Academic Press). 12.2 Fast Fourier Transform (FFT) HowmuchcomputationisinvolvedincomputingthediscreteFouriertransform (12.1.7) of Npoints? For many years, until the mid-1960s, the standard answer was this: Define Was the complex number W≡e2πi/N(12.2.1 ) Then (12.1.7) can be written as Hn=N−1/summationdisplay k=0Wnkhk (12.2.2 ) In other words, the vector of hk’s is multiplied by a matrix whose (n, k )th element is the constant Wto the power n×k. The matrix multiplication produces a vector resultwhosecomponentsarethe Hn’s. Thismatrixmultiplicationevidentlyrequires N2complex multiplications, plus a smaller number of operations to generate the required powers of W. So, the discrete Fourier transform appears to be an O(N2) process. These appearances are deceiving! The discrete Fourier transform can,in fact, be computed in O(Nlog 2N)operations with an algorithm called the fast Fourier transform ,o rFFT. The difference between Nlog2Nand N2is immense. With N=1 06,forexample,itisthedifferencebetween,roughly,30secondsofCPU timeand2weeksofCPUtimeonamicrosecondcycletimecomputer. Theexistence ofanFFTalgorithmbecamegenerallyknownonlyinthemid-1960s,fromtheworkofJ.W.CooleyandJ.W.Tukey. Retrospectively,wenowknow(see [1])thatefficient methods for computing the DFT had been independently discovered, and in some cases implemented,byas manyas adozenindividuals,startingwithGauss in1805! One“rediscovery”oftheFFT,thatofDanielsonandLanczosin1942,provides one of the clearest derivations of the algorithm. Danielson and Lanczos showed that a discrete Fourier transform of length Ncan be rewritten as the sum of two discrete Fouriertransforms, each of length N/2. One of the two is formedfromthe 12.2 FastFourierTransform(FFT) 499Sample 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).even-numbered points of the original N, the other from the odd-numbered points. The proof is simply this: Fk=N−1/summationdisplay j=0e2πijk/Nfj =N/2−1/summationdisplay j=0e2πik (2j)/Nf2j+N/2−1/summationdisplay j=0e2πik (2j+1)/Nf2j+1 =N/2−1/summationdisplay j=0e2πikj/ (N/2)f2j+WkN/2−1/summationdisplay j=0e2πikj/ (N/2)f2j+1 =Fe k+WkFo k(12.2.3 ) In the last line, Wis the same complex constant as in (12.2.1), Fe kdenotes the kth componentoftheFouriertransformoflength N/2formedfromtheevencomponents of the original fj’s, while Fo kis the correspondingtransform of length N/2formed from the odd components. Notice also that kin the last line of (12.2.3)varies from 0toN, not just to N/2. Nevertheless, the transforms Fe kand Fo kare periodic in k with length N/2. So each is repeated through two cycles to obtain Fk. Thewonderfulthingaboutthe Danielson-LanczosLemma isthatitcanbeused recursively. Having reduced the problem of computing Fkto that of computing Fe kand Fo k, we can do the same reduction of Fe kto the problem of computing the transform of itsN/4even-numbered input data and N/4odd-numbered data. In other words, we can define Fee kand Feo kto be the discrete Fourier transforms of the points which are respectively even-even and even-odd on the successive subdivisions of the data. Although there are ways of treating other cases, by far the easiest case is the one in which the original Nis an integer power of 2. In fact, we categorically recommendthatyou onlyuseFFTswith Napoweroftwo. Ifthelengthofyourdata setisnotapoweroftwo,paditwithzerosuptothenextpoweroftwo. (Wewillgive more sophisticated suggestions in subsequent sections below.) With this restriction onN, it is evident that we can continue applying the Danielson-Lanczos Lemma until we have subdividedthe data all the way downto transformsof length 1. What is the Fouriertransformof lengthone? It is just theidentityoperationthatcopiesitsoneinputnumberintoitsoneoutputslot! Inotherwords,foreverypatternof log 2N e’s and o’s, thereis a one-pointtransformthat is just oneof the inputnumbers fn Feoeeoeo ···oee k =fnforsome n (12.2.4 ) (Ofcoursethisone-pointtransformactuallydoesnotdependon k,sinceitisperiodic inkwith period 1.) Thenexttrickis tofigureoutwhichvalueof ncorrespondstowhichpatternof e’s and o’s in equation (12.2.4). The answer is: Reverse the pattern of e’s and o’s, then let e=0and o=1, and you will have, in binary the value of n. Do you see whyitworks? Itisbecausethesuccessivesubdivisionsofthedataintoevenandodd aretestsofsuccessivelow-order(leastsignificant)bitsof n. Thisideaof bitreversal can be exploited in a very clever way which, along with the Danielson-Lanczos 500 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).000 001010011100101110111000 001010011100101110111000 001010011100 101 110111 (a) (b) Figure 12.2.1. Reordering an array (here of length 8) by bit reversal, (a) between two arrays, versus (b) in place. Bit reversal reordering is a necessary part of the fast Fourier transform (FFT) algorithm. Lemma, makes FFTs practical: Suppose we take the original vector of data fj and rearrange it into bit-reversed order (see Figure 12.2.1), so that the individual numbers are in the order not of j, but of the number obtained by bit-reversing j. ThenthebookkeepingontherecursiveapplicationoftheDanielson-LanczosLemma becomes extraordinarily simple. The points as given are the one-point transforms. We combineadjacentpairstogettwo-pointtransforms,thencombineadjacentpairsof pairs to get 4-point transforms, and so on, until the first and second halves of the whole data set are combined into the final transform. Each combination takes of order Noperations, and there are evidently log 2Ncombinations, so the whole algorithmis of order Nlog2N(assuming, as is the case, that the process of sorting into bit-reversed order is no greater in order than Nlog2N). This, then, is the structure of an FFT algorithm: It has two sections. The first sectionsortsthedataintobit-reversedorder. Luckilythistakesnoadditionalstorage, sinceitinvolvesonlyswappingpairsofelements. (If k1isthebitreverseof k2,then k2is the bit reverse of k1.) The second section has an outer loop that is executed log2Ntimes and calculates, in turn, transforms of length 2,4,8,...,N. For each stage of this process, two nested inner loops range over the subtransforms alreadycomputedandtheelementsofeachtransform,implementingtheDanielson-Lanczos Lemma. The operation is made more ef ficient by restricting external calls for trigonometricsines and cosines to the outer loop, where they are made only log 2N times. Computation of the sines and cosines of multiple angles is through simple recurrence relations in the inner loops (cf. 5.5.6). The FFT routine given below is based on one originally written by N. M. Brenner. The input quantities are the number of complex data points ( nn), the data array ( data), and isign, which should be set to either ±1and is the sign of iin the exponential of equation (12.1.7). When isignis set to −1, the routine thus calculates the inverse transform (12.1.9) —except that it does not multiply by the normalizingfactor 1/Nthat appearsin that equation. You can do that yourself. Notice that the argument nnis the number of complexdata points, although 12.2 FastFourierTransform(FFT) 501Sample 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 avoid the use of complex arithmetic because of the inef ficient implementations found on many computers. The actual length of the real array ( data) is 2 times nn, with each complex value occupying two consecutive locations. In other words, data(1) is the real part of f0,data(2) is the imaginary part of f0, and so on up todata(2*nn-1) , which is the real part of fN−1, and data(2*nn) , which is the imaginary part of fN−1. The FFT routine returns the Fn’s packed in exactly the same fashion, as nncomplex numbers. Therealandimaginarypartsofthezerofrequencycomponent F0arein data(1) anddata(2);thesmallestnonzeropositivefrequencyhasrealandimaginarypartsin data(3) anddata(4);the smallest(inmagnitude)nonzeronegativefrequencyhas real and imaginary parts in data(2*nn-1) anddata(2*nn) . Positive frequencies increasing in magnitude are stored in the real-imaginary pairs data(5), data(6) up to data(nn-1), data(nn) . Negative frequencies of increasing magnitude are stored in data(2*nn-3), data(2*nn-2) down to data(nn+3), data(nn+4) . Finally,thepair data(nn+1), data(nn+2) containtherealandimaginarypartsof theonealiasedpointthatcontainsthemostpositiveandthemostnegativefrequency.You should try to develop a familiarity with this storage arrangement of complex spectra, also shown in Figure 12.2.2, since it is the practical standard. SUBROUTINE four1(data,nn,isign) INTEGER isign,nn REAL data(2*nn) Replaces data(1:2*nn) by its discrete Fourier transform, if isign is input as 1; or replaces data(1:2*nn) bynntimes its inverse discrete Fourier transform, if isign is input as −1. data is a complex array of length nn or, equivalently, a real array of length 2*nn .nn MUST be an integer power of 2 (this is not checked for!). INTEGER i,istep,j,m,mmax,nREAL tempi,tempr DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp Double precision for the trigonomet- ric recurrences. n=2*nn j=1 do 11i=1,n,2 This is the bit-reversal section of the routine. if(j.gt.i)then tempr=data(j) Exchange the two complex numbers. tempi=data(j+1) data(j)=data(i) data(j+1)=data(i+1)data(i)=temprdata(i+1)=tempi endif m=nn 1 if ((m.ge.2).and.(j.gt.m)) then j=j-m m=m/2 goto 1endif j=j+m enddo 11 mmax=2 Here begins the Danielson-Lanczos section of the routine. 2 if (n.gt.mmax) then Outer loop executed log2nntimes. istep=2*mmaxtheta=6.28318530717959d0/(isign*mmax) Initialize for the trigonometric recur- rence. wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0wi=0.d0 do 13m=1,mmax,2 Here are the two nested inner loops. 502 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).1 234real imagrealimagt = 0 t = ∆ real imagrealimagt = (N − 2)∆ t = (N − 1)∆real array of length 2 N1 234real imagrealimagf = 0 f = N − 1 NN + 1 N + 2 N + 3 N + 4real imagrealimagrealimagf = f = ± (combination) f = −real array of length 2 N1 N∆ N/2 − 1 N∆ 1 2∆ 2N − 1 2Nreal imagf = − 1 N∆2N − 3 2N − 2 2N − 1 2NN/2 − 1 N∆ (b) (a) Figure 12.2.2. Input and output arrays for FFT. (a) The input array contains N(a power of 2) complex time samples in a real array of length 2N, with real and imaginary parts alternating. (b) The output array contains the complex Fourier spectrum at Nvalues of frequency. Real and imaginary parts again alternate. The array starts with zero frequency, works up to the most positive frequency (which is ambiguous with the most negative frequency). Negative frequencies follow, from the second-mostnegative up to the frequency just below zero. do 12i=m,n,istep j=i+mmax This is the Danielson-Lanczos formula: tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j) data(j)=data(i)-tempr data(j+1)=data(i+1)-tempidata(i)=data(i)+temprdata(i+1)=data(i+1)+tempi enddo 12 wtemp=wr Trigonometric recurrence. wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi enddo 13 mmax=istep goto 2 Not yet done. endif All done. returnEND (Adoubleprecisionversionof four1,named dfour1, is usedbythe routine mpmul in§20.6. You can easily make the conversion, or else get the converted routine from the Numerical Recipes diskette.) 12.2 FastFourierTransform(FFT) 503Sample 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).OtherFFT Algorithms WeshouldmentionthatthereareanumberofvariantsonthebasicFFTalgorithm given above. As we have seen, that algorithm first rearranges the input elements into bit-reverse order, then builds up the output transform in log2Niterations. In the literature, this sequence is called a decimation-in-time orCooley-Tukey FFT algorithm. It is also possible to derive FFT algorithms that first go through a set of log2Niterations on the input data, and rearrange the outputvalues into bit-reverse order. Thesearecalled decimation-in-frequency orSande-Tukey FFTalgorithms. For some applications,such as convolution( §13.1),one takes a data set into the Fourier domainandthen,aftersomemanipulation,backoutagain. Inthesecasesitispossible to avoid all bit reversing. You use a decimation-in-frequencyalgorithm (without itsbit reversing)to get into the “scrambled ”Fourier domain, do your operations there, and then use an inverse algorithm (without itsbit reversing) to get back to the time domain. While elegant in principle, this procedure does not in practice save much computationtime, since the bit reversalsrepresentonlya small fractionof an FFT ’s operations count, and since most useful operations in the frequencydomain requirea knowledge of which points correspond to which frequencies. Another class of FFTs subdivides the initial data set of length Nnot all the way down to the trivial transform of length 1, but rather only down to some other small powerof 2, forexample N=4,base-4FFTs ,o r N=8,base-8FFTs . These small transforms are then done by small sections of highly optimized coding which take advantage of special symmetries of that particular small N. For example, for N=4, the trigonometric sines and cosines that enter are all ±1or0, so many multiplications are eliminated, leaving largely additions and subtractions. Thesecan be faster than simpler FFTs by some signi ficant, but not overwhelming, factor, e.g., 20 or 30 percent. TherearealsoFFTalgorithmsfordatasetsoflength Nnotapoweroftwo. They work by using relations analogous to the Danielson-Lanczos Lemma to subdivide the initial problem into successively smaller problems, not by factors of 2, but by whatever small prime factors happen to divide N. The larger that the largest prime factor of Nis, the worse this method works. If Nis prime, then no subdivision is possible, and the user (whether he knows it or not) is taking a slowFourier transform, of order N 2instead of order Nlog2N. Our advice is to stay clear of such FFT implementations, with perhaps one class of exceptions, the Winograd Fourier transform algorithms . Winogradalgorithms are in some ways analogousto the base-4 and base-8 FFTs. Winograd has derived highly optimized codings for taking small- Ndiscrete Fourier transforms, e.g., for N=2 ,3,4,5,7,8,11,13,16. The algorithms also use a new and clever way of combining the subfactors. The methodinvolvesa reorderingofthedatabothbeforethehierarchicalprocessingand after it, but it allows a signi ficant reduction in the number of multiplications in the algorithm. Forsomeespeciallyfavorablevaluesof N,the Winogradalgorithmscan be significantly (e.g., up to a factor of 2) faster than the simpler FFT algorithms of the nearest integer power of 2. This advantage in speed, however, must beweighed against the considerablymore complicated data indexinginvolvedin these transforms,and the fact that the Winogradtransformcannot be done “in place.” Finally, an interesting class of transforms for doing convolutions quickly are number theoretic transforms. These schemes replace floating-point arithmetic with 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 satisfiesFN−n*= Fn. Since this complex-valued array has real values for F0 and FN/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, theuseofthefullcomplexFFTalgorithmforrealdataisinef ficient,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.