The Fast Fourier Transform and Its Applications
(1969)

The advent of the fast Fourier transform method has greatly extended our ability to implement Fourier methods on digital computers. A description of the alogorithm and its programming is given here and followed by a theorem relating its operands, the finite sample sequences, to the continuous functions they often are intended to approximate. An analysis of the error due to discrete sampling over finite ranges is given in terms of aliasing. Procedures for computing Fourier integrals, convolutions and lagged products are outlined. HISTORICAL BACKGROUND T HE FAST Fourier transform algorithm has an interesting history which has been described in [3]. Time does not permit repeating this history here in detail. The essentials are, however, that until the recent publication of fast Fourier transform methods, computer programs were using up hundreds of hours of computer time with procedures requiring something proportional to N2 operations to compute Fourier transforms of N data points. It is not surprising then that the "new" methods requiring a number of operations proportional to N log N received considerable attention and led to revisions in computer programs and in problemn-solving techniques using Fourier methods. It was discovered later that the base 2 form of the fast Fourier transform algorithm had been published many years ago by Runge and Konig [10] and by Stumpff [12], [13]. These are authors whose works are widely read and their papers certainly were used by those computing Fourier series. How then could these important algorithms have gone unnoticed? The answer is that the papers of Runge, K6nig, and Stumpif described primarily how one could use symmetries of the sine-cosine functions to reduce the amount of computation by factors of 4, 8, or even more. Relatively small portions of these papers mentioned the successive doubling algorithm which permitted one to take two Fourier analyses of N-point samples of data and combine them in N operations to obtain an analysis of a 2N-point sampling of the same data. Successive application of this algorithm obviously yields an N-point Fourier analysis in 10g2 N doublings, and therefore, takes N log2 N operations. Thus, while the computational method using symmetries reduced the proportionality factor in the KN2 operations required to transform an N-point sequence, the method based on the doubling algorithm took a number of operations proportional to N log2 N. Manuscript received August 16, 1968. The authors are with the IBM Watson Research Center, Yorktown Heights, N.Y. It is most likely that with the relatively small values of N used in preelectronic computer days, the former methods were easier to use and took fewer operations. Consequently, the methods requiring N log N operations were neglected. With the arrival of electronic computers capable of doing calculations of Fourier transforms with large values of N, the N log N methods were overlooked and the well-known hand calculator methods were programmed for the computers. Perhaps there is something to be learned from this experience, namely, that there may exist numerical methods in the older literature which should be reappraised whenever computing devices undergo radical changes. FAST FOURIER ALGORITHMS In the interest of coherent presentation, the definitions and procedures utilized herein will be developed in blocks and shown as figures. Thus, the discrete Fourier series is defined in Fig. 1 with A (n) being a sequence which gives the complex Fourier amplitudes as a function of frequency n. The X(j), j=0, 1, * * *, N-1 are regarded here as a complex sequence, and in a problem may represent a sampling of a signal at N sampling points. Also, (2iri' 2r 2r WN= exp .= cos-+ isinN1 N N is the principle Nth root of unity and if we substitute the expression for WN in terms of sines and cosines, we obtain the perhaps more familiar sine-cosine Fourier series. Complex Fourier series is employed for ease of notation and derivation of formulas. One should note next the inversion formula in Fig. 1, giving the A (n)'s in terms of the X(j)'s. Since A (n) is also a Fourier series, an algorithm or a program for computing the A (n)'s from the X(j)'s can be used to compute the X(j)'s from the A(n)'s. Since WNN= 1, the exponent of WN is to be interpreted modulo N. This leads to an essential property of the sequences X(j) and A (n), i.e., that they are periodic functions of j and n, respectively, with period N. It is shown in Fig. 2 that when N is a product, N= rs, the Fourier series can be calculated in a two-stage process. This is done just as though the sequences A (n) and X(j) were defined on two-dimensional rXs arrays with the array indices (ji, jo) and (n1, no) being defined as shown. When we substitute for j and n in WNin, and jn is reduced modulo N, it is found that the series can be 27 Authorized licensed use limited to: Peking University. Downloaded on March 24, 2009 at 01:10 from IEEE Xplore. Restrictions apply. IEEE TRANSACTIONS ON EDUCATION, MARCH 1969 Discrete Fourier Series: N-1 X(j) E A(n)wjn n,O N where WN exp(27riIN) Inversion Formula: N-1 N j.0 N Let us write: X(tjt4-A(n) Periodicity: WNN * I, Wj'N+N.