c12-2 (Numerical Recipes in C)

504Chapter 12.Fast Fourier TransformThe discrete form of Parseval’s theorem isN−1X|hk |2 =k=0N−11 X|Hn |2N n=0(12.1.10)CITED REFERENCES AND FURTHER READING:Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall).Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (NewYork: Academic Press).12.2 Fast Fourier Transform (FFT)How much computation is involved in computing the discrete Fourier transform(12.1.7) of N points? For many years, until the mid-1960s, the standard answerwas this: Define W as the complex numberW ≡ e2πi/N(12.2.1)Then (12.1.7) can be written asHn =N−1XW nk hk(12.2.2)k=0In other words, the vector of hk ’s is multiplied by a matrix whose (n, k)th elementis the constant W to the power n × k.

The matrix multiplication produces a vectorresult whose components are the Hn ’s. This matrix multiplication evidently requiresN 2 complex multiplications, plus a smaller number of operations to generate therequired powers of W . So, the discrete Fourier transform appears to be an O(N 2 )process. These appearances are deceiving! The discrete Fourier transform can,in fact, be computed in O(N log2 N ) operations with an algorithm called the fastFourier transform, or FFT. The difference between N log2 N and N 2 is immense.With N = 106 , for example, it is the difference between, roughly, 30 seconds of CPUtime and 2 weeks of CPU time on a microsecond cycle time computer. The existenceof an FFT algorithm became generally known only in the mid-1960s, from the workof J.W.

Cooley and J.W. Tukey. Retrospectively, we now know (see [1]) that efficientmethods for computing the DFT had been independently discovered, and in somecases implemented, by as many as a dozen individuals, starting with Gauss in 1805!One “rediscovery” of the FFT, that of Danielson and Lanczos in 1942, providesone of the clearest derivations of the algorithm. Danielson and Lanczos showedthat a discrete Fourier transform of length N can be rewritten as the sum of twodiscrete Fourier transforms, each of length N/2. One of the two is formed from theSample 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).There are also discrete analogs to the convolution and correlation theorems (equations12.0.9 and 12.0.11), but we shall defer them to §13.1 and §13.2, respectively.12.2 Fast Fourier Transform (FFT)505even-numbered points of the original N , the other from the odd-numbered points.The proof is simply this:Fk =N−1Xe2πijk/N fjj=0=j=0XXN/2−1e2πik(2j)/N f2j +=XN/2−1e2πikj/(N/2)f2j + W kj=0Fke(12.2.3)j=0N/2−1=e2πik(2j+1)/N f2j+1e2πikj/(N/2)f2j+1j=0+WkFkoIn the last line, W is the same complex constant as in (12.2.1), Fke denotes the kthcomponent of the Fourier transform of length N/2 formed from the even componentsof the original fj ’s, while Fko is the corresponding transform of length N/2 formedfrom the odd components.

Notice also that k in the last line of (12.2.3) varies from0 to N , not just to N/2. Nevertheless, the transforms Fke and Fko are periodic in kwith length N/2. So each is repeated through two cycles to obtain Fk .The wonderful thing about the Danielson-Lanczos Lemma is that it can be usedrecursively. Having reduced the problem of computing Fk to that of computingFke and Fko , we can do the same reduction of Fke to the problem of computingthe transform of its N/4 even-numbered input data and N/4 odd-numbered data.In other words, we can define Fkee and Fkeo to be the discrete Fourier transformsof the points which are respectively even-even and even-odd on the successivesubdivisions of the data.Although there are ways of treating other cases, by far the easiest case is theone in which the original N is an integer power of 2.

In fact, we categoricallyrecommend that you only use FFTs with N a power of two. If the length of your dataset is not a power of two, pad it with zeros up to the next power of two. (We will givemore sophisticated suggestions in subsequent sections below.) With this restrictionon N , it is evident that we can continue applying the Danielson-Lanczos Lemmauntil we have subdivided the data all the way down to transforms of length 1.

Whatis the Fourier transform of length one? It is just the identity operation that copies itsone input number into its one output slot! In other words, for every pattern of log2 Ne’s and o’s, there is a one-point transform that is just one of the input numbers fnFkeoeeoeo···oee = fnfor some n(12.2.4)(Of course this one-point transform actually does not depend on k, since it is periodicin k with period 1.)The next trick is to figure out which value of n corresponds to which pattern ofe’s and o’s in equation (12.2.4).

The answer is: Reverse the pattern of e’s and o’s,then let e = 0 and o = 1, and you will have, in binary the value of n. Do you seewhy it works? It is because the successive subdivisions of the data into even and oddare tests of successive low-order (least significant) bits of n. This idea of bit reversalcan be exploited in a very clever way which, along with the Danielson-LanczosSample 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).XN/2−1506Chapter 12.Fast Fourier Transform000000001001001010010010011011011100100100101101101110110110111111111(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 fjand rearrange it into bit-reversed order (see Figure 12.2.1), so that the individualnumbers are in the order not of j, but of the number obtained by bit-reversing j.Then the bookkeeping on the recursive application of the Danielson-Lanczos Lemmabecomes extraordinarily simple.

The points as given are the one-point transforms.We combine adjacent pairs to get two-point transforms, then combine adjacent pairsof pairs to get 4-point transforms, and so on, until the first and second halves ofthe whole data set are combined into the final transform. Each combination takesof order N operations, and there are evidently log2 N combinations, so the wholealgorithm is of order N log2 N (assuming, as is the case, that the process of sortinginto bit-reversed order is no greater in order than N log2 N ).This, then, is the structure of an FFT algorithm: It has two sections. The firstsection sorts the data into bit-reversed order. Luckily this takes no additional storage,since it involves only swapping pairs of elements.

(If k1 is the bit reverse of k2 , thenk2 is the bit reverse of k1 .) The second section has an outer loop that is executedlog2 N times and calculates, in turn, transforms of length 2, 4, 8, . . ., N . For eachstage of this process, two nested inner loops range over the subtransforms alreadycomputed and the elements of each transform, implementing the Danielson-LanczosLemma.

The operation is made more efficient by restricting external calls fortrigonometric sines and cosines to the outer loop, where they are made only log2 Ntimes. Computation of the sines and cosines of multiple angles is through simplerecurrence 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 dataarray (data[1..2*nn]), and isign, which should be set to either ±1 and is the signof i in the exponential of equation (12.1.7). When isign is set to −1, the routinethus calculates the inverse transform (12.1.9) — except that it does not multiply bythe normalizing factor 1/N that appears in that equation.

You can do that yourself.Notice that the argument nn is the number of complex data points. The actualSample 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).00012.2 Fast Fourier Transform (FFT)507#include <math.h>#define SWAP(a,b) tempr=(a);(a)=(b);(b)=temprvoid four1(float data[], unsigned long nn, int isign)Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as 1; or replacesdata[1..2*nn] by nn times 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.