15-451 Algorithms 04/01 * Fast Fourier Transform (FFT) (review session today 7pm here) =========================================================================== Today: continue on with Fast Fourier Transform. Review where we got to last time. How to do inverse. Fast implementation. Connection to "usual" Fourier transform. Last time: Have a polynomial A(x) = a_0 + a_1*x + a_2*x^2 + ... + a_{m-1}*x^{m-1} Want to evaluate it at m points of our choosing in time O(m log m). Idea: let w be a principal mth root of unity: w = sin(2pi/m) + i cos(2pi/m) Evaluate A at 1, w, w^2, w^3, ..., w^{m-1}. Notice that w^{m/2} = -1, so w^{m/2+j} = -w^j, which means we can divide into even and odd like we discussed last time. A(x) = A_even(x^2) + x*A_odd(x^2) A(-x) = A_even(x^2) - x*A_odd(x^2) Recursively, we'll be evaluating A_even and A_odd on 1, w^2, (w^2)^2, (w^2)^3,...,(w^2)^{m/2-1}. Since w^2 is the principal m/2 root of unity, we have the same problem but smaller inductively, so the divide-and-conquer recursion works. Final algorithm runs in time O(m log m). Parallel FFT Circuit ==================== FFT can also be done efficiently in hardware on a parallel circuit. We're not focusing on parallel computation in this course, but in case of FFT, I think drawing the parallel cirsuit helps to give a feel for what's going on. --> show picture from book. Go through example too. A(x) = a_0 + a_1x + ... + a_7x^7. A_e(x) = a_0 + a_2x + a_4x^2 + a_6x^3. Goal: eval at 1, w^2, w^4, w^6. Equivalently: 1, w^2, -1, -w^2 A_ee(x) = a_0 + a_4(x). Eval at 1 and w^4 = -1. Recursively get a_0, a_4 back. Output V_ee = (a_0 + 1*a_4, a_0 - 1*a_4) A_eo(x) = a_2 + a_6(x). Eval at 1,-1. Output V_eo = (a_2 + a_6, a_2 - a_6) Output V_e = (V_ee[0] + 1*V_eo[0], V_ee[1] + w^2*V_eo[1], V_ee[0] - 1*V_eo[0], V_ee[1] - w^2*V_eo[1]) --> Note: this is called a butterfly network. Modular FFT =========== One thing about FFT is we need to store w as a floating point number and do floating point calculations. So, technically we should be worrying about roundoff error and numerical stability, which all works out fine but it's painful. We got into this mess because we needed some w such that w^m = 1 but w^j != 1 for 0 < j < m. A neat alternative is to use modular arithmetic. For instance, say we just want to find the coefficients modulo 17. The cool thing about this is now we can just use "2" as an 8th root of unity. In particular, the powers of 2 are 1,2,4,8,16,15,13,9 (equivalently: 1,2,4,8,-1,-2,-4,-8). Then 2^8=1. A plug: Iain will talk more about this in the optional rec session this evening, in context of going through examples of FFT and also because next week we'll be doing number-theoretic algorithms. So, if you're weak on modular arithmetic, come by the recitation this evening. The DFT and its inverse ======================= This is the Discrete Fourier Transform. Given a vector A = [a_0, a_1, ..., a_m-1] we compute F(A) = [y_0, y_1, ..., y_m-1] where y_j = a_0 + a_1*w^j + a_2*w^{2j} + ... + a_{m-1}*w^{(m-1)j} Now, one way to think of this is that we're implicitly doing a matrix-vector product: +---------------------------------------+ +-----+ +-----+ | 1 1 1 1 ... 1 | | a_0 | | y_0 | | | | | | | | 1 w w^2 w^3 ... w^{m-1} | | a_1 | | y_1 | | | | | | | | 1 w^2 w^4 w^6 ... w^{2(m-1)} | | a_2 | | y_2 | | | | | = | | | 1 w^3 w^6 w^9 ... w^{3(m-1)} | | a_3 | | y_3 | | | | | | | | ..... ..... | | ... | | ... | | | | | | | | 1 w^{-1} w^{-2} w^{-3} ... w | |a_m-1| |y_m-1| +---------------------------------------+ +-----+ +-----+ (Note w^{m-1} = w^{-1} since w^m = 1) I say "implicitly" because we don't even have time to write down the matrix. But, this view is helpful in thinking about how to do the next thing we need, which is the INVERSE transform F^{-1}. To do the inverse transform, what we want is to multiply by the inverse of the F matrix. As it turns out, this inverse looks very much like F itself. In particular, - notice that w^{-1} is also a principal mth root of unity. Let's define \bar{F} to be the fourier transform matrix using w^{-1} instead of w. Then, Claim: F^{-1} = (1/m) * \bar{F} Proof: What is the i,j entry of \bar{F} * F? It is: 1 + w^{j-i} + w^{2j-2i} + w^{3j-3i} + ... + w^{(m-1)j - (m-1)i} If i=j, then these are all = 1, so the sum is m. Then when we divide by m we get 1. If i!=j, then the claim is these all cancel out and we get zero. Maybe easier to see if we let z = w^{j-i}, so then the sum is: 1 + z + z^2 + z^3 + z^4 + ... + z^{m-1}. Then can see these cancel by picture. E.g, do m=8 and j-i=3, then j-i=2. Or can use the formula for summations. Since F^{-1} has the same form as F, we can also compute it in time O(n log n) with the same algorithms, just feeding it w^{-1} = cos(2pi/m) - isin(2pi/m) instead of w = cos(2pi/m) + isin(2pi/m). -- break-- CONNECTION TO USUAL FOURIER TRANSFORM ------------------------------------- We've built up the FFT without saying in what sense this is a "Fourier transform". The usual Fourier transform works like this: You're given some function f -- so the meaning of f(x) is the value of function f at point x. The fourier transform is a new function F, where F(z) is the correlation of f with the sine wave of frequency z. F(z) = Integral of f(x)e^{-2*pi*i*zx}dx In the discrete version, say we have sampled f at points j=0,1,2,3,...,m-1. So, we are represnting f as a vector of length m. If we think of the a_i's as representing these values, then the FFT is computing the correlations with these different frequencies. So this is sort of backwards compared to how we've been thinking about it. In fact, some definitions call the DFT what we've been calling the *inverse* DFT. In other words, m-1 m-1 F(z) = (1/m) SUM f(j) e^{-(2*pi*i/m)jz} = (1/m)SUM f(j) w^{-jz} j=0 j=0 This makes more sense in a way since to get the correlation we should be dividing by m anyway.