Web Appendix M - The Fast Fourier Transform

M. J. Roberts - 2/18/07Web Appendix M - The Fast FourierTransformM.1 IntroductionThe forward DFT is defined byX k N F 1 x n en 0 j 2 nk / N F.A straightforward way of computing the DFT would be by the following algorithm(written in MATLAB) which directly implements the operations indicated above.%(Acquire the input data in an array, x, with NF elements.).%%Initialize the DFT array to a column vector of zeros.%X zeros(NF,1) ;%%Compute the Xk ’s in a nested, double for loop.%for k 0:NF-1for n 0:NF-1X(k 1) X(k 1) x(n 1)*exp(-j*2*pi*n*k/NF) ;endend.(By the way, one should never actually write this program in MATLAB because the DFTis already built in to MATLAB as an intrinsic function called fft.)The computation of a DFT using this algorithm requires N 2 complex multiply-addoperations. Therefore the number of computations increases as the square of the numberof elements in the input vector which is being transformed. In 1965 James Cooley andJohn Tukey popularized an algorithm which is much more efficient in computing time forlarge input arrays whose length is an integer power of 2. This algorithm for computing theDFT is the so-called fast Fourier transform or just FFT.M-1

M. J. Roberts - 2/18/07James W. Cooley (1926-)John Wilder Tukey (1915-2000)James Cooley received his Ph.D. in applied mathematics from Columbia University in 1961.Cooley was a pioneer in the digital signal processing field, having developed, with John Tukey, the fastFourier transform. He developed the FFT through mathematical theory and applications, and has helpedmake it more widely available by devising algorithms for scientific and engineering applications.John Tukey received his Ph.D. from Princeton in mathematics in 1939. He worked at Bell Labsfrom 1945 to 1970. He developed new techniques in data analysis. He developed graphing and plottingmethods that now appear in standard statistics texts. He wrote many publications on time series analysisand other aspects of digital signal processing that are now very important in engineering and science. Hedeveloped, along with James Cooley the fast Fourier transform algorithm. He is credited with havingcoined , as a contraction of “binary digit”, the word “bit”, the smallest unit of information used by acomputer.M.2 The Decimation-in-Frequency FFTThere are two commonly-used FFT algorithms called the decimation-in-time(DIT) algorithm and the decimation-in-frequency (DIF) algorithm. We will consider firstthe DIF algorithm. We will see in this development that the FFT algorithm is a goodexample of the problem-solving motto “divide and conquer” which is so effective in LTIsystem analysis. At every stage of calculation in the FFT algorithm we defer the actualDFT calculation by splitting the data into two halves and finding the DFT of the halvesseparately. We do this until we get down to multiple 2-point DFT’s which cannot befurther subdivided. The DIF FFT algorithm can also be considered an example of “topdown” problem solving in which we attack the overall problem and, as necessary, assignparts of the analysis as sub-problems until we get to sub-problems that are simple enoughto solve directly.The DFT formula for an N -point transform isN 1X k x n e j 2 nk / N .n 0M-2(M.1)

M. J. Roberts - 2/18/07Similarly to what we did in the development of the DTFS in Chapter 4, it is convenient touse the notation WN e j 2 / N . (In chapter 4 it was defined as WN e j 2 / N . Because ofthe minus sign in the exponent of e in (M.1) it is more convenient here to define it withthe sign of the exponent changed.) ThenN 1X k x n WNnk .(M.2)n 0W can be thought of as a vector in the complex plane of length 1 and angle 2 / N(Figure M-1).Im(W)Re(W)-2π/N1Figure M-1 W as a vector in the complex planeWe will restrict this development to N 2 , with being a positive integer becauseonly in that case is the FFT algorithm optimally efficient. To compute all N values ofX k , 0 k N directly from (M.2) requires N N 1 complex additions and N 2complex multiplications. (These numbers are based on a simple calculation treating every()calculation the same. Some of the multiplications are by WN0 which is obviously 1. Inthose cases we could reduce the number of multiplications but, for large N which is themost practical case, the number of multiplications would not change much, would be morecomplicated to calculate and writing a computer program to take advantage of theelimination of these multiplications would complicate the program. So we will stick withthis, admittedly oversimplified, estimate.)The FFT algorithm is very intricate and recursive so it helps to start with a smallexample and then to expand the concepts to the general case. We will start with thesmallest possible DFT N 2 , 1 , WN e j 2 / 2 1 .X k 1 x n Wn0 0nk2and X 0 W 0 W 0 x 0 1 1 x 0 2 20 1 X 1 W2 W2 x 1 1 1 x 1 We can realize a 2-point DFT as a 2-input-2-output system (Figure M-2).M-3

M. J. Roberts - 2/18/072-Point DFTx[0]X[0]x[1]X[1]-1Figure M-2 A 2-point DFTNow we will step up to the next level, a 4-point DFT.313n 0n 0n 2X k x n W4nk x n W4nk x n W4nk1X k x n Wnk4n 011 x n 2 W4()n 2 kn 0()X k x n x n 2 W42k W4nkn 0W42k e( ) j 2 2k / 4( ) e j k 1k1(( ))X k x n 1 x n 2 W4nkn 0kLet q be any integer. Then 12q2nqX 2q x n 1 x n 2 W4 x n x n 2 W42nq n 0 n 0 1W42nq () j 2 2nq / 4(( )1) e j 2 nq / 2 W2nq1()(X 2q x n x n 2 W2nq 2-point DFT of x n x n 2 n 0)Similarly,1()()X 2q 1 x n x n 2 W4nW2nq 2-point DFT of x n x n 2 W4nn 0M-4

M. J. Roberts - 2/18/074-Point FTW41x[3]X[3]-1Figure M-3 A 4-point DIF DFT computed using two 2-point DFT’sNow consider the general case, an N-point DFT.N 1N / 2 1n 0n 0X k x n WNnk X k N / 2 1 n 0X k WNNk / 2 e() j 2 Nk / 2 / Nx n WNnk x n WNnk N / 2 1 n 0N 1 n N / 2)n N / 2 k ( x n x n N / 2 WNk / 2Nn 0( )x n WNnkx n N / 2 WN(N / 2 1 e j k 1X k )WnkNkN / 2 1 n 0( x n ( 1) x n N / 2 )WknkNLet q be any integer. ThenX 2q N / 2 1 n 0 2q x n 1 x n N / 2 WN2nq 1 ( x nN / 2 1n 0( )) x n N / 2 WN2nqM-5

M. J. Roberts - 2/18/07WN2nq () j 2 2nq / N e( j 2 nq / N / 2) WNnq/ 2 ( x n x n N / 2 )WN / 2 1X 2q nqN /2(n 0)())() N / 2 -point DFT of x n x n N / 2 Similarly,(X 2q 1 N / 2 -point DFT of x n x n N / 2 WNnN-Point DFTx[0]X[0]X[2].N/2-pointDFT.x[1]x[N/2 - 1]X[N - 2]WN0-1WN/2-pointDFTWNN/2-1X[3]Odd.-1.x[N/2 1]X[1]1N.x[N/2]x[N - 1]EvenX[N - 1]-1Figure M-4 An N-point DIF DFT computed using two N/2-point DFT’sThe system in Figure M-4 illustrates the recursive nature of the FFT algorithm. Insteadof actually computing an N-point DFT directly, split the data into two halves, form twolinear combinations of those N/2-point data sets and the compute two N/2-point DFT’s.The even-indexed DFT data points X 0 , X 2 , , X N 2 are computed in the tophalf and the odd-indexed data points X 1 , X 3 , , X N 1 are computed in thebottom half. N complex additions and N/2 complex multiplications are needed to preparethe N input data for the two N/2-point DFT’s. Then each of the two N/2-point DFT’s willrequire N/2 complex additions and N/4 complex multiplications to prepare for the fourN/4-point DFT’s. This subdivision continues until we get to N/2 2-point DFT’s, each ofwhich requires two complex additions and one complex multiplication. Suppose N iseight. Then the numbers of complex additions and complex multiplications are4 A 2 M 4 2 A M 24 A 12 M8A 4M 2 8-point level4-point level2-point levelM-6

M. J. Roberts - 2/18/07where A is the number of complex additions and M is the number of complexmultiplications. Generalizing to N points, at each level there are N complex additions andN/2 complex multiplications and there are levels. Therefore, in general, A N andM N / 2 . Table M-1 summarizes the number of arithmetic operations for severalpowers of two for both the original DFT formula, (M.2), and the DIF FFT algorithm.Table M-1 Numbers of additions and multiplications and ratios for several N’s NADFTM DFTAFFTM FFTADFT / AFFTM DFT / M 56.8102.364113.8204.88 256652806553620489 512 261632 262144 460810 1024 1047552 1048576 10240So, by using this technique we significantly reduce the number of complexadditions and multiplications, especially for large N. It is important to emphasize thatthese speed improvement factors do not apply if is not an integer. For this reason,practically all actual DFT analysis is done with the FFT using data vector lengths whichare an integer power of 2. (In MATLAB if the input vector is an integer power of 2 inlength, the algorithm used in the MATLAB function fft is the FFT algorithm justdiscussed. If it is not an integer power of 2 in length, the DFT is still computed but thespeed suffers because a less efficient algorithm must be used.)The output data from the FFT algorithm don’t come out in natural order but reordering them is simple and fast. The scrambling is illustrated for N 8 in Figure M-5with the natural order on the left and the scrambled order on the right.012345670246135704261537Figure M-5 Scrambling of the output data order for N 8M-7

M. J. Roberts - 2/18/07The process of unscrambling the output data order is not as complicated as it may look.It can be done very quickly and efficiently by expressing the N indices as -bit binarynumbers, and then simply reversing the bit order (Figure 00010110001101011111000010100110001011101111Figure M-6 Bit-order reversal to correctly identify the indices of the output dataSo, for example, the second data point in an 8-point output array is actually X 4 notX 1 .M.3 The Decimation-in-Time FFTProbably the easiest way to introduce the DIT FFT algorithm is to take the 4-pointDIF FFT and re-order the input and output data while keeping all the connections thesame (Figure M-7).4-point DIF -point DIT -8

M. J. Roberts - 2/18/07Figure M-7 Converting from a 4-point DIF FFT to a 4-point DIT FFTThe DIT FFT is the dual of the DIF FFT. It begins by reordering the input data by bitreversing the indices and then builds up from 2-point DFT’s to the final N-point DFTinstead of decomposing an N-point FFT into multiple 2-point DFT’s. If the DIF FFT isan example of “top-down” problem solving, then the DIT FFT must be an example of“bottom-up” problem solving. We begin the development of the DIT FFT by separatingthe input data into its even-indexed and odd-indexed values.N 1X k x n WNnk x n WnkNn 0n even x n WnkNn oddWe can make the transformations n 2q in the even summation and n 2q 1 in theodd summation yieldingX k N / 2 1 q 0x 2q WN2qk N / 2 1 q 0x 2q 1 WN2qkWNk .Using the previously derived WN2nq WNnq/ 2 ,X k N / 2 1 q 0x 2q WqkN /2 WkNN / 2 1 q 0x 2q 1 WNqk/ 2 .Let the even-indexed points in x be x e n and let the odd-indexed points in x be x o n .Then()X k N / 2 point DFT of x e n WNk N / 2 point DFT of x o n .In computing the values of X k for N / 2 k N we need values of the two N/2-pointDFT’s which have been computed in the range 0 k N / 2 . Since the DFT is periodicthe N/2-point DFT at any index k is the same as it is at k mN / 2 where m is anyinteger. Therefore the needed values at indices N / 2 k N can come directly fromindices 0 k N / 2 which have already been computed. The form of a general N-pointDIT DFT is illustrated in Figure M-8.M-9