c12-0 (Numerical Recipes in C)

Fast Fourier Transform12.0 IntroductionA very large class of important computational problems falls under the generalrubric of “Fourier transform methods” or “spectral methods.” For some of theseproblems, the Fourier transform is simply an efficient computational tool foraccomplishing certain common manipulations of data. In other cases, we haveproblems for which the Fourier transform (or the related “power spectrum”) is itselfof intrinsic interest. These two kinds of problems share a common methodology.Largely for historical reasons the literature on Fourier and spectral methods hasbeen disjoint from the literature on “classical” numerical analysis. Nowadays there isno justification for such a split. Fourier methods are commonplace in research and weshall not treat them as specialized or arcane.

At the same time, we realize that manycomputer users have had relatively less experience with this field than with, say,differential equations or numerical integration. Therefore our summary of analyticalresults will be more complete. Numerical algorithms, per se, begin in §12.2. Variousapplications of Fourier transform methods are discussed in Chapter 13.A physical process can be described either in the time domain, by the values ofsome quantity h as 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 complexnumber 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) and H(f) as beingtwo different representations of the same function. One goes back and forth betweenthese two representations by means of the Fourier transform equations,Z∞H(f) =Zh(t)e2πift dt−∞∞(12.0.1)−2πiftH(f)eh(t) =df−∞If t is measured in seconds, then f in equation (12.0.1) is in cycles per second,or Hertz (the unit of frequency).

However, the equations work with other units too.If h is a function of position x (in meters), H will be a function of inverse wavelength(cycles per meter), and so on. If you are trained as a physicist or mathematician, youare probably more used to using angular frequency ω, which is given in radians persec. The relation between ω and f, H(ω) and H(f) isω ≡ 2πfH(ω) ≡ [H(f)] f=ω/2π496(12.0.2)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Chapter 12.49712.0 Introductionand equation (12.0.1) looks like thisZ∞h(t)eiωt dtH(ω) =−∞Z(12.0.3)∞−iωtH(ω)edω−∞We were raised on the ω-convention, but we changed! There are fewer factors of2π to remember if you use the f-convention, especially when we get to discretelysampled data in §12.1.From equation (12.0.1) it is evident at once that Fourier transformation is alinear operation. The transform of the sum of two functions is equal to the sum ofthe transforms.

The transform of a constant times a function is that same constanttimes the transform of the function.In the time domain, function h(t) may happen to have one or more specialsymmetries It might be purely real or purely imaginary or it might be even,h(t) = h(−t), or odd, h(t) = −h(−t). In the frequency domain, these symmetrieslead to relationships between H(f) and H(−f). The following table gives thecorrespondence between symmetries in the two domains:If . . .h(t) is realh(t) is imaginaryh(t) is evenh(t) is oddh(t) is real and evenh(t) is real and oddh(t) is imaginary and evenh(t) is imaginary and oddthen .

. .H(−f) = [H(f)]*H(−f) = −[H(f)]*H(−f) = H(f) [i.e., H(f) is even]H(−f) = −H(f) [i.e., H(f) is odd]H(f) is real and evenH(f) is imaginary and oddH(f) is imaginary and evenH(f) is real and oddIn subsequent sections we shall see how to use these symmetries to increasecomputational efficiency.Here are some other elementary properties of the Fourier transform. (We’ll usethe “⇐⇒” symbol to indicate transform pairs.) Ifh(t) ⇐⇒ H(f)is such a pair, then other transform pairs areh(at) ⇐⇒f1H( )|a|a“time scaling”(12.0.4)1 th( ) ⇐⇒ H(bf)|b| b“frequency scaling”(12.0.5)h(t − t0 ) ⇐⇒ H(f) e2πift0h(t) e−2πif0 t ⇐⇒ H(f − f0 )“time shifting”“frequency shifting”(12.0.6)(12.0.7)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).1h(t) =2π498Chapter 12.Fast Fourier TransformWith two functions h(t) and g(t), and their corresponding Fourier transformsH(f) and G(f), we can form two combinations of special interest. The convolutionof the two functions, denoted g ∗ h, is defined byZg∗h≡∞g(τ )h(t − τ ) dτ(12.0.8)Note that g ∗ h is a function in the time domain and that g ∗ h = h ∗ g. It turns outthat the function g ∗ h is one member of a simple transform pairg ∗ h ⇐⇒ G(f)H(f)“Convolution Theorem”(12.0.9)In other words, the Fourier transform of the convolution is just the product of theindividual Fourier transforms.The correlation of two functions, denoted Corr(g, h), is defined byZCorr(g, h) ≡∞g(τ + t)h(τ ) dτ(12.0.10)−∞The correlation is a function of t, which is called the lag.

It therefore lies in the timedomain, and it turns out to be one member of the transform pair:Corr(g, h) ⇐⇒ G(f)H*(f)“Correlation Theorem”(12.0.11)[More generally, the second member of the pair is G(f)H(−f), but we are restrictingourselves to the usual case in which g and h are real functions, so we take the liberty ofsetting H(−f) = H*(f).] This result shows that multiplying the Fourier transformof one function by the complex conjugate of the Fourier transform of the other givesthe Fourier transform of their correlation.

The correlation of a function with itself iscalled its autocorrelation. In this case (12.0.11) becomes the transform pairCorr(g, g) ⇐⇒ |G(f)|2“Wiener-Khinchin Theorem”(12.0.12)The total power in a signal is the same whether we compute it in the timedomain or in the frequency domain. This result is known as Parseval’s theorem:ZTotal Power ≡∞−∞Z2|h(t)| dt =∞−∞2|H(f)| df(12.0.13)Frequently one wants to know “how much power” is contained in the frequencyinterval between f and f + df.

In such circumstances one does not usuallydistinguish between positive and negative f, but rather regards f as varying from 0(“zero frequency” or D.C.) to +∞. In such cases, one defines the one-sided powerspectral density (PSD) of the function h asPh (f) ≡ |H(f)|2 + |H(−f)| 20≤f <∞(12.0.14)so that the total power is just the integral of Ph (f) from f = 0 to f = ∞. When the2function h(t) is real, then the two terms in (12.0.14) are equal, so Ph (f) = 2 |H(f)| .Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.