f12-1

494 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).PSD-per-unit-time converges to finite values at all frequencies exceptthose where h(t)has a discrete sine-wave (or cosine-wave) component of finite amplitude. At those frequencies, it becomes a delta-function, i.e., a sharp spike, whose width gets narrower and narrower, but whose area converges to be the mean square amplitude of the discrete sine or cosine component at that frequency. We have by now stated all of the analytical formalismthat we will need in this chapter with one exception: In computational work, especially with experimental data, we are almost never given a continuous function h(t)to work with, but are given,rather,a list of measurementsof h(ti)for a discrete set of ti’s. Theprofound implicationsofthis seeminglyunimportantfact are thesubject ofthe nextsection. CITED REFERENCES AND FURTHER READING: Champeney,D.C.1973, FourierTransformsandTheirPhysicalApplications (NewYork:Academic Press). Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork: Academic Press). 12.1 Fourier Transform of Discretely Sampled Data In the most common situations, function h(t)is sampled (i.e., its value is recorded)atevenlyspacedintervalsintime. Let ∆denotethetimeintervalbetween consecutive samples, so that the sequence of sampled values is hn=h(n∆) n=...,−3,−2,−1,0,1,2,3,... (12.1.1 ) The reciprocal of the time interval ∆is called the sampling rate ;i f∆is measured in seconds, for example, then the sampling rate is the number of samples recorded per second. SamplingTheorem andAliasing For any sampling interval ∆, there is also a special frequency fc, called the Nyquist critical frequency , given by fc≡1 2∆(12.1.2 ) Ifa sinewave oftheNyquistcriticalfrequencyis sampledat its positivepeakvalue, then the next sample will be at its negative trough value, the sample after that at the positive peak again, and so on. Expressed otherwise: Critical sampling of a sine wave is two sample points per cycle. One frequently chooses to measure time in units of the sampling interval ∆. In this case the Nyquist critical frequency is just the constant 1/2. TheNyquistcriticalfrequencyisimportantfortworelated,butdistinct,reasons. Oneis goodnews,andtheotherbadnews. Firstthe goodnews. Itis theremarkable 12.1FourierTransformofDiscretelySampledData 495Sample 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).fact known as the sampling theorem : If a continuous function h(t), sampled at an interval ∆,happenstobe bandwidthlimited tofrequenciessmallerinmagnitudethan fc,i.e.,if H(f)=0forall |f|≥fc,thenthefunction h(t)iscompletelydetermined by its samples hn. In fact, h(t)is given explicitly by the formula h(t)=∆+∞/summationdisplay n=−∞hnsin[2πfc(t−n∆)] π(t−n∆)(12.1.3 ) This is a remarkable theorem for many reasons, among them that it shows that the “information content” of a bandwidth limited function is, in some sense, infinitely smaller than that of a general continuous function. Fairly often, one is dealing with a signal that is known on physical grounds to be bandwidth limited (or at least approximately bandwidth limited). For example, the signal may have passed through an amplifier with a known, finite frequency response. In this case, thesampling theorem tells us that the entire information content of the signal can be recordedbysamplingitatarate ∆ −1equaltotwicethemaximumfrequencypassed by the amplifier (cf. 12.1.2). Nowthebadnews. Thebadnewsconcernstheeffectofsamplingacontinuous function that is notbandwidth limited to less than the Nyquist critical frequency. In that case, it turns out that all of the power spectral density that lies outside of the frequency range −fc<f<f cis spuriously moved into that range. This phenomenonis called aliasing. Any frequencycomponentoutside of the frequency range (−fc,fc)isaliased(falsely translated) into that range by the very act of discrete sampling. You can readily convince yourself that two waves exp(2 πif 1t) and exp(2 πif 2t)give the same samples at an interval ∆if and only if f1and f2differ by a multiple of 1/∆, which is just the width in frequency of the range (−fc,fc). There is little that you can do to remove aliased power once you have discretelysampleda signal. The wayto overcomealiasing is to (i) knowthe natural bandwidth limit of the signal — or else enforce a known limit by analog filtering of the continuous signal, and then (ii) sample at a rate sufficiently rapid to give atleast two points per cycle of the highest frequencypresent. Figure 12.1.1illustrates these considerations. To put the best face on this, we can take the alternative point of view: If a continuousfunctionhasbeencompetentlysampled,then,whenwecometoestimate its Fourier transform from the discrete samples, we can assume(or rather we might as wellassume) that its Fourier transform is equal to zero outside of the frequency rangein between −f candfc. Thenwe lookto theFouriertransformto tell whether thecontinuousfunction hasbeencompetentlysampled(aliasingeffectsminimized). We do this by looking to see whether the Fourier transform is already approaching zero as the frequency approaches fcfrom below, or −fcfrom above. If, on the contrary, the transform is going towards some finite value, then chances are thatcomponentsoutsideoftherangehavebeenfoldedbackoverontothecritical range. Discrete FourierTransform WenowestimatetheFouriertransformofafunctionfromafinitenumberofits sampled points. Suppose that we have Nconsecutive sampled values hk≡h(tk),t k≡k∆,k =0,1,2,...,N −1( 12.1.4 ) 496 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).h(t) t (a) f0H(f) (b) (c)aliased Fourier transform true Fourier transform 0H(f) 1 2∆1 2∆−f∆ T Figure 12.1.1. The continuous function shown in (a) is nonzero only for a finite interval of time T. It follows that its Fourier transform, whose modulus is shown schematically in (b), is not bandwidth limited but has finite amplitude for all frequencies. If the original function is sampled with a sampling interval ∆, as in (a), then the Fourier transform (c) is de fined only between plus and minus the Nyquist critical frequency. Power outside that range is folded over or “aliased”into the range. The effect can be eliminated only by low-pass filtering the original function before sampling . so that the sampling interval is ∆. To make things simpler, let us also suppose that Nis even. If the function h(t)is nonzero only in a finite interval of time, then that whole interval of time is supposedto be containedin the range of the Npoints given. Alternatively,ifthefunction h(t)goesonforever,thenthesampledpointsare supposed to be at least “typical”of what h(t)looks like at all other times. With Nnumbers of input, we will evidently be able to produce no more than Nindependent numbers of output. So, instead of trying to estimate the Fourier transform H(f)at all values of fin the range −fctofc, let us seek estimates only at the discrete values fn≡n N∆,n =−N 2,...,N 2(12.1.5 ) Theextremevaluesof nin(12.1.5)correspondexactlytothelowerandupperlimits of the Nyquist critical frequency range. If you are really on the ball, you will have noticed that there are N+1, not N, values of nin (12.1.5); it will turn out that the two extreme values of nare not independent (in fact they are equal), but all the others are. This reduces the count to N. 12.1FourierTransformofDiscretelySampledData 497Sample 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).Theremainingstepis to approximatetheintegralin (12.0.1)byadiscretesum: H(fn)=/integraldisplay∞ −∞h(t)e2πif ntdt≈N−1/summationdisplay k=0hke2πif ntk∆=∆N−1/summationdisplay k=0hke2πikn/N (12.1.6 ) Here equations (12.1.4) and (12.1.5) have been used in the final equality. The final summation in equation (12.1.6) is called the discrete Fourier transform of the N points hk. Let us denote it by Hn, Hn≡N−1/summationdisplay k=0hke2πikn/N(12.1.7 ) ThediscreteFouriertransformmaps Ncomplexnumbers(the hk’s)into Ncomplex numbers (the Hn’s). It does not depend on any dimensional parameter, such as the time scale ∆. The relation (12.1.6) between the discrete Fourier transform of a set ofnumbersandtheircontinuousFouriertransformwhentheyareviewedassamples of a continuous function sampled at an interval ∆can be rewritten as H(fn)≈∆Hn (12.1.8 ) where fnis given by (12.1.5). Uptonowwehavetakentheviewthattheindex nin(12.1.7)variesfrom −N/2 toN/2(cf. 12.1.5). Youcan easily see, however,that (12.1.7)is periodicin n, with period N. Therefore, H−n=HN−nn=1,2,.... With this conversionin mind, onegenerallylets the ninHnvaryfrom 0toN−1(onecompleteperiod). Then n andk(inhk) vary exactly over the same range, so the mapping of Nnumbers into Nnumbersis manifest. Whenthis conventionis followed,youmust rememberthat zero frequency corresponds to n=0, positive frequencies 0<f<f ccorrespond tovalues 1≤n≤N/2−1,whilenegativefrequencies −fc<f< 0correspondto N/2+1≤n≤N−1. Thevalue n=N/2correspondsto bothf=fcandf=−fc. ThediscreteFouriertransformhassymmetrypropertiesalmostexactlythesame as the continuous Fourier transform. For example, all the symmetries in the table following equation (12.0.3) hold if we read hkforh(t),HnforH(f), and HN−n forH(−f). (Likewise, “even”and“odd”intimerefertowhetherthevalues hkatk andN−kare identical or the negative of each other.) The formula for the discrete inverseFourier transform, which recovers the set ofhk’s exactly from the Hn’s is: hk=1 NN−1/summationdisplay n=0Hne−2πikn/N(12.1.9 ) Notice that the only differences between (12.1.9) and (12.1.7) are (i) changing the sign in the exponential, and (ii) dividing the answer by N. This means that a routineforcalculatingdiscreteFouriertransformscanalso,withslightmodi fication, calculate the inverse transforms. 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: De fineWas 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 Nlog2NandN2is 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