Fast Fourier Transform: Theory And Algorithms

Transcription

Fast Fourier Transform:Theory and AlgorithmsLecture 8Vladimir Stojanović6.973 Communication System Design – Spring 2006Massachusetts Institute of TechnologyCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].

Discrete Fourier Transform – A review Definition {Xk} is periodic Since {Xk} is sampled, {xn} must also be periodicFrom a physical point of view, both arerepeated with period NRequires O(N2) operations6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].2

Fast Fourier Transform History Twiddle factor FFTs (non-coprime sub-lengths) 1805 Gauss Predates even Fourier’s work on transforms!1903 Runge1965 Cooley-Tukey1984 Duhamel-Vetterli (split-radix FFT)FFTs w/o twiddle factors (coprime sub-lengths) 1960 Good’s mapping application of Chinese Remainder Theorem 100 A.D.1976 Rader – prime length FFT1976 Winograd’s Fourier Transform (WFTA)1977 Kolba and Parks (Prime Factor Algorithm – PFA)6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].3

Divide and conquer Divide and conquer always has less computationsSuppose all Il sets have same numberof elements N1 so, N N1*N2, r N2Each inner-most sum takes N12 multiplicationsThe outer sum will need N2 multiplications per output pointN2*N for the whole sum (for all output points) Hence, total number of multiplications6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].4

Generalizations The inner-most sum has to represent a DFT Only possible if the subset (possibly permuted) All main classes of FFTs can be cast in the above formBoth sums have same periodicity (Good’s mapping) No permutations (i.e. twiddle factors)All the subsets have same number of elements m N/r Has the same periodicity as the initial sequence(m,r) 1 – i.e. m and r are coprimeIf not, then inner sum is one stap of radix-r FFTIf r 3, subsets with N/2, N/4 and N/4 elements Split-radix algorithm6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].5

The cost of mapping The goal for divide and conquerDifferent types balance mapping with subproblemcostE.g. in radix-2 subproblems are trivial (only sum and differences)Mapping requires twiddle factors (large number ofmultiplies)E.g. in prime-factor algorithm Subproblems are DFTs with coprime lengths (costly)Mapping trivial (no arithmetic operations)6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].6

FFTs with twiddle factors Reintroduced by Cooley-Tukey ‘65Start from general divide and conquerKeep periodicity compatible withperiodicity of the input sequenceUse decimationalmost N1 DFTs of size N26.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].7

Cooley-Tukey FFT contd.can be taken mod N2Yn1 ,k Yn1 ,k21. N1 DFTs of length N22. N multplications with twiddle factors3. N2 DFTs of length N16.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].8

Step 1: Evaluate N1 DFTs of length N2Step 2: N multiplications with twiddle factorsStep 3: Evaluate N2 DFTs of length N1Vector xi mapped to matrix xn1,n2 (N1xN2)Compute N1 DFTs of length N2 on each rowPoint-to-point multiply with twiddle factorsCompute N2 DFTs of length N1 on the columns6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].9

2-D view of Cooley-Tukey mapping N 15 (N1 3, N2 5)x0 x1.x13 x14x0 x3 x6 x9 x12x1 x4 x7 x10 x13x2 x5 x8 x11 x6x4DFT-3x12x11DFT-5x10Figure by MIT OpenCourseWare. Cannot exchange the order of DFTs Because of twiddle multiplyDifferent mapping for N1 5, N2 3 Not just transposex0x1x2x3x4x5 x10x6 x11x7 x12x8 x13x9 x14Figure by MIT OpenCourseWare.6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].10

Radix 2 and radix 4 algorithms Lengths as powers of 2 or 4 are most popularAssume N 2n N1 2, N2 2n-1 (divides input sequence into even andodd samples – decimation in time – DIT)“Butterfly”(sum or difference followed orpreceeded by a twiddle factor multiply) Xm and XN/2 m outputs of N/2 2-pt DFTs on outputs of2, N/2-pt DFTs weighted with twiddle factors6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].11

DIT radix-2 implementations Several different ways Reorder the input data Input samples for innerDFTs in subsequentlocationsResults in bit-reversedinput, in-order output DITSelectively computeDFTs on evens andodds Perform in-placecomputationOutput in bit-reversedorder (X3 in position six(011- 110))X0DFTDFTX1N 4X2{x2i}N 2DFTN 2X3X4DFTDFTX5N 4w18X6{x2i 1}w2X7DivisionDFT ofinto even and N / 2odd numberedsequencesN 28DFTw38N 2Multiplicationby twiddlefactorsX0X4X1X5X2X6X3X7DFT of2Figure by MIT OpenCourseWare.Which type is this implementation?6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].12

Decimation in frequency (DIF) radix-2 implementation If reverse the role of N1 and N2, get DIF X2k1 N1 N/2, N2 2N / 2-1n1 0X2k1 1 n kW N1/ 12(xn xN / 2 n ),11N / 2-1n kW N1/ 12 Wn1(xn -xN / 2 n ).11Nn1 0X0DFTDFTX0X1N 2N 4X2X2DFT{x2k}X4X3N 2X4DFTX5N 2X6DFTX7N 2w18w28X6DFTX1N 4X3{x2k 1}X5w3DFT of Multiplication2by twiddlefactors6.973 Communication System DesignX78DFT ofN/2Figure by MIT OpenCourseWare.Cite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].13

Duality DIT - DIFX0DFTDFTX1N 4X2{x2i}N 2DFTN 2X3X4DFTDFTX5N 4w18X6{x2i 1}w2X78w38DFT of MultiplicationDivisionby twiddleinto even and N / 2factorsodd numberedsequencesN 2DFTN 2DFT of2X0X0X4X1X1X2X0DFTDFTN 2DFTN 2X5X3X2X4DFTX6X5N 2X3X6DFTX7X7N 2N 4X2{ }X4x2kw1X68DFTw28X1N 4X3{ }X5x2k 1w3X78DFT of Multiplication2by twiddlefactorsDFT ofN/2. are.Figure by MIT OpenCourseW Which one is DIT (DIF)?How can we get one from another?6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].14

Complexity of radix-2 FFTs DFT of length N replaced by two length-N/2 At the cost of N complex multiplications (twiddle) Iterate the scheme log2N-1 times And N complex additions (2pt DFTs)Obtain trivial transforms (length 2) of the length-N/2 DFTsTwiddle multiplies (WNi) Complex multiply – 3 real mult 3 real addIf i is multiple of N/4, no arithmetic operation required (why?)4 butterflies (one general, 3 special cases)6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].15

Radix-4 N 4n, N1 4, N2 N/4 Cost of length-4 DFT 4 DFTs of length N/43N/4 twiddle multipliesN/4 DFTs of length 4No multiplicationOnly 16 real T4x12x15DFT4DFT4DFT4Reduces the number of stages to log4N X0X12X1X13X2X14X3X15Figure by MIT OpenCourseWare.Radix-8 can reduce number of operations even more6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].16

Mixed-radix and Split-radix Mixed-radix Diferent radices in differentstagesSplit-radix X0DFT8DFT2R2X14X1x8DFT8Different radices in the samestage x0Simultaneously on differentparts of the transformCan achieve lowest numberof adds and multiplies forlength 2n X13X3S-RX15Figure by MIT OpenCourseWare.6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].17

Split-radix (DIF SRFFT) Look at DIF radix-2 X2k1 don’t have twiddlesX2k1 N / 2-1n1 0X2k 1 1nkW N1 / 12(xn xN / 2 n ),11N / 2-1n1 0nkW N1/ 12 Wn1(xn -xN / 2 n ).11NX0DFTDFTX0X1N 2N 4X2X2DFT{x2k}X4X3N 2X4DFTX5N 2X6DFTX7N 2w18w28Even samples X2k in DIF shouldbe computed separately fromother samples With same algorithm (recursively) asthe original sequenceNo general rule for odd samples DFTX1N 4X3{x2k 1}X5w83DFT of Multiplication2by twiddlefactors X6Radix-4 is more efficient than radix-2Higher radices are inefficientX2k 1N / 2-1n1 0X4k 1 1N / 4-1n1 0nkW N1/ 14Wn1NX[(xn -xN / 2 n )11 j(xn1 N / 4-xn1 3N / 4)],X7DFT ofN/2nkW N1/ 12(xn xN / 2 n ),11X4k 3 1N / 4-1n1 0Wn1k1 3n1N / 4WNFigure by MIT OpenCourseWare.X[(xn xN / 2 n )11-j(xn1 N / 4-xn1 3N / re by MIT OpenCourseWare.Cite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].

Split-radix (DIT SRFFT) Dual to DIF SRFFT Considers separately subsets {x2i}, {x4i 1} and {x4i 3} Redundancy in Xk, Xk N/4, Xk N/2, Xk 3N/4 computation6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].19

FFTs without twiddle factors Divide and conquer requirements N-long DFT computed from DFTs with lengths that arefactors of N (allows the inner sum to be a DFT) Provided that subsets Il guarantee periodic xiWhen N factors into co-prime factors N N1*N2 Starting from any xi form subset with compatible periodicity (theperiodicity of the subset divides the periodicity of the set){xi N1n2 n2 1,., N 2 1} or {xi N2 n1 n1 1,., N1 1}Both subsets have only one common point xiAllows a rearrangement of the input (periodic) vector into amatrix with a periodicity in both dimensions (rows andcolumns), both periodicities being comatible with the initial one6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].20

Good’s mapping FFTs without twiddle factors all based on the samemapping Turns original transform into a set of small DFTs withcoprime lengths{xi N1n2 n2 1,., N 2 1} or {xi N2 n1 n1 1,., N1 1}equivalent to01 2 3 4 5 6 7 8 9 10 11 12 13 14Good's mapping0369 1258 11 14 210 13 147. are.Figure by MIT OpenCourseW This mapping is one-to-one if N1 and N2 are coprimeAll congruences modulo N1 obtained For a given congruence modulo N2 and vice versa6.973 Communication System DesignCite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology.Downloaded on [DD Month YYYY].21

Just another arrangement of CRT Chinese Remainder Theorem (CRT) If we know the residue of some number k modulotwo coprime numbers N1 and N2k Nk NIt is possible to reconstruct k N NLet k N k1 k N k20 1 2 3 4 5 6 7 8 9 10 11 12 13 14Then k N N N1t1k2 N 2t2 k1 N1 11N2221t1 N12 1 and t2 N 2N1 1t1 multiplicative inverse of N1 mod N2t2 multiplicative inverse of N2 mod N1t1, t2 always exist since N1, N2 coprime (N1,N2) 1Good's mappingWhat are t1, t2 for N1 3, N2 5?CRT mapping Reversing N1 and N2 203658 11 14 210 13 109 12476 12 3910 1 7 13 45 11 2 8 14. are.Figure by MIT OpenCourseWResults in transposed mapping6.973Cite as: Vladimir Stojanovic, course materials for 6.973 Communication System Design, Spring 2006.MIT OpenCourseWare (http://ocw.mit.edu/), Massachuse

Step 2: N multiplications with twiddle factors . DIT radix-2 implementations Several different ways Reorder the input data Input samples for inner DFTs in subsequent locations Results in bit-reversed input, in-order output DIT Selectively compute DFTs on evens and odds Perform in-place computation Output in bit-reversed order (X3 in position six (011- 110)) 6.973 Communication System Design .