c12-1 (779557), страница 2
Текст из файла (страница 2)
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).(a)50312.1 Fourier Transform of Discretely Sampled DataThe remaining step is to approximate the integral in (12.0.1) by a discrete sum:ZH(fn ) =∞−∞h(t)e2πifn t dt ≈N−1Xhk e2πifn tk ∆ = ∆k=0N−1Xhk e2πikn/Nk=0(12.1.6)Hn ≡N−1Xhk e2πikn/N(12.1.7)k=0The discrete Fourier transform maps N complex numbers (the hk ’s) into N complexnumbers (the Hn ’s). It does not depend on any dimensional parameter, such as thetime scale ∆.
The relation (12.1.6) between the discrete Fourier transform of a setof numbers and their continuous Fourier transform when they are viewed as samplesof a continuous function sampled at an interval ∆ can be rewritten asH(fn ) ≈ ∆Hn(12.1.8)where fn is given by (12.1.5).Up to now we have taken the view that the index n in (12.1.7) varies from−N/2 to N/2 (cf. 12.1.5). You can easily see, however, that (12.1.7) is periodic inn, with period N .
Therefore, H−n = HN−n n = 1, 2, . . .. With this conversionin mind, one generally lets the n in Hn vary from 0 to N − 1 (one completeperiod). Then n and k (in hk ) vary exactly over the same range, so the mappingof N numbers into N numbers is manifest.
When this convention is followed,you must remember that zero frequency corresponds to n = 0, positive frequencies0 < f < fc correspond to values 1 ≤ n ≤ N/2 − 1, while negative frequencies−fc < f < 0 correspond to N/2 + 1 ≤ n ≤ N − 1. The value n = N/2corresponds to both f = fc and f = −fc .The discrete Fourier transform has symmetry properties almost exactly the sameas the continuous Fourier transform. For example, all the symmetries in the tablefollowing equation (12.0.3) hold if we read hk for h(t), Hn for H(f), and HN−nfor H(−f). (Likewise, “even” and “odd” in time refer to whether the values hk at kand N − k are identical or the negative of each other.)The formula for the discrete inverse Fourier transform, which recovers the setof hk ’s exactly from the Hn ’s is:hk =N−11 XHn e−2πikn/NN n=0(12.1.9)Notice that the only differences between (12.1.9) and (12.1.7) are (i) changing thesign in the exponential, and (ii) dividing the answer by N .
This means that aroutine for calculating discrete Fourier transforms can also, with slight modification,calculate the inverse transforms.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).Here equations (12.1.4) and (12.1.5) have been used in the final equality. The finalsummation in equation (12.1.6) is called the discrete Fourier transform of the Npoints hk . Let us denote it by Hn ,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..