Introduction to the Fast-Fourier Transform (FFT) Algorithm

1 Introduction to the Fast-FourierTransform (FFT) RamalingamDepartment of electrical EngineeringIIT Ramalingam (EE Dept., IIT Madras)Intro to FFT1 / 30 The Discrete Fourier Transform (DFT)DFT of anN-point sequencexn,n= 0,1,2,..,N 1 isdefined asXk=N 1 n=0xne j2 kNnk= 0,1,2, ,N 1 AnN-point sequence yields anN-point transformXkcan be expressed as aninner product:Xk=[1e j2 kNe j2 j2 kN(N 1)] 1 Ramalingam (EE Dept., IIT Madras)Intro to FFT2 / 30 The Discrete Fourier Transform (DFT)Notation:WN=e j2 N. Hence,Xk=[ (N 1)kN] 1 By varyingkfrom 0 toN 1 and combining theNinnerproducts, we get the following:X=WxWis anN Nmatrix, called as the DFT Matrix Ramalingam (EE Dept.)

2 , IIT Madras)Intro to FFT3 / 30 The DFT MatrixW= 111 11 WNW2N WN 1N1W2NW4N W2(N 1) 1NW2(N 1)N W(N 1)(N 1)N N NThe notationWNis used if we want to make the size of theDFT matrix Ramalingam (EE Dept., IIT Madras)Intro to FFT4 / 30 How Many Complex Multiplications Are Required?Each inner product requiresNcomplex multiplicationsThere areNinner productsHence we requireN2multiplicationsHowever, the first row and first column are all 1s, andshouldnot be counted as multiplicationsThere are 2N 1 such instancesHence, the number of complex multiplications isN2 2N+ 1, , (N 1) Ramalingam (EE Dept., IIT Madras)Intro to FFT5 / 30 How Many Complex Additions Are Required?

3 Each inner product requiresN 1 complex additionsThere areNinner productsHence we requireN(N 1) complex Ramalingam (EE Dept., IIT Madras)Intro to FFT6 / 30 Total Operation CountNo. of complex multiplications: (N 1)2No. of complex additions:N(N 1)The operation count for multiplications and additions assumesthatWkNhas been computed offline and is available in memoryIf pre-computed values ofWkNare not available, then theoperation count will increaseWe will assume that all the requiredWkNhave beenpre-computed and are Ramalingam (EE Dept., IIT Madras)Intro to FFT7 / 30 Operation Count Makes DFT ImpracticalFor largeN,(N 1)2 N2N(N 1) N2 Hence both multiplications and additions areO(N2)IfN= 103, thenO(N2) = 106, , a million!

4 This makes the straightforward method slow and impracticaleven for a moderately long Ramalingam (EE Dept., IIT Madras)Intro to FFT8 / 30 The Divide and Conquer ApproachSupposeNis even and we split the sequence into two sequence has N/2 pointsSuppose we compute theN2point DFT of each sequenceMultiplications : 2 (N2)2=N22 Suppose we are able to combine the individual DFT results toget the originally required DFTSome computational overhead will be consumed to combinethe two resultsIfN22+ overhead<N2, then this approach will reduce theoperation Ramalingam (EE Dept., IIT Madras)Intro to FFT9 / 30 The Divide and Conquer ApproachLetN= 8 Straightforward implementation requires,approximately, 64multiplicationsThe divide and conquer approach requires,approximately,2 (82)2+ overhead, , 32 + overhead multiplicationsQuestions:Can the two DFTs be combined to get the original DFT ?

5 If so, how ? What is the overhead involved ?Will 32 + overhead be less than 64 ? Ramalingam (EE Dept., IIT Madras)Intro to FFT10 / 30 The Decimation in Time (DIT) AlgorithmFrom{xn}form two sequences as follows:{gn}={x2n} {hn}={x2n+1}{gn}contains the even-indexed samples, while{hn}containsthe odd-indexed samplesThe DFT of{xn}isXk=N 1 n=0xnWnkN=N2 1 r=0x2rW(2r)kN+N2 1 r=0x2r+1W(2r+1)kN=N2 1 r=0grW(2r)kN+WkNN2 1 r=0hrW(2r) Ramalingam (EE Dept., IIT Madras)Intro to FFT11 / 30 The Decimation in Time (DIT) AlgorithmBut,W2rkN=e j2 N(2rk)=e j2 N/2(rk)=WrkN/2and henceXk=N2 1 r=0grWrkN/2+WkNN2 1 r=0hrWrkN/2=Gk+WkNHkk= 0,1,..,N 1{Gk}and{Hk}areN2point DFTsThe overhead for combining the twoN2point DFTs is themultiplicative factorWkNfork= 0,1.

6 ,N 1 WkNis called twiddle factor Ramalingam (EE Dept., IIT Madras)Intro to FFT12 / 30 The Decimation in Time (DIT) AlgorithmTheN/2 point DFTs{Gk}and{Hk}are periodic with periodN/2Gk+N2=GkHk+N2=HkWk+N2N= WkNHence, ifXk=Gk+WkNHk, thenXk+N2=Gk WkNHkWkNHkneeds to be computed only once fork= 0 toN2 1 Thus, the multiplication overhead due to the twiddle factors Ramalingam (EE Dept., IIT Madras)Intro to FFT13 / 30 Butterfly DiagramWkN WkNGkHkXkXk+N2Xk=Gk+WkNHkXk+N2=Gk+N2+Wk+ N2 NHk+N2=Gk Ramalingam (EE Dept., IIT Madras)Intro to FFT14 / 30 The Decimation in Time (DIT) AlgorithmFigure Flowgraph of Decimation in Time Algorithm forN= 8 (Oppenheim and Schafer,Discrete-Time SignalProcessing, 3rd edition, Pearson Education, 2010, p.)

7 726) Ramalingam (EE Dept., IIT Madras)Intro to FFT15 / 30 Divide and Conquer Results in Savings!ForN= 8, the straightforward approach requires,approximately, 64 multiplicationsThe Divide and Conquer approach, after the first stage,requires 32 + 4 = 36 multiplicationsThus, this approach clearly reduces the number of additionsand multiplications Ramalingam (EE Dept., IIT Madras)Intro to FFT16 / 30 Reusing the Divide and Conquer StrategyThe same idea can be applied for calculating theN2pointDFT of the sequences{gr}and{hr}Computational savings can be obtained by dividing{gr}and{hr}intotheirodd- and even-indexed halvesThis idea can be applied recursively log2 Ntimes ifNis apower of 2 Such algorithms are called radix 2 algorithmsIfN= 2 , then the final stage sequences are all of length 2 For a 2-point sequence{p0,p1}, the DFT coefficients areP0=p0+p1P1=p0 Ramalingam (EE Dept.

8 , IIT Madras)Intro to FFT17 / 30 DIT Flowgraph forN= 8 Figure Flowgraph of Decimation in Time Algorithm forN= 8 (Oppenheim and Schafer,Discrete-Time SignalProcessing, 3rd edition, Pearson Education, 2010, p. 730) Ramalingam (EE Dept., IIT Madras)Intro to FFT18 / 30 Overall Operation CountThe direct method requiresN2multiplicationsAfter the first split,N2 2(N2)2+N2N2is due to thetwiddle factorsAfter the second split,(N2)2 2(N4)2+N4 Hence,N2 2(N2)2+N2 first stage 4(N4)2+N2+N2 second stageGeneralizing, if there are log2 Nstages, the number ofmultiplications needed will be,approximately, Ramalingam (EE Dept., IIT Madras)Intro to FFT19 / 30 Overall Operation CountIfWk+N2N= WkNis not considered, the overhead count willbeNand notN2In this case,N2 2(N2)2+N first stage 4(N4)2+N+N second stageHence the overall multiplication count will beNlog2 NForN= 1024N2= 1,048,576 Nlog2N= 10,240 Savings of two orders of magnitude!

9 Ramalingam (EE Dept., IIT Madras)Intro to FFT20 / 30 Input Sequence OrderRecall that, forN= 8, the first split requires the data to bearranged as follows:x0,x2,x4,x6,x1,x3,x5,x7In the second and final split, the data appear in the followingorder:x0,x4,x2,x6,x1,x5,x3,x7 The final order is said to be in bit reversed form:OriginalBinary FormReversed Ramalingam (EE Dept., IIT Madras)Intro to FFT21 / 30An Algorithm For Sequence ReversalConsider the card sequence7, 8, 9, 10, J, Q, K, AFirst, reverse pairwise:8, 7, 10, 9, Q, J, A, KThen swap the adjacent pairs:10, 9, 8, 7, A, K, Q, JFinally, swap the two groups of 4 (each group is half theoriginal size):A, K, Q, J, 10, 9, 8, 7 Done!

10 Ramalingam (EE Dept., IIT Madras)Intro to FFT22 / 30 How To Use It For Bit ReversalThe first step of swapping of bits pairwise can be done withbitwise AND/OR and bit shift operatorsPick out the even and odd bits by using masksABCDEFGH & 01010101 = 0B0D0F0 HABCDEFGH & 10101010 = A0C0E0G0 Left shift the first result and right shift the second resultB0D0F0H00A0C0E0 GBitwise OR the above resultsB0D0F0H0 0A0C0E0G = BADCFEHGP airwise bit swapping accomplished! Ramalingam (EE Dept., IIT Madras)Intro to FFT23 / 30C Code For Bit Reversalunsigned reverse_bits(unsigned input){//works on 32-bit machineinput = (input & 0x55555555) << 1 | (input & 0xAAAAAAAA) >> 1;input = (input & 0x33333333) << 2 | (input & 0xCCCCCCCC) >> 2;input = (input & 0x0F0F0F0F) << 4 | (input & 0xF0F0F0F0) >> 4;input = (input & 0x00FF00FF) << 8 | (input & 0xFF00FF00) >> 8;input = (input & 0x0000 FFFF) << 16 | (input & 0xFFFF0000) >> 16;return input;}Bit reversal for the entire array can take a large overhead ifperformed inefficientlyThere are several efficient algorithms for sorting an array inbit-reversed orderBit reversal on uniprocessorsby Alan H.