快速傅里葉變換（Fast Fourier Transform）是信號處理與數據分析領域里最重要的算法之一。沒有正規計算機科學課程背景的我，使用這個算法多年，但這周我卻突然想起自己從沒思考過為什么FFT能如此快速地計算離散傅里葉變換。我打開一本老舊的算法書，欣賞了JW Cooley 和 John Tukey 在1965年的文章中，以看似簡單的計算技巧來講解這個東西。
本文的目標是，深入Cooley-Tukey FFT 算法，解釋作為其根源的“對稱性”，并以一些直觀的python代碼將其理論轉變為實際。我希望這次研究能使數據科學家（例如我），對這個算法的背景原理有更全面的認識。
FFT（快速傅里葉變換）本身就是離散傅里葉變換（Discrete Fourier Transform）的快速算法，使算法復雜度由原本的O(N^2) 變為 O(NlogN)，離散傅里葉變換DFT，如同更為人熟悉的連續傅里葉變換，有如下的正、逆定義形式：
xn 到 Xk 的轉化就是空域到頻域的轉換，這個轉換有助于研究信號的功率譜，和使某些問題的計算更有效率。事實上，你還可以查看一下我們即將推出的天文學和統計學的圖書的第十章（這里有一些圖示和python代碼）。作為一個例子，你可以查看下我的文章《用python求解薛定諤方程》，是如何利用FFT將原本復雜的微分方程簡化。
正因為FFT在那么多領域里如此有用，python提供了很多標準工具和封裝來計算它。NumPy 和 Scipy 都有經過充分測試的封裝好的FFT庫，分別位于子模塊 numpy.fft 和 scipy.fftpack 。我所知的最快的FFT是在 FFTW包中 ，而你也可以在python的pyFFTW 包中使用它。
雖然說了這么遠，但還是暫時先將這些庫放一邊，考慮一下怎樣使用原始的python從頭開始計算FFT。
簡單起見，我們只關心正變換，因為逆變換也只是以很相似的方式就能做到。看一下上面的DFT表達式，它只是一個直觀的線性運算：向量x的矩陣乘法，
矩陣M可以表示為
這么想的話，我們可以簡單地利用矩陣乘法計算DFT：
Python| 123456789 | import numpy as npdef DFT_slow(x): """Compute the discrete Fourier Transform of the 1D array x""" x = np.asarray(x, dtype=float) N = x.shape[0] n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) |
對比numpy的內置FFT函數，我們來對結果進行仔細檢查，
Python| 12 | x=np.random.random(1024)np.allclose(DFT_slow(x),np.fft.fft(x)) |
輸出：
| 1 | True |
現在為了驗證我們的算法有多慢，對比下兩者的執行時間
Python| 123 | %timeitDFT_slow(x) %timeitnp.fft.fft(x) |
輸出：
| 12 | 10 loops, best of 3: 75.4 ms per loop10000 loops, best of 3: 25.5 μs per loop |
使用這種簡化的實現方法，正如所料，我們慢了一千多倍。但問題不是這么簡單。對于長度為N的輸入矢量，FFT是O(N logN)級的，而我們的慢算法是O(N^2)級的。這就意味著，FFT用50毫秒能干完的活，對于我們的慢算法來說，要差不多20小時！ 那么FFT是怎么提速完事的呢？答案就在于他利用了對稱性。
算法設計者所掌握的最重要手段之一，就是利用問題的對稱性。如果你能清晰地展示問題的某一部分與另一部分相關，那么你就只需計算子結果一次，從而節省了計算成本。
Cooley 和 Tukey 正是使用這種方法導出FFT。 首先我們來看下
的值。根據上面的表達式，推出：
對于所有的整數n，exp[2π i n]=1。
最后一行展示了DFT很好的對稱性：
簡單地拓展一下：
同理對于所有整數 i 。正如下面即將看到的，這個對稱性能被利用于更快地計算DFT。
DFT 到 FFT：
利用對稱性 Cooley 和 Tukey 證明了，DFT的計算可以分為兩部分。從DFT的定義得：
我們將單個DFT分成了看起來相似兩個更小的DFT。一個奇，一個偶。目前為止，我們還沒有節省計算開銷，每一部分都包含(N/2)?N的計算量，總的來說，就是N^2 。
技巧就是對每一部分利用對稱性。因為 k 的范圍是 0≤k<N ， 而 n 的范圍是 0≤n<M≡N/2 ， 從上面的對稱性特點來看，我們只需對每個子問題作一半的計算量。O(N^2) 變成了 O(M^2) 。
但我們不是到這步就停下來，只要我們的小傅里葉變換是偶倍數，就可以再作分治，直到分解出來的子問題小到無法通過分治提高效率，接近極限時，這個遞歸是 O(n logn) 級的。
這個遞歸算法能在python里快速實現，當子問題被分解到合適大小時，再用回原本那種“慢方法”。
Python| 123456789101112131415 | defFFT(x): """A recursive implementation of the 1D Cooley-Tukey FFT""" x=np.asarray(x,dtype=float) N=x.shape[0] ifN%2>0: raiseValueError("size of x must be a power of 2") elifN<=32: # this cutoff should be optimized returnDFT_slow(x) else: X_even=FFT(x[::2]) X_odd=FFT(x[1::2]) factor=np.exp(-2j*np.pi*np.arange(N)/N) returnnp.concatenate([X_even+factor[:N/2]*X_odd, X_even+factor[N/2:]*X_odd]) |
現在我們做個快速的檢查，看結果是否正確：
Python| 12 | x = np.random.random(1024)np.allclose(FFT(x), np.fft.fft(x)) |
| 1 | True |
然后與“慢方法”的運行時間對比下：
Python| 12345 | %timeit DFT_slow(x) %timeit FFT(x) %timeit np.fft.fft(x) |
| 123 | 10loops,bestof3:77.6msperloop100loops,bestof3:4.07msperloop10000loops,bestof3:24.7μsperloop |
現在的算法比之前的快了一個數量級。而且，我們的遞歸算法漸近于 O(n logn) 。我們實現了FFT 。
需要注意的是，我們還沒做到numpy的內置FFT算法，這是意料之中的。numpy 的 fft 背后的FFTPACK算法 是以 Fortran 實現的，經過了多年的調優。此外，我們的NumPy的解決方案，同時涉及的Python堆棧遞歸和許多臨時數組的分配，這顯著地增加了計算時間。
還想加快速度的話，一個好的方法是使用Python/ NumPy的工作時，盡可能將重復計算向量化。我們是可以做到的，在計算過程中消除遞歸，使我們的python FFT更有效率。
注意上面的遞歸FFT實現，在最底層的遞歸，我們做了N/32次的矩陣向量乘積。我們的算法會得益于將這些矩陣向量乘積化為一次性計算的矩陣-矩陣乘積。在每一層的遞歸，重復的計算也可以被向量化。因為NumPy很擅長這類操作，我們可以利用這一點來實現向量化的FFT。
Python| 12345678910111213141516171819202122232425262728 | def FFT_vectorized(x): """A vectorized, non-recursive version of the Cooley-Tukey FFT""" x = np.asarray(x, dtype=float) N = x.shape[0] if np.log2(N) % 1 > 0: raise ValueError("size of x must be a power of 2") # N_min here is equivalent to the stopping condition above, # and should be a power of 2 N_min = min(N, 32) # Perform an O[N^2] DFT on all length-N_min sub-PRoblems at once n = np.arange(N_min) k = n[:, None] M = np.exp(-2j * np.pi * n * k / N_min) X = np.dot(M, x.reshape((N_min, -1))) # build-up each level of the recursive calculation all at once while X.shape[0] < N: X_even = X[:, :X.shape[1] / 2] X_odd = X[:, X.shape[1] / 2:] factor = np.exp(-1j * np.pi * np.arange(X.shape[0]) / X.shape[0])[:, None] X = np.vstack([X_even + factor * X_odd, X_even - factor * X_odd]) return X.ravel() |
| 12 | x=np.random.random(1024)np.allclose(FFT_vectorized(x),np.fft.fft(x)) |
| 1 | True |
因為我們的算法效率更大幅地提升了，所以來做個更大的測試（不包括DFT_slow）
Python| 123456 | x=np.random.random(1024*16)%timeitFFT(x) %timeitFFT_vectorized(x) %timeitnp.fft.fft(x) |
| 123 | 10loops,bestof3:72.8msperloop100loops,bestof3:4.11msperloop1000loops,bestof3:505μsperloop |
我們的實現又提升了一個級別。這里我們是以 FFTPACK中大約10以內的因數基準，用了僅僅幾十行 Python + NumPy代碼。雖然沒有相應的計算來證明， Python版本是遠優于 FFTPACK源，這個你可以從這里瀏覽到。
那么 FFTPACK是怎么獲得這個最后一點的加速的呢？也許它只是一個詳細的記錄簿， FFTPACK花了大量時間來保證任何的子計算能夠被復用。我們這里的numpy版本涉及到額外的內存的分配和復制，對于如Fortran的一些低級語言就能夠很容易的控制和最小化內存的使用。并且Cooley-Tukey算法還能夠使其分成超過兩部分（正如我們這里用到的Cooley-Tukey FFT基2算法），而且，其它更為先進的FFT算法或許也可以能夠得到應用，包括基于卷積的從根本上不同的方法（例如Bluestein的算法和Rader的算法）。結合以上的思路延伸和方法，就可使陣列大小即使不滿足2的冪，FFT也能快速執行。
我希望他們提供大量的基于FFT數據分析是怎樣進行的背景，盡管純Python函數在實際中并不適用。作為數據科學家，我們可以暫且引進能夠包含基本應用工具的黑盒子，是由我們有更多算法思想的同事參入所構建，但是我堅定的相信：當我們對所應用的數據的底層算法有更深的理解，我們也將會成為更優秀的從業者。
新聞熱點
疑難解答
