Fast Fourier Transform Jean Baptiste Joseph Fourier 1768
Fast Fourier Transform Jean Baptiste Joseph Fourier (1768 -1830) These lecture slides are adapted from CLRS. Princeton University • COS 423 • Theory of Algorithms • Spring 2002 • Kevin Wayne
Fast Fourier Transform Applications. n n n Perhaps single algorithmic discovery that has had the greatest practical impact in history. Optics, acoustics, quantum physics, telecommunications, systems theory, signal processing, speech recognition, data compression. Progress in these areas limited by lack of fast algorithms. History. n Cooley-Tukey (1965) revolutionized all of these areas. n Danielson-Lanczos (1942) efficient algorithm. n Runge-König (1924) laid theoretical groundwork. n Gauss (1805, 1866) describes similar algorithm. n Importance not realized until advent of digital computers. 2
Polynomials: Coefficient Representation Degree n polynomial. Addition: O(n) ops. Evaluation: O(n) using Horner's method. Multiplication (convolution): O(n 2). 3
Polynomials: Point-Value Representation Degree n polynomial. n Uniquely specified by knowing p(x) at n different values of x. y = p(x) yj xj x 4
Polynomials: Point-Value Representation Degree n polynomial. Addition: O(n). Multiplication: O(n), but need 2 n points. Evaluation: O(n 2) using Lagrange's formula. 5
Best of Both Worlds Can we get "fast" multiplication and evaluation? Representation Multiplication Evaluation coefficient O(n 2) O(n) point-value O(n) O(n 2) FFT O(n log n) ! Yes! Convert back and forth between two representations. coefficient multiplication O(n 2) evaluation FFT interpolation inverse FFT O(n log n) point-value multiplication O(n) 6
Converting Between Representations: Naïve Solution Evaluation (coefficient to point-value). n n Given a polynomial p(x) =a 0 + a 2 x 1 +. . . + an-1 xn-1, choose n distinct points {x 0, x 1, . . . , xn-1 } and compute yk = p(xk), for each k using Horner's method. O(n 2). Interpolation (point-value to coefficient). n n n Given n distinct points {x 0, x 1, . . . , xn-1 } and yk = p(xk), compute the coefficients {a 0, a 1, . . . , an-1 } by solving the following linear system of equations. Note Vandermonde matrix is invertible iff xk are distinct. O(n 3). 7
Fast Interpolation: Key Idea Key idea: choose {x 0, x 1, . . . , xn-1 } to make computation easier! n Set xk = xj? 8
Fast Interpolation: Key Idea Key idea: choose {x 0, x 1, . . . , xn-1 } to make computation easier! n n Set xk = xj? Use negative numbers: set xk = -xj so that (xk )2 = (xj )2. – set xk = -xn/2 + k E = peven(x 2) = a 0 + a 2172 + a 4174 + a 6176 +. . . + an-217 n-2 – O = x podd(x 2) = a 117 + a 3173 + a 5175 + a 7177 +. . . + an-117 n - 1 – y 0 = E + O, yn/2 = E - O – 9
Fast Interpolation: Key Idea Key idea: choose {x 0, x 1, . . . , xn-1 } to make computation easier! n n n Set xk = xj? Use negative numbers: set xk = -xj so that (xk )2 = (xj )2. – set xk = -xn/2 + k Use complex numbers: set x k = k where is nth root of unity. – (xk )2 = (-xn/2 + k )2 – (xk )4 = (-xn/4 + k )4 – (xk )8 = (-xn/8 + k )8 10
Roots of Unity An nth root of unity is a complex number z such that zn = 1. n = e 2 i / n = principal n th root of unity. n e i t = cos t + i sin t. n i 2 = -1. n There are exactly n roots of unity: k, k = 0, 1, . . . , n-1. 2 = i 3 1 4 = -1 0 = 1 7 5 6 = -i 11
Roots of Unity: Properties L 1: Let be the principal nth root of unity. If n > 0, then n/2 = -1. n Proof: = e 2 i / n n/2 = e i = -1. (Euler's formula) L 2: Let n > 0 be even, and let and be the principal nth and (n/2)th roots of unity. Then ( k ) 2 = k. n Proof: ( k ) 2 = e (2 k)2 i / n = e (k) 2 i / (n / 2) = k. L 3: Let n > 0 be even. Then, the squares of the n complex n th roots of unity are the n/2 complex (n/2) th roots of unity. n Proof: If we square all of the nth roots of unity, then each (n/2) th root is obtained exactly twice since: – L 1 k + n / 2 = - k – thus, ( k + n / 2) 2 = ( k ) 2 – L 2 both of these = k – k + n / 2 and k have the same square 12
Divide-and-Conquer Given degree n polynomial p(x) = a 0 + a 1 x 1 + a 2 x 2 +. . . + an-1 xn-1. n n n Assume n is a power of 2, and let be the principal n th root of unity. Define even and odd polynomials: – peven(x) : = a 0 + a 2 x 1 + a 4 x 2 + a 6 x 3 +. . . + an-2 xn/2 - 1 – podd(x) : = a 1 + a 3 x 1 + a 5 x 2 + a 7 x 3 +. . . + an-1 xn/2 - 1 – p(x) = peven(x 2) + x podd(x 2) Reduces problem of – evaluating degree n polynomial p(x) at 0, 1, . . . , n-1 to – n evaluating two degree n/2 polynomials at: ( 0)2, ( 1)2, . . . , ( n-1)2. L 3 peven(x) and podd(x) only evaluated at n/2 complex (n/2) th roots of unity. 13
FFT Algorithm FFT (n, a 0, a 1, a 2, . . . , an-1) if (n == 1) return a 0 // n is a power of 2 e 2 i / n (e 0, e 1, e 2, . . . , en/2 -1) FFT(n/2, a 0, a 2, a 4, . . . , an-2) (d 0, d 1, d 2, . . . , dn/2 -1) FFT(n/2, a 1, a 3, a 5, . . . , an-1) for k = 0 to n/2 - 1 y k e k + k d k yk+n/2 ek - k dk O(n) complex multiplies if we pre-compute k. return (y 0, y 1, y 2, . . . , yn-1) 14
Recursion Tree a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 a 0 , a 2 , a 4 , a 6 a 0 , a 4 a 0 a 4 a 1 , a 3 , a 5 , a 7 a 2 , a 6 a 2 a 6 a 3 , a 7 a 1 , a 5 a 1 a 5 a 3 a 7 "bit-reversed" order 15
Proof of Correctness Proof of correctness. Need to show yk = p( k ) for each k = 0, . . . , n-1, where is the principal n th root of unity. n n Base case. n = 1. Algorithm returns y 0 = a 0 0. Induction step. Assume algorithm correct for n / 2. – let be the principal (n/2) th root of unity – ek = peven( k ) = peven( 2 k ) by Lemma 2 – dk = podd( k ) = podd( 2 k ) by Lemma 2 – recall p(x) = peven(x 2) + x podd(x 2) 16
Best of Both Worlds Can we get "fast" multiplication and evaluation? Representation Multiplication Evaluation coefficient O(n 2) O(n) point-value O(n) O(n 2) FFT O(n log n) ! Yes! Convert back and forth between two representations. coefficient multiplication O(n 2) evaluation FFT interpolation inverse FFT O(n log n) point-value multiplication O(n) 17
Inverse FFT Forward FFT: given {a 0, a 1, . . . , an-1 } , compute {y 0, y 1, . . . , yn-1 }. Inverse FFT: given {y 0, y 1, . . . , yn-1 } compute {a 0, a 1, . . . , an-1 }. 18
Inverse FFT Great news: same algorithm as FFT, except use -1 as "principal" nth root of unity (and divide by n). 19
Inverse FFT: Proof of Correctness Summation lemma. Let be a primitive n th root of unity. Then n n If k is a multiple of n then k = 1. Each nth root of unity k is a root of xn - 1 = (x - 1) (1 + x 2 +. . . + xn-1), if k 1 we have: 1 + k(2) +. . . + k(n-1) = 0. Claim: Fn and Fn-1 are inverses. 20
Inverse FFT: Algorithm IFFT (n, a 0, a 1, a 2, . . . , an-1) if (n == 1) return a 0 // n is a power of 2 e-2 i/n (e 0, e 1, e 2, . . . , en/2 -1) FFT(n/2, a 0, a 2, a 4, . . . , an-2) (d 0, d 1, d 2, . . . , dn/2 -1) FFT(n/2, a 1, a 3, a 5, . . . , an-1) for k = 0 to n/2 - 1 yk (ek + k dk ) / n yk+n/2 (ek - k dk ) / n return (y 0, y 1, y 2, . . . , yn-1) 21
Best of Both Worlds Can we get "fast" multiplication and evaluation? Representation Multiplication Evaluation coefficient O(n 2) O(n) point-value O(n) O(n 2) FFT O(n log n) ! Yes! Convert back and forth between two representations. coefficient multiplication O(n 2) evaluation FFT O(n log n) interpolation inverse FFT O(n log n) point-value multiplication O(n) 22
Integer Arithmetic Multiply two n-digit integers: a = an-1. . . a 1 a 0 and b = bn-1. . . b 1 b 0. n Form two degree n polynomials. n Note: a = p(10), b = q(10). n Compute product using FFT in O(n log n) steps. n Evaluate r(10) = a b. n Problem: O(n log n) complex arithmetic steps. Solution. n n Strassen (1968): carry out arithmetic to suitable precision. – T(n) = O(n T(log n)) T(n) = O(n log n (log n)1+ ) Schönhage-Strassen (1971): use modular arithmetic. – T(n) = O(n log log n) 23
Extra Slides Princeton University • COS 423 • Theory of Algorithms • Spring 2002 • Kevin Wayne
Polynomials: Representation and Addition Degree n polynomial. n Coefficient representation. n Point-value representation. Addition. n O(n) for either representation. 25
Polynomials: Evaluation and Multiplication (convolution). n n Brute force. O(n 2). Use 2 n points. O(n). Evaluation. n n Horner's method. O(n). Lagrange's formula. O(n 2). 26
- Slides: 26