The Fast Fourier Transform

MASSACHUSETTS INSTITUTE OF TECHNOLOGY DEPARTMENT OF MECHANICAL ENGINEERING 2.161 Signal Processing - Continuous and Discrete The Fast Fourier Transform Introduction: Although the DFT was known for many decades, its utility was severely limited because of the computational burden. The calculation of the DFT of an input sequence of an N length sequence {fn } Fm = N −1 � fn e−j 2πmn N , m = 0, . . . , N − 1 (1) n=0 requires N complex multiplications to compute each on the N values, Fm , for a total of N 2 multi­ plications. Early digital computers had neither fixed-point nor floating-point hardware multipliers, and multiplication was performed by binary shift-and-add software algorithms. Multiplication was therefore a computationally “expensive” and time consuming operation, rendering machine com­ putation of the DFT impractical for common usage. In 1965 Cooley and Tukey [1] of IBM labs published a ground-breaking paper that presented a computational algorithm for the DFT that required just a small fraction of the complex multiplica­ tions in Eq. (1). The Fast Fourier Transform (FFT) has revolutionized digital signal processing by allowing practical fast frequency domain implementation of processing algorithms. (It is interesting to note in retrospect that the techniques used in the FFT can be found all the way back to Gauss, although the significance was not recognized until Cooley and Tukey). There are many variations on the FFT algorithm. This handout examines just one of them: the Radix-2 FFT with decimation in time, which is probably the most commonly used FFT. The term Radix-2 refers to the limitation that the sample length N must be an integer power of 2, while decimation in time means that the sequence {fn } must be re-ordered before applying the algorithm. The Radix-2 FFT with Decimation in Time: We start by writing the DFT of Eq. (1) as Fm = N −1 � fn e−j 2πmn N n=0 = N −1 � fn WNmn , m = 0, . . . , N − 1 (2) n=0 where WN = e−j2π/N . We also note the following periodic and symmetry properties of WN k(N −n) • WN = WN−kn = WNkn (complex conjugate symmetry), k(n+N ) • WNkn = WN n−N/2 • WNn = −WN n(k+N ) = WN (periodicity in n and k), for n ≥ N/2. The FFT recognizes that these properties render many of the N 2 complex multiplications in Eq. (1) redundant. We start by assuming that the input sequence length N is even. We then ask whether any computational efficiency might be gained from splitting the calculation of {Fm } into two subcomputations, each of length N/2, involving the even samples, {f2n }, and odd samples {f2n+1 } for 1 n = 0 . . . N/2 − 1: Fm = = = P� −1 n=0 P� −1 n=0 P� −1 f2n WN2mn + P� −1 m(2n+1) f2n+1 WN n=0 f2n WN2mn + WNm f2n WPmn + WNm n=0 P� −1 f2n+1 WN2mn n=0 P� −1 f2n+1 WPmn n=0 = Am + WNm Bm , m = 0, . . . , N − 1. (3) where P = N/2, {Am } is a DFT of length N/2, based on the even sample points, and similarly {Bm } is a DFT of length N/2 based on the odd sample points of {fn }. We also note from the properties of the DFT that both {Am } and {Bm } are periodic with period N/2, that is Am+N/2 = Am , and Bm+N/2 = Bm so that Fm = Am + WNm Bm Fm = Am−N/2 − m−N/2 WN Bm−N/2 for m = 0 . . . (N/2 − 1), and (4) for m = N/2 . . . (N − 1) (5) Equations (4) and (5) show that a DFT of length N may be synthesized by combining two shorter DFTs from an even/odd decomposition of the original data set. For example, if N = 8, F3 and F7 are simply related: F3 = A3 + W83 B3 F7 = A7 + W87 B7 = A3 − W83 B3 so that N/2 multiplications are required to combine the two sets. Each of the two shorter DFTs requires (N/2)2 complex multiplications, therfore the total required is N 2 /2+N/2 < N 2 , for N > 2, indicating a computational saving. A modified discrete form of Mason’s signal-flow graph is commonly used to display the algo­ rithmic structure of Eq. (3). Figure 1 shows the signal-flow graph - consisting of a network of nodes connected by line segments. The algorithm works from left to right, and each right-hand node is assigned a value that is the weighted sum of the connected left-hand nodes, where the indicated weight n is the exponent of WN . If no weight is indicated, it is assumed to be unity (or equivalent to WN0 ). Thus the output of the step shown in Fig. 1 is c = a + WNn b. a a + W n N n b b Figure 1: Signal-flow graph notation. With this notation the combining of the two length N/2 DFTs is illustrated in Fig. 2 for N=8. Each right-hand node is one of the Fm , formed by the appropriate combination of Am and Bm as developed in Eqs. (3), (4) and (5). If N is divisible by 4, the process may be repeated, and each length N/2 DFT may be formed by decimating the two N/2 sequences into even and odd components, forming the length N/4 DFTs, 2 A 4 - p o in t D F T 0 fr o m e v e n s a m p le s { f 0 , f 2 , f 4 , f6 } 1 A 2 0 A B 3 { f 1 , f 3 , f 5 , f7 } B F4 = A 4 0 5 1 F5 = A 6 B B 2 7 3 F6 = A F7 = A + W 5 + W 6 7 + W B 2 8 = A 5 B 7 = A 4 B 6 8 3 B 5 8 2 B 4 8 1 B 3 8 + W 4 0 1 8 + W 3 B 8 + W 2 F3 = A 0 8 + W 1 F2 = A 2 3 + W 0 F1 = A 1 A 4 - p o in t D F T fr o m o d d s a m p le s F0 = A = A 6 B 7 = A - W 0 - W 1 - W 2 3 - W 0 8 8 1 B 8 B 0 2 1 8 B 3 B 2 3 Figure 2: Signal-flow graph representation of combining two length-4 DFT sequences, derived from even and odd sample sets, into a length-8 sequence as defined in Eq. (3). and combining these back into a length N/2 DFT, as is shown for N = 8 in Fig. 3. Notice that all weights in the figure are expressed by convention as exponents of W8 . In general, if the length of the data sequence is an integer power of 2, that is N = 2q for integer q, the DFT sequence {Fm } may be formed by adding additional columns to the left and halving the length of the DFT at each step, until the length is two. For example if N = 256 = 28 a total of seven column operations would be required. The final step is to evaluate the N/2 length-2 DFTs. Each one may be written F0 = f0 + W20 f1 = f0 + f1 F1 = f0 + W21 f1 = f0 − f1 , which is simply the sum and difference of the two sample points. No complex multiplications are necessary. The 2-point DFT is shown in signal-flow graph form in Fig. 4, and is known as the FFT butterfly. The complete FFT algorithm for N = 8 is formed by adding a column of 2-point DFTs to the left of Fig. 3, as shown in Fig. 5. We note that if N = 2q , there will be q = log2 (N ) columns in the signal-flow graph, and after the sum and difference to form the 2-point DFTs there will be log2 (N ) − 1 column operations, each involving N/2 complex multiplications, giving a total of N/2 (log2 (N ) − 1) ≤ N 2 . We will address the issue of computational savings in more detail later. Input Bit-Reversal: Notice that the algorithm of Fig. 5 requires that the input sequence {fn } be re-ordered in the left-hand column to accomplish the even-odd decomposition at each step. For this reason this form of the FFT is known as the FFT with decimation in time through input bit reversal. The term input bit reversal refers to a simple algorithm to determine the position k of the sample fn in the re-ordered sequence: 1. Express the index n as a N -bit binary number. 2. Reverse the order of the binary digits (bits) in the binary number. 3. Translate the bit-reversed number back into decimal, to create the position in the sequence k. 3 4 2 - p o in t D F T s 2 4 - p o in t D F T s 2 - p o in t D F T fr o m s a m p le s 8 - p o in t D F T 0 { f 0 , f4 } 2 { f 2 , f6 } F 1 F2 2 6 F 3 3 2 - p o in t D F T fr o m s a m p le s F 4 0 { f 1 , f5 } F 5 6 4 2 - p o in t D F T fr o m s a m p le s 4 5 2 { f 3 , f7 } 0 1 4 2 - p o in t D F T fr o m s a m p le s F 0 6 7 F F 6 7 Figure 3: Two steps to combining four length-2 DFT sequences into a length-8 sequence. f0 f F0 = f0 + W 0 F 1 = f0 + W 1 1 f1 = f0 + f1 0 2 2 1 f 1 = f0 - f1 Figure 4: The 2-point DFT butterfly. 4 2 - p o in t D F T s 8 - p o in t d a ta s e t f0 f 2 4 - p o in t D F T s 0 0 f 2 f 6 2 f1 f 2 6 4 F 3 F F 4 5 5 6 4 0 7 F2 4 0 4 f3 1 3 0 f5 F 2 6 4 0 1 4 0 F 0 4 4 8 - p o in t D F T 7 F F 6 7 Figure 5: Complete radix-2 with decimation in time (input bit reversal) FFT algorithm for N = 8. 4 For example, the re-ordered position of f37 in a data sequence of length N = 256 = 28 is found from bit reversal 3710 = 001001012 −→ 101001002 = 16410 so that f37 would be positioned at k = 164 in the decimated input sequence. Table 1 shows the complete re-ordering process for N = 8. Input position n: Bit reversal Modified position k: 0 (000)2 ↓ (000)2 0 1 (001)2 ↓ (100)2 4 2 (010)2 ↓ (010)2 2 3 (011)2 ↓ (110)2 6 4 (100)2 ↓ (001)2 1 5 (101)2 ↓ (101)2 5 6 (110)2 ↓ (011)2 3 7 (111)2 ↓ (111)2 7 Table 1: Determination of input bit-reversed sequence order for N = 8. A sample FFT routine with input bit-reversal, written as a MATLAB script fftx(), is presented in the Appendix. This is not intended as a substitute for MATLAB’s built-in FFT functions, and is intended only as a tutorial example. As noted in the introduction, the algorithmic structure presented here is just one of several FFT structures. Some of these take the time domain data in the correct order and produce the DFT values in bit-reversed order. Other algorithms are based on frequency decomposition, and matrix factoring operations. The Inverse Fast Fourier Transform (IFFT): The inverse FFT is defined as fn = −1 2πmn 1 N� F m ej N , N m=0 n = 0, . . . , N − 1 (6) While the IFFT can be implemented in the same manner as the FFT described above, it is possible to use a forward FFT routine to compute the IFFT as follows: Since the conjugate of a product is the product of the conjugates, if we take the complex conjugate of both sides we have fn = −1 2πmn 1 N� F m e−j N . N m=0 The right-hand side is recognized as the DFT of F m and can be computed using a forward FFT, such as described above. The complete IDFT may therefore be computed by conjugating the output, that is � � −1 2πmn 1 N� fn = F m e−j N , N m=0 n = 0, . . . , N − 1 The steps are: 1. Conjugate the data set {Fm }. 2. Compute the forward FFT. 3. Conjugate the result and divide by N . A tutorial MATLAB inverse FFT routine, ifftx(), is presented in the Appendix. 5 (7) Computational Savings of the FFT: As expressed above the computational requirements (in terms of complex multiplications) is MFFT = (N/2) log2 (N ) if the initial 2-point DFTs are implemented with exponentials. The number of complex multiplications for the direct DFT computation is MDFT = N 2 ) We can therefore define a speed improvement factor MFFT /MDFT as is shown in Table 2. It should be realized,however, that in practice these dramatic improvements shown in Table 2 will not be realized because of the significant improvement in arithmetic processing speeds, particularly with regard to multiplications. N 4 8 16 32 64 128 256 512 1024 2048 4096 MDFT 16 64 256 1,024 4,096 16,384 65,536 262,144 1,048,576 4,194,304 16,777,216 MFFT 4 12 32 80 192 448 1024 2,304 5,120 11,264 24,576 MFFT /MDFT 0.25 0.188 0.125 0.0781 0.0469 0.0273 0.0156 0.00879 0.00488 0.00268 0.00146 Table 2: Computational speed improvement from the use of the FFT . Efficient Use of the FFT: The computational efficiency of an FFT function may be increased by several methods, depending on the application: • If the algorithm is to be applied to a fixed length data set repetitively, the table of complex exponentials WNn (sines and cosines) may be pre-computed and stored. • If it is known that {fn } is real, then {Fm } is conjugate symmetric about its mid-point, so that only N/2 values have to be computed. • It is possible to compute the FFT of a two data sequences of length N in one length N FFT by creating a complex sequence {fn } with one data set in the real part and the other in the imaginary part. Let {cn } and {dn } be two real data sets of length N , and form them into a single complex sequence fn = cn + jdn 0 ≤ n ≤ N − 1. Then since the DFT is a linear operation Fm = Cm + jDm . The two sequences can be expressed in terms of fn as follows: � 1� cn = fn + f n 2 � 1 � dn = fn − f n 2j and therefore � �� 1� Cm = DFT {fn } + DFT f n 2 � �� 1 � Dm = DFT {fn } − DFT f n 2j 6 � � but the DFT of f n = {FN −m } so that the two DFTs may be recovered from {Fm }: Cm = Dm = � 1� Fm + F N −m 2 � 1 � Fm − F N −m 2j Because this method requires real data, it will not generally be applicable to computing an inverse FFT. • Similarly, it is possible to compute the DFT of a length 2N data sequence by splitting into two length N sequences and placing each in the real and imaginary parts of a length N complex sequence. Perform an even/odd decomposition into two N -point sequences: cn = f2n even samples, dn = f2n+1 odd samples. Form the complex N -point sequence {fn } = {cn } + j {dn }, and use the method described above, so that Cm = Dm = � 1� Fm + F N −m 2 � 1 � Fm − F N −m . 2j Finally, combine the two N -point DFTs into a single 2N -point DFT using Eqs. (3), Fm = Cm + W2kN Dm m = 0, 1, . . . , N − 1 W2kN Dm m = 0, 1, . . . , N − 1 Fm+N = Cm − Again, this method will generally not be applicable to inverse FFT calculations because the Fm are usually complex. MATLAB FFT Routines: MATLAB contains several DFT routines based on the FFT, including multidimensional FFTs. These routines are not restricted to radix-2 lengths (and will even handle prime lengths) by adjusting the internal algorithm. In particular, the following four routines are useful for one dimensional transforms: fft(x) is the discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is applied to each column. For N-D arrays, the FFT operation operates on the first nonsingleton dimension. f(x,n) is the N-point FFT, padded with zeros if x has less than N points and truncated if it has more. fft(X,[],DIM) or fft(X,N,DIM) applies the FFT operation across the dimension DIM. ifft(x) is the inverse discrete Fourier transform of X. ifft(X,N) is the N -point inverse transform. ifft(X,[],DIM) or ifft(X,N,DIM) is the inverse discrete Fourier transform of X across the dimension DIM. fftshift(X) is useful for visualizing the Fourier transform with the zero-frequency component in the middle of the spectrum. The routine shifts the zero-frequency component to center of 7 spectrum. For vectors, fftshift(X) swaps the left and right halves of X. For matrices, fftshift(X) swaps the first and third quadrants and the second and fourth quadrants. For N -dimensional arrays, fftshiftT(X) swaps ”half-spaces” of X along each dimension. fftshift(X,DIM) applies the fftshift operation along the dimension DIM. ifftshift(X) undoes the effects of fftshift. For vectors, ifftshift(X) swaps the left and right halves of X. For matrices, ifftshift(X) swaps the first and third quadrants and the second and fourth quadrants. For N -dimensional arrays, ifftshift(X) swaps ”half-spaces” of X along each dimension. ifftshift(X,DIM) applies the ifftshift operation along the dimension DIM. See MATLaB’s help and doc facilities for more information and descriptions of multidimensional transforms. References: [1 ] Cooley, J. W. and J. W. Tukey, ”An Algorithm for the Machine Computation of the Complex Fourier Series,”Mathematics of Computation, Vol. 19, April 1965, pp. 297-301. 8 Appendix: MATLAB Demonstration FFT Functions: %--------------------------------------------------------------­ % *** 2.161 Signal Processing - Continuous and Discrete *** % % fftx - Tutorial FFT routine to demonstrate the Radix-2 FFT with % decimation in time. % % Fout = fftx(f, k) % Arguments: % f - time domain data set ()real or complex % k - size of the data set k = log_2(N) % % Author: D. Rowell % Revision: 1.0 10-1-2007 % %--------------------------------------------------------------­ function Fout = fftx(f,ksize) N = 2^ksize; % Move the input data to the output array Fout=f; % % Perform the "bit-reversed" re-ordering of the input data % MR = 0; for M = 1:N-1 L = N/2; while MR + L > N-1 L = L/2; end MR = mod(MR,L) + L; if MR >= M temp = Fout(M+1); %swap the data points Fout(M+1) = Fout(MR+1); Fout(MR+1) = temp; end end % % Now perform the column operations % L = 1; while L < N ISTEP = 2*L; for K = 1:L W = exp(-i*pi*(K-1)/L); for P = K:ISTEP:N Q = P + L; % P & Q are the two points to be updated WBm = W*Fout(Q); Fout(Q) = Fout(P) - WBm; % Q identifies the "odd" block. Fout(P) = Fout(P) + WBm; % P identifies "even" block, end end L = ISTEP; end 9 %--------------------------------------------------------------­ % *** 2.161 Signal Processing - Continuous and Discrete *** % % ifftx - Tutorial inverse FFT routine to demonstrate the Radix-2 % inverseFFT using a forward FFT routine. % % fout = ifftx(F, k) Arguments: % % F - frequency domain data set % k - size of the data set k = log_2(N) % % Author: D. Rowell % Revision: 1.0 10-1-2007 % %--------------------------------------------------------------­ function fout = ifftx(F,ksize) N = 2^ksize; % % Conjugate the input data fout=conj(F); % Take the forward FFT, conjugate and divide by N fout =demofft(fout,ksize); fout=conj(fout)/N; The following MATLAB commands take the FFT and then the IFFT of a 16-point data set and plot the results: % Form a 16-point data seqience t =[0:.375:15*.375]; f = exp(-t).*sin(t); % Compute the FFT and then the IFFT F = fftx(f,4); f1 = real(ifftx(F,4)); % Plot the results plot(t, f, t, f1, ’o’); Demonstration of 2.161 Sample FFT Functions 0.35 Input data Output data 0.3 f(t) = exp(−t)sin(t) 0.25 0.2 0.15 0.1 0.05 0 −0.05 0 1 2 3 4 5 6 time (s) 10