f12-0

Sample 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).Chapter 12. Fast Fourier Transform 12.0 Introduction A verylargeclass of importantcomputationalproblemsfalls underthe general rubric of “Fourier transform methods” or “spectral methods.” For some of these problems, the Fourier transform is simply an efficient computational tool for accomplishing certain common manipulations of data. In other cases, we haveproblemsfor whichthe Fouriertransform(or the related“powerspectrum”)is itself of intrinsic interest. These two kindsof problemssharea commonmethodology. LargelyforhistoricalreasonstheliteratureonFourierandspectralmethodshas been disjoint from the literature on “classical” numerical analysis. Nowadays there isnojustificationforsuchasplit. Fouriermethodsarecommonplaceinresearchandwe shall not treat them as specialized or arcane. At the same time, we realize that manycomputerusershavehadrelativelylessexperiencewiththisfieldthanwith,say, differentialequationsornumericalintegration. Thereforeoursummaryofanalyticalresultswillbemorecomplete. Numericalalgorithms,perse,beginin §12.2. Various applications of Fourier transform methods are discussed in Chapter 13. A physicalprocesscan bedescribedeither inthe time domain ,bythevaluesof some quantity has a function of time t, e.g., h(t), or else in the frequency domain , where the process is specified by giving its amplitude H(generally a complex number indicating phase also) as a function of frequency f, that is H(f), with −∞ <f< ∞. For many purposes it is useful to think of h(t)andH(f)as being twodifferent representations ofthesamefunction. Onegoesbackandforthbetween these two representations by means of the Fourier transform equations, H(f)=/integraldisplay ∞ −∞h(t)e2πiftdt h(t)=/integraldisplay∞ −∞H(f)e−2πiftdf(12.0.1 ) Iftis measuredin seconds, then fin equation(12.0.1)is in cycles per second, orHertz(theunitoffrequency). However,theequationsworkwithotherunitstoo. Ifhis a functionof position x(in meters), Hwill be a functionof inversewavelength (cyclespermeter),andsoon. Ifyouaretrainedasaphysicistormathematician,you are probablymore usedto using angularfrequency ω, whichis givenin radiansper sec. The relation between ωandf,H(ω)andH(f)is ω≡2πf H (ω)≡[H(f)] f=ω/2π(12.0.2 ) 490 12.0 Introduction 491Sample 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).and equation (12.0.1) looks like this H(ω)=/integraldisplay∞ −∞h(t)eiωtdt h(t)=1 2π/integraldisplay∞ −∞H(ω)e−iωtdω(12.0.3 ) We were raised on the ω-convention, but we changed! There are fewer factors of 2πto remember if you use the f-convention, especially when we get to discretely sampled data in §12.1. From equation (12.0.1) it is evident at once that Fourier transformation is a linearoperation. The transform of the sum of two functions is equal to the sum of the transforms. The transform of a constant times a function is that same constant times the transform of the function. In the time domain, function h(t)may happen to have one or more special symmetries It might be purely real orpurely imaginary or it might be even, h(t)=h(−t),o rodd,h(t)=−h(−t). In the frequencydomain, these symmetries lead to relationships between H(f)andH(−f). The following table gives the correspondence between symmetries in the two domains: If... then... h(t)is real H(−f)=[H(f)]* h(t)is imaginary H(−f)=−[H(f)]* h(t)is even H(−f)=H(f)[i.e., H(f)is even] h(t)is odd H(−f)=−H(f)[i.e., H(f)is odd] h(t)is realandeven H(f)is real andeven h(t)is realandodd H(f)is imaginaryandodd h(t)is imaginaryandeven H(f)is imaginaryandeven h(t)is imaginaryandodd H(f)is real andodd In subsequent sections we shall see how to use these symmetries to increase computational efficiency. Herearesomeotherelementarypropertiesofthe Fouriertransform. (We’ll use the “⇐⇒” symbol to indicate transform pairs.) If h(t)⇐⇒ H(f) is such a pair, then other transform pairs are h(at)⇐⇒1 |a|H(f a) “timescaling” (12.0.4 ) 1 |b|h(t b)⇐⇒ H(bf) “frequencyscaling” (12.0.5 ) h(t−t0)⇐⇒ H(f)e2πift 0“timeshifting” (12.0.6 ) h(t)e−2πif 0t⇐⇒ H(f−f0)“frequencyshifting” (12.0.7 ) 492 Chapter12. Fast FourierTransformSample 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).With two functions h(t)andg(t), and their corresponding Fourier transforms H(f)andG(f), wecanformtwocombinationsofspecialinterest. The convolution of the two functions, denoted g∗h, is defined by g∗h≡/integraldisplay∞ −∞g(τ)h(t−τ)dτ (12.0.8 ) Note that g∗his a function in the time domain and that g∗h=h∗g. It turns out that the function g∗his one member of a simple transform pair g∗h⇐⇒ G(f)H(f)“ConvolutionTheorem” (12.0.9 ) In other words, the Fourier transform of the convolution is just the product of the individual Fourier transforms. Thecorrelation of two functions, denoted Corr (g,h), is defined by Corr(g,h)≡/integraldisplay∞ −∞g(τ+t)h(τ)dτ (12.0.10 ) Thecorrelationisafunctionof t,whichiscalledthe lag. Itthereforeliesinthetime domain, and it turns out to be one member of the transform pair: Corr(g,h)⇐⇒ G(f)H*(f)“CorrelationTheorem” (12.0.11 ) [Moregenerally,thesecondmemberofthepairis G(f)H(−f),butwearerestricting ourselvestotheusualcaseinwhich gandharerealfunctions,sowetakethelibertyof setting H(−f)=H*(f).] This result shows that multiplyingthe Fouriertransform ofonefunctionbythecomplexconjugateoftheFouriertransformoftheothergivestheFouriertransformoftheircorrelation. Thecorrelationofafunctionwithitself is called its autocorrelation . In this case (12.0.11)becomes the transformpair Corr(g,g)⇐⇒ | G(f)| 2“Wiener-KhinchinTheorem” (12.0.12 ) Thetotal power in a signal is the same whether we compute it in the time domainor in the frequencydomain. This result is knownas Parseval’s theorem : TotalPower ≡/integraldisplay∞ −∞|h(t)|2dt=/integraldisplay∞ −∞|H(f)|2df (12.0.13 ) Frequentlyonewantstoknow“howmuchpower”iscontainedinthefrequency interval between fandf+df. In such circumstances one does not usually distinguish between positive and negative f, but rather regards fas varying from 0 (“zero frequency”or D.C.) to +∞. In such cases, one defines the one-sided power spectral density (PSD) of the function has Ph(f)≡|H(f)|2+|H(−f)|20≤f<∞ (12.0.14 ) sothatthetotalpoweris justtheintegralof Ph(f)from f=0tof=∞. Whenthe function h(t)isreal,thenthetwotermsin(12.0.14)areequal,so Ph(f)=2|H(f)|2. 12.0 Introduction 493Sample 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)2 (a) (b) (c)fPh(f) (one-sided) 0 −f Ph(f) (two-sided)t f 0 Figure 12.0.1. Normalizations of one- and two-sided power spectra. The area under the square of the function, (a), equals the area under its one-sided power spectrum at positive frequencies, (b), and also equals the area under its two-sided power spectrum at positive and negative frequencies, (c). Be warned that one occasionally sees PSDs de fined without this factor two. These, strictly speaking, are called two-sided power spectral densities , but some books are not careful about stating whether one- or two-sided is to be assumed. We will always use the one-sided density given by equation (12.0.14). Figure 12.0.1 contrasts the two conventions. If the function h(t)goes endlessly from −∞ <t< ∞, then its total power and power spectral density will, in general, be in finite. Of interest then is the (one- or two-sided) power spectral density per unit time . This is computed by taking a long, but finite, stretch of the function h(t), computing its PSD [that is, the PSD of a function that equals h(t)in thefinite stretch but is zero everywhere else], and thendividingtheresultingPSDbythelengthofthestretchused. Parseval ’stheorem in this case states that the integral of the one-sided PSD-per-unit-time over positive frequency is equal to the mean square amplitude of the signal h(t). Youmightwell worryabouthowthePSD-per-unit-time,whichis a functionof frequency f,convergesas one evaluatesit usinglongerandlongerstretches ofdata. Thisinterestingquestionisthecontentofthesubjectof “powerspectrumestimation, ” and will be considered below in §13.4–§13.7. A crude answer for now is: The 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