6.046J Lecture 5: Fast Fourier Transform

Lecture 5 Fast Fourier Transform Supplemental reading in CLRS: Chapter 30 The algorithm in this lecture, known since the time of Gauss but popularized mainly by Cooley and Tukey in the 1960s, is an example of the divide-and-conquer paradigm. Actually, the main uses of the fast Fourier transform are much more ingenious than an ordinary divide-and-conquer strategy— there is genuinely novel mathematics happening in the background. Ultimately, the FFT will allow us to do n computations, each of which would take Ω(n) time individually, in a total of Θ(n lg n) time. 5.1 Multiplication To motivate the fast Fourier transform, let’s start with a very basic question: How can we efficiently multiply two large numbers or polynomials? As you probably learned in high school, one can use essentially the same method for both: 385 ×426 2310 770 +1540 164010 (3x2 + 8x + 5) × (4x2 + 2x + 6) 18x2 + 48x + 30 6x3 + 16x2 + 10x 4 12x + 32x3 + 20x2 12x4 + 38x3 + 54x2 + 58x + 30 Of these two, polynomials are actually the easier ones to work with. One reason for this is that multiplication of integers requires carrying, while multiplication of polynomials does not. To make a full analysis of this extra cost, we would have to study the details of how large integers can be stored and manipulated in memory, which is somewhat complicated.1 So instead, let’s just consider multiplication of polynomials. Suppose we are trying to multiply two polynomials p, q of degree at most n with complex coefficients. In the high-school multiplication algorithm (see Figure 5.1), each row of the diagram is 1 More ideas are required to implement efficient multiplication of n-bit integers. In a 1971 paper, Schönhage and ³ ´ ∗ Strassen exhibited integer multiplication in O(n lg n · lg lg n) time; in 2007, Fürer exhibited O n lg n · 2O(lg n) time. The iterated logarithm lg∗ n is the smallest k such that lg lg · · · lg n (k times) ≤ 1. It is not known whether integer multiplication can be achieved in O(n lg n) time, whereas by the end of this lecture we will achieve multiplication of degree-n polynomials in O(n lg n) time. n+1 ◦ ◦ ◦ ◦ ◦ ◦ ··· ◦ × ◦ ◦ ··· ◦ ◦ ◦ ◦ ◦ ··· ··· ··· ◦ ◦ n+1 ◦ ··· ◦ ◦ 2n + 1 Figure 5.1. The high-school polynomial multiplication algorithm. obtained by multiplying p by a monomial. Thus any single row would take Ω(n) time to compute individually. There are n + 1 rows, so the computation takes Ω(n2 ) time. The key to improving this time is to consider alternative ways of representing p and q. Fix some N > n. We’ll choose N later, but for now all that matters is that N = O(n). A polynomial of degree at most N − 1 is uniquely determined by its values at N points. So, instead of storing the coefficients of p and q, we could represent p and q as the lists ­ ® ­ ® p z = p(z0 ), . . . , p(z N −1 ) and q z p(z0 ), . . . , p(z N −1 ) for any distinct complex numbers z0 , . . . , z N −1 . In this representation, computing pq is very cheap— the list ­ ® (pq) z = p(z0 )q(z0 ), . . . , p(z N −1 )q(z N −1 ) can be computed in O(n) time. If N > deg pq, then pq is the unique polynomial whose FFT is (pq) z . Algorithm: M ULTIPLY(p, q) 1. Fix some N = O(deg p + deg q) such that N > deg p + deg q. 2. Compute the sample vectors p z = 〈 p(z0 ), . . . , p(z N −1 )〉 and q z = 〈 q(z0 ), . . . , q(z N −1 )〉. 3. Compute the unique polynomial r such that r z = p z · q z , where · denotes entrywise multiplication.∗ 4. Return r. ∗ For example, 〈1, 2, 3〉 · 〈4, 5, 6〉 = 〈4, 10, 18〉. Thus, the cost of finding the coefficients of pq is equal to the cost of converting back and forth between the coefficients representation and the sample-values representation (see Figure 5.2). We are free to choose values of z0 , . . . , z N −1 that make this conversion as efficient as possible. 5.1.1 Efficiently Computing the FFT −1 In the fast Fourier transform, we choose z0 , . . . , z N −1 to be the Nth roots of unity, 1, ω N , ω2N , . . . , ω N , N 2π i/N where ω N = e . We make this choice because roots of unity enjoy the useful property that ω2N = ω N/2 . Lec 5 – pg. 2 of 6 PN × PN FFT / CN × CN Entrywise mult. (fast) Polynomial mult. (slow)  PN o FFT−1  CN Figure 5.2. Commutative diagram showing the cost of multiplication on either side of a fast Fourier transform. As we will see, the fastest way to get from the top-left to the bottom-left is through the FFT. (We can choose N to be a power of 2.) This allows us to divide the task of computing an FFT of size N into computing two FFTs of size N/2, and then combining them in a certain way.2 ¡ ¢ Algorithm: FFT N, p(x) = a N −1 x N −1 + a N −2 x N −2 + · · · + a 1 x + a 0 1 B As part of the specification, we assume N is a power of 2 if N = 1 then return a 0 4 Let 2 3 peven (y) ← a 0 + a 2 y + · · · + a N −2 y N/2−1 , podd (y) ← a 1 + a 3 y + · · · + a N −1 y N/2−1 be the even and odd parts of p 5 B Thus ¡ ¢ ¡ ¢ p(x) = peven x2 + xpodd x2 (5.1) even peven ) ω N/2 ← FFT ¡(N/2, p ¢ odd odd 7 p ωN/2 ← FFT N/2, p 8 B That is, 6 D ³ ´E even even even N/2−1 peven = p (1), p ω , . . . , p ω ( ) N/2 ω N/2 N/2 D ³ ´E odd odd odd odd N/2−1 p ωN/2 = p (1), p (ω N/2 ) , . . . , p ω N/2 . 9 −1 B Because ω N/2 = ω2N , we can calculate the vector p ωN = p(1), p (ω N ) , . . . , p ω N N ­ ¡ ¢® very quickly using (5.1): ­ ® −1 10 ω ← 1, ω N , . . . , ω N B the left side is a bold omega N odd 11 return peven + ω · p , where · denotes entrywise multiplication ω N/2 ω N/2 2 I suppose the “conquer” stage is when we recursively compute the smaller FFTs (but of course, each of these smaller FFTs begins with its own “divide” stage, and so on). After the “conquer” stage, the answers to the smaller problems are combined into a solution to the original problem. Lec 5 – pg. 3 of 6 even Above, we compute p ωN by computing podd ω N/2 and p ω N/2 and combining them in Θ(N) time. Thus, the running time of an FFT of size N satisfies the recurrence T(N) = 2T(N/2) + Θ(N). This recurrence is solved in CLRS as part of the Master Theorem in §4.5. The solution is T(N) = Θ(N lg N). 5.1.2 Computing the Inverse FFT Somewhat surprisingly, the inverse FFT can be computed in almost exactly the same way as the FFT. In this section we will see the relation between the two transforms. If you don’t have a background in linear algebra, you can take the math on faith. Let P N be the vector space of polynomials of degree at most N − 1 with complex coefficients. Then the FFT is a bijective linear map P N → C N . If we use the ordered basis 1, x, . . . , x N −1 for P N and the standard basis for C N , then the matrix of the FFT is   1 1 1 ··· 1 1 ω −1  ³ ´ ω2N · · · ωN N   N ij .  A = ωN = .. .. .. .. ..  0≤ i, j ≤ N −1  . . . .  . 2 (N −1) 2(N −1) N −1 · · · ωN 1 ωN ωN j The jth column of A −1 contains the coefficients of the polynomial which, when evaluated on ω N , i gives 1 and, when evaluated on ω N (i 6= j), gives zero. This polynomial is Q Q i 6= j ³ i 6= j ¡ i x − ωN j ¢ i ωN − ωN ´. We will not show the details here, but after about a page of calculation you can find that A −1 = 1 N  1 1 1  . A= N   .. 1 1 ω− N .. . (N −1) 1 ω− N 1 2 ω− N .. . ··· ··· .. . 2(N −1) ω− N ···  1 (N −1)  ω−  N , ..  .  2 (N −1) ω− N where the bar indicates entrywise complex conjugation. You can certainly check that the i, j entry of 1 N A A is ( NX −1 NX −1 1, i = j `(i − j) i` 1 − j` 1 ωN · N ωN = N ωN = 0, i 6= j, `=0 `=0 and similarly for 1 N A A. Thus FFT−1 (v) = 1 N FFT ¡ ¢ v. (5.2) Again, this is surprising. A priori, we would consider P N and C N to be “different worlds”—like two different implementations of the same abstract data structure, the coefficients representation and Lec 5 – pg. 4 of 6 the sample-values representation of polynomials serve essentially the same role but have different procedures for performing operations. Yet, (5.2) tells us that inverting an FFT involves creating a new polynomial whose coefficients are the sample values of our original polynomial. This of course has everything to do with the fact that roots of unity are special; it would not have worked if we had −1 as our sample points. not chosen 1, ω N , . . . , ω N N ¡ ¢ Algorithm: I NVERSE -FFT N, v = 〈v0 , . . . , v N −1 〉 ­ ® 1 v ← v0 , . . . , v N −1 2 Let p v (x) be the polynomial v N −1 x N −1 + · · · + v1 x + v0 3 return the value of ¡ ¢ 1 N FFT N, p v With all the details in place, the M ULTIPLY algorithm looks like this: Algorithm: FFT-M ULTIPLY(p, q) 1. Let N be the smallest power of 2 such that N − 1 ≥ deg p + deg q. 2. Compute v = FFT(N, p) and w = FFT(N, q), and let u = v · w, where · denotes entrywise multiplication. 3. Compute and return I NVERSE -FFT(N, u). ¡ ¢ It runs in Θ(N lg N) = Θ (deg p + deg q) lg(deg p + deg q) time. 5.2 Convolution Let p(x) = a n−1 x n−1 + · · · + a 1 x + a 0 and q(x) = b n−1 x n−1 + · · · + b 1 x + b 0 . The coefficient of x k in pq is X a ` b k−` , `∈Z2n−1 where the subscripts are interpreted modulo 2n − 1. This sort of operation is known as convolution, and is written as ∗. Convolution can be applied quite generally: if f and g are any functions Z2n−1 → C, then one common definition of the convolution operation is ( f ∗ g)(k) = X 1 f (`)g(k − `). 2n − 1 `∈Z2n−1 For each i, computing ( f ∗ g)(i) requires linear time, but the fast Fourier transform allows us to compute ( f ∗ g)(i) for all i = 0, . . . , 2n − 2 in a total of Θ(n lg n) time. Convolution appears frequently, which is part of the reason that the FFT is useful. The above notion of convolution can easily be generalized to allow f and g to be functions from any group G to any ring in which |G | is a unit. There is also a continuous version of convolution which involves an integral. This continuous version is useful in signal compression. The idea is to let f : R → R be an arbitrary signal and let g be a smooth bump function which spikes up near the origin and is zero elsewhere. (See Figure 5.3.) Then the convolution of f and g is a new function f ∗ g defined by Z ∞ ( f ∗ g)(t) = f (τ)g(t − τ) dτ. −∞ Lec 5 – pg. 5 of 6 Noisy signal f Bump function g Convolution f ∗ g Figure 5.3. A complicated, noisy signal f , a bump function g, and the convolution f ∗ g. The value ( f ∗ g)(t) is essentially an average of the values of f at points close to t. Thus, the funtion f ∗ g inherits smoothness from g while still carrying most of the information from f . As a bonus, if f comes with some unwanted random noise, then g will have much less noise. For example, f may be the raw output of a recording device, and computing f ∗ g may be the first step towards encoding f into a compressed audio format such as MP3. The point is that smooth functions are much easier to describe concisely, so creating a smooth version of f is useful in (lossy) data compression. Lec 5 – pg. 6 of 6