500 Chapter 12.Fast Fourier Transform PSD-per-unit-time converges to finite values at all frequencies except those 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 formalism that 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 measurements of h(ti)for a discrete set ofti's.The profound implications of this seemingly unimportant fact are the subject of the next section. CITED REFERENCES AND FURTHER READING: Champeney,D.C.1973,Fourier Transforms and Their Physical Applications(New York:Academic Press). Elliott,D.F.,and Rao,K.R.1982,Fast Transforms:Algorithms,Analyses,Applications(New York: Academic Press). 12.1 Fourier Transform of Discretely Sampled 西 令 Press. Data In the most common situations,function h(t)is sampled (i.e.,its value is 是8e recorded)at evenly spaced intervals in time.Let A denote the time interval between consecutive samples,so that the sequence of sampled values is 6 hm=h(n△)n=.,-3,-2,-1,0,1,2,3,… (12.1.1) The reciprocal of the time interval A is called the sampling rate;if A is measured in seconds,for example,then the sampling rate is the number of samples recorded per second. Sampling Theorem and Aliasing 4 E喜 Numerical Recipes 10.621 43106 For any sampling interval A,there is also a special frequency fe,called the Nyquist critical frequency,given by (outside f三2△ (12.1.2) If a sine wave of the Nyquist critical frequency is sampled at its positive peak value, 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 A.In this case the Nyquist critical frequency is just the constant 1/2. The Nyquist critical frequency is important for two related,but distinct,reasons. One is good news,and the other bad news.First the good news.It is the remarkable
500 Chapter 12. Fast Fourier Transform Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). PSD-per-unit-time converges to finite values at all frequencies except those 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 formalism that 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 measurements of h(t i) for a discrete set of ti’s. The profound implications of this seemingly unimportant fact are the subject of the next section. CITED REFERENCES AND FURTHER READING: Champeney, D.C. 1973, Fourier Transforms and Their Physical Applications (New York: Academic Press). Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (New York: 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) at evenly spaced intervals in time. Let ∆ denote the time interval between 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; if ∆ is measured in seconds, for example, then the sampling rate is the number of samples recorded per second. Sampling Theorem and Aliasing For any sampling interval ∆, there is also a special frequency f c, called the Nyquist critical frequency, given by fc ≡ 1 2∆ (12.1.2) If a sine wave of the Nyquist critical frequency is sampled at its positive peak value, 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. The Nyquist critical frequency is important for two related, but distinct, reasons. One is good news, and the other bad news. First the good news. It is the remarkable
12.1 Fourier Transform of Discretely Sampled Data 501 fact known as the sampling theorem:If a continuous function h(t).sampled at an interval A,happens to be bandwidth limited to frequencies smaller in magnitude than fe,i.e.,if H(f)=0 for all |f>fe,then the function h(t)is completely determined by its samples hn.In fact,h(t)is given explicitly by the formula sin[2rfe(t-n△)] h(t)=△ hn (12.1.3) π(t-n△) n=- 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,the sampling theorem tells us that the entire information content of the signal can be recorded by sampling it at a rate A-equal to twice the maximum frequency passed by the amplifier (cf.12.1.2). Now the bad news.The bad news concerns the effect of sampling a continuous function that is not bandwidth limited to less than the Nyquist critical frequency. 9 In that case,it turns out that all of the power spectral density that lies outside of the frequency range -fe<f fe is spuriously moved into that range.This phenomenon is called aliasing.Any frequency component outside of the frequency range (-fe,fe)is aliased (falsely translated)into that range by the very act of discrete sampling.You can readily convince yourself that two waves exp(2mifit) a%29 9 and exp(2mif2t)give the same samples at an interval A if and only if fi and f2 differ by a multiple of 1/A,which is just the width in frequency of the range (-fe,fe).There is little that you can do to remove aliased power once you have a冰 OF SCIENTIFIC 6 discretely sampled a signal.The way to overcome aliasing is to (1)know the 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 at least two points per cycle of the highest frequency present.Figure 12.1.1 illustrates these considerations. To put the best face on this,we can take the alternative point of view:If a Numerica 10621 continuous function has been competently sampled,then,when we come to estimate its Fourier transform from the discrete samples,we can assume(or rather we might 431 as well assume)that its Fourier transform is equal to zero outside of the frequency Recipes range in between-fe and fe.Then we look to the Fourier transform to tell whether 腿 the continuous function has been competently sampled(aliasing effects minimized). We do this by looking to see whether the Fourier transform is already approaching North zero as the frequency approaches fe from below,or -fe from above.If,on the contrary,the transform is going towards some finite value,then chances are that components outside of the range have been folded back over onto the critical range. Discrete Fourier Transform We now estimate the Fourier transform of a function from a finite number of its sampled points.Suppose that we have N consecutive sampled values hk≡h(tk),tk≡k△，k=0,1,2,..,N-1 (12.1.4)
12.1 Fourier Transform of Discretely Sampled Data 501 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). fact known as the sampling theorem: If a continuous function h(t), sampled at an interval ∆, happens to be bandwidth limited to frequencies smaller in magnitude than fc, i.e., if H(f)=0 for all |f| ≥ fc, then the function h(t) is completely determined by its samples hn. In fact, h(t) is given explicitly by the formula h(t)=∆ +∞ n=−∞ hn sin[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, the sampling theorem tells us that the entire information content of the signal can be recorded by sampling it at a rate ∆−1 equal to twice the maximum frequency passed by the amplifier (cf. 12.1.2). Now the bad news. The bad news concerns the effect of sampling a continuous function that is not bandwidth 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<fc is spuriously moved into that range. This phenomenon is called aliasing. Any frequency component outside of the frequency range (−fc, fc) is aliased (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πif2t) give the same samples at an interval ∆ if and only if f1 and f2 differ 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 discretely sampled a signal. The way to overcome aliasing is to (i) know the 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 at least two points per cycle of the highest frequency present. Figure 12.1.1 illustrates these considerations. To put the best face on this, we can take the alternative point of view: If a continuous function has been competently sampled, then, when we come to estimate its Fourier transform from the discrete samples, we can assume (or rather we might as well assume) that its Fourier transform is equal to zero outside of the frequency range in between −fc and fc. Then we look to the Fourier transform to tell whether the continuous function has been competently sampled (aliasing effects minimized). We do this by looking to see whether the Fourier transform is already approaching zero as the frequency approaches f c from below, or −fc from above. If, on the contrary, the transform is going towards some finite value, then chances are that components outside of the range have been folded back over onto the critical range. Discrete Fourier Transform We now estimate the Fourier transform of a function from a finite number of its sampled points. Suppose that we have N consecutive sampled values hk ≡ h(tk), tk ≡ k∆, k = 0, 1, 2,...,N − 1 (12.1.4)
502 Chapter 12.Fast Fourier Transform h(0 (a) H http://www.r. Permission is read able files .com or call granted for (including this one) 19881992 11-800.872 (b) from NUMERICAL RECIPES IN C: /Cambridge aliased Fourier transform true Fourier transform (North America server computer, tusers to make one paper UnN电.t THE ART 2△ strictly proh Programs (c) 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 to dir limited but has finite amplitude for all frequencies.If the original function is sampled with a sampling interval A,as in (a),then the Fourier transform (c)is defined 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. v@cam so that the sampling interval is A.To make things simpler,let us also suppose that N is even.If the function h(t)is nonzero only in a finite interval of time,then that whole interval of time is supposed to be contained in the range of the N points 乳 1881892 OF SCIENTIFIC COMPUTING(ISBN .Further reproduction, Numerical Recipes 10-621 43108 given.Alternatively,if the function h(t)goes on forever,then the sampled points are supposed to be at least"typical"of what h(t)looks like at all other times. (outside With N numbers of input,we will evidently be able to produce no more than N independent numbers of output.So,instead of trying to estimate the Fourier North Software. transform H(f)at all values of f in the range -fe to fe,let us seek estimates only at the discrete values visit website fn=N△' n=- (12.1.5) machine 292 The extreme values of n in(12.1.5)correspond exactly to the lower and upper limits 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 n in (12.1.5):it will turn out that the two extreme values of n are not independent(in fact they are equal),but all the others are.This reduces the count to N
502 Chapter 12. Fast Fourier Transform Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). h(t) t (a) f 0 H( f ) (b) (c) aliased Fourier transform true Fourier transform 0 H( 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 defined 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 N is even. If the function h(t) is nonzero only in a finite interval of time, then that whole interval of time is supposed to be contained in the range of the N points given. Alternatively, if the function h(t) goes on forever, then the sampled points are supposed to be at least “typical” of what h(t) looks like at all other times. With N numbers of input, we will evidently be able to produce no more than N independent numbers of output. So, instead of trying to estimate the Fourier transform H(f) at all values of f in the range −fc to fc, let us seek estimates only at the discrete values fn ≡ n N∆, n = −N 2 ,..., N 2 (12.1.5) The extreme values of n in (12.1.5) correspond exactly to the lower and upper limits 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 n in (12.1.5); it will turn out that the two extreme values of n are not independent (in fact they are equal), but all the others are. This reduces the count to N
12.1 Fourier Transform of Discretely Sampled Data 503 The remaining step is to approximate the integral in(12.0.1)by a discrete sum: N-1 N-1 H(fn)= h(t)e2≈hke2mifn△=△hx e2mikn/N k=0 k=0 (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, N-1 Hn≡ hk e2nikn/N (12.1.7) k=0 素 The discrete Fourier transform maps N complex numbers(the hk's)into N complex numbers(the Hn's).It does not depend on any dimensional parameter,such as the time scale A.The relation (12.1.6)between the discrete Fourier transform of a set of numbers and their continuous Fourier transform when they are viewed as samples of a continuous function sampled at an interval A can be rewritten as H(fn)≈△Hn (12.1.8) 鸟6 9 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 in n,with IENTIFIC period N.Therefore,Hn =HN-n n 1,2,....With this conversion in mind, one generally lets the n in Hn vary from 0 to N-1 (one complete period).Then n and k (in h)vary exactly over the same range,so the mapping of N numbers into N numbers is manifest.When this convention is followed,you must remember that zero frequency corresponds to n =0,positive frequencies 0<f<fe correspond to values 1 nN/2-1,while negative frequencies-fe<f<0 correspond to N/2+1<n N-1.The value nN/2 corresponds to both f fcand f=-fe. The discrete Fourier transform has symmetry properties almost exactly the same Recipes Numerica 10.621 as the continuous Fourier transform.For example,all the symmetries in the table 431 following equation (12.0.3)hold if we read hk for h(t),Hn for H(f),and HN-n (outside Recipes for H(-f).(Likewise,"even"and"odd"in time refer to whether the values h k at k and N-k are identical or the negative of each other.) The formula for the discrete imverse Fourier transform,which recovers the set North of hk's exactly from the Hn's is: N-1 hk一 e-2xikn/N (12.1.9) =0 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 routine for calculating discrete Fourier transforms can also,with slight modification, calculate the inverse transforms
12.1 Fourier Transform of Discretely Sampled Data 503 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The remaining step is to approximate the integral in (12.0.1) by a discrete sum: H(fn) = ∞ −∞ h(t)e2πifnt dt ≈ N −1 k=0 hk e2πifntk∆=∆ N −1 k=0 hk e2π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 k=0 hk e2πikn/N (12.1.7) The discrete Fourier transform maps N complex numbers (the h k’s) into N complex 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 of numbers and their continuous Fourier transform when they are viewed as samples of a continuous function sampled at an interval ∆ can be rewritten as H(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 in n, with period N. Therefore, H−n = HN−n n = 1, 2,... . With this conversion in mind, one generally lets the n in Hn vary from 0 to N − 1 (one complete period). Then n and k (in hk) vary exactly over the same range, so the mapping of N numbers into N numbers is manifest. When this convention is followed, you must remember that zero frequency corresponds to n = 0, positive frequencies 0 <f<f c 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/2 corresponds to both f = f c and f = −fc. The discrete Fourier transform has symmetry properties almost exactly the same as the continuous Fourier transform. For example, all the symmetries in the table following equation (12.0.3) hold if we read h k for h(t), Hn for H(f), and HN−n for H(−f). (Likewise, “even” and “odd” in time refer to whether the values h k at k and N − k are identical or the negative of each other.) The formula for the discrete inverse Fourier transform, which recovers the set of hk’s exactly from the Hn’s is: hk = 1 N N −1 n=0 Hn e−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 routine for calculating discrete Fourier transforms can also, with slight modification, calculate the inverse transforms.
504 Chapter 12.Fast Fourier Transform The discrete form of Parseval's theorem is N-1 N-1 (12.1.10) k=0 刀=0 There are also discrete analogs to the convolution and correlation theorems(equations 12.0.9 and 12.0.11),but we shall defer them to $13.1 and $13.2,respectively. 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 (New York: Academic Press). 、意轮 12.2 Fast Fourier Transform(FFT) RECIPES I 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 answer 平空> Press. was this:Define W as the complex number W三e2m/N (12.2.1) Then (12.1.7)can be written as IENTIFIC 61 N-1 Hn= wnkhk (12.2.2) k=0 In other words,the vector of hk's is multiplied by a matrix whose(n,k)th element is the constant W to the power n x k.The matrix multiplication produces a vector result whose components are the A,'s.This matrix multiplication evidently requires Numerica 10.621 N2 complex multiplications,plus a smaller number of operations to generate the 431 required powers of W.So,the discrete Fourier transform appears to be an O(N2) Recipes 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 fast (outside Fourier transform,or FFT.The difference between Nlog2 N and N2 is immense. North With N =106,for example,it is the difference between,roughly,30 seconds of CPU time and 2 weeks of CPU time on a microsecond cycle time computer.The existence of an FFT algorithm became generally known only in the mid-1960s,from the work of J.W.Cooley and J.W.Tukey.Retrospectively,we now know (see [11)that efficient methods for computing the DFT had been independently discovered.and in some cases 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,provides one of the clearest derivations of the algorithm.Danielson and Lanczos showed that a discrete Fourier transform of length N can be rewritten as the sum of two discrete Fourier transforms,each of length N/2.One of the two is formed from the
504 Chapter 12. Fast Fourier Transform Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The discrete form of Parseval’s theorem is N −1 k=0 |hk| 2 = 1 N N −1 n=0 |Hn| 2 (12.1.10) There are also discrete analogs to the convolution and correlation theorems (equations 12.0.9 and 12.0.11), but we shall defer them to §13.1 and §13.2, respectively. 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 (New York: 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 answer was this: Define W as the complex number W ≡ e2πi/N (12.2.1) Then (12.1.7) can be written as Hn = N −1 k=0 Wnkhk (12.2.2) In other words, the vector of hk’s is multiplied by a matrix whose (n, k)th element is the constant W to the power n × k. The matrix multiplication produces a vector result whose components are the Hn’s. This matrix multiplication evidently requires N2 complex multiplications, plus a smaller number of operations to generate the required 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 fast Fourier transform, or FFT. The difference between N log2 N and N2 is immense. With N = 106, for example, it is the difference between, roughly, 30 seconds of CPU time and 2 weeks of CPU time on a microsecond cycle time computer. The existence of an FFT algorithm became generally known only in the mid-1960s, from the work of J.W. Cooley and J.W. Tukey. Retrospectively, we now know (see [1]) that efficient methods for computing the DFT had been independently discovered, and in some cases 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, provides one of the clearest derivations of the algorithm. Danielson and Lanczos showed that a discrete Fourier transform of length N can be rewritten as the sum of two discrete Fourier transforms, each of length N/2. One of the two is formed from the