- 編注：這篇譯文由@unblock 和@jingliang 共同完成。
- 再次推薦：《如果看了此文你還不懂傅立葉變換，那就過來掐死我吧【完整版】》
快速傅立葉變換（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：
|
1 2 3 4 5 6 7 8 9 |
import numpy as np def 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函式，我們來對結果進行仔細檢查，
|
1 2 |
x = np.random.random(1024) np.allclose(DFT_slow(x), np.fft.fft(x)) |
輸出：
|
1 |
True |
現在為了驗證我們的演算法有多慢，對比下兩者的執行時間
|
1 2 3 |
%timeit DFT_slow(x) %timeit np.fft.fft(x) |
輸出：
|
1 2 |
10 loops, best of 3: 75.4 ms per loop 10000 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裡快速實現，當子問題被分解到合適大小時，再用回原本那種“慢方法”。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
def FFT(x): """A recursive implementation of the 1D Cooley-Tukey FFT""" x = np.asarray(x, dtype=float) N = x.shape[0] if N % 2 > 0: raise ValueError("size of x must be a power of 2") elif N <= 32: # this cutoff should be optimized return DFT_slow(x) else: X_even = FFT(x[::2]) X_odd = FFT(x[1::2]) factor = np.exp(-2j * np.pi * np.arange(N) / N) return np.concatenate([X_even + factor[:N / 2] * X_odd, X_even + factor[N / 2:] * X_odd]) |
現在我們做個快速的檢查，看結果是否正確：
|
1 2 |
x = np.random.random(1024) np.allclose(FFT(x), np.fft.fft(x)) |
|
1 |
True |
然後與“慢方法”的執行時間對比下：
|
1 2 3 4 5 |
%timeit DFT_slow(x) %timeit FFT(x) %timeit np.fft.fft(x) |
|
1 2 3 |
10 loops, best of 3: 77.6 ms per loop 100 loops, best of 3: 4.07 ms per loop 10000 loops, best of 3: 24.7 µs per loop |
現在的演算法比之前的快了一個數量級。而且，我們的遞迴演算法漸近於 O(n logn) 。我們實現了FFT 。
需要注意的是，我們還沒做到numpy的內建FFT演算法，這是意料之中的。numpy 的 fft 背後的FFTPACK演算法 是以 Fortran 實現的，經過了多年的調優。此外，我們的NumPy的解決方案，同時涉及的Python堆疊遞迴和許多臨時陣列的分配，這顯著地增加了計算時間。
還想加快速度的話，一個好的方法是使用Python/ NumPy的工作時，儘可能將重複計算向量化。我們是可以做到的，在計算過程中消除遞迴，使我們的python FFT更有效率。
向量化的NumPy
注意上面的遞迴FFT實現，在最底層的遞迴，我們做了N/32次的矩陣向量乘積。我們的演算法會得益於將這些矩陣向量乘積化為一次性計算的矩陣-矩陣乘積。在每一層的遞迴，重複的計算也可以被向量化。因為NumPy很擅長這類操作，我們可以利用這一點來實現向量化的FFT。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
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() |
|
1 2 |
x = np.random.random(1024) np.allclose(FFT_vectorized(x), np.fft.fft(x)) |
|
1 |
True |
因為我們的演算法效率更大幅地提升了，所以來做個更大的測試（不包括DFT_slow）
|
1 2 3 4 5 6 |
x = np.random.random(1024 * 16) %timeit FFT(x) %timeit FFT_vectorized(x) %timeit np.fft.fft(x) |
|
1 2 3 |
10 loops, best of 3: 72.8 ms per loop 100 loops, best of 3: 4.11 ms per loop 1000 loops, best of 3: 505 µs per loop |
我們的實現又提升了一個級別。這裡我們是以 FFTPACK中大約10以內的因數基準，用了僅僅幾十行 Python + NumPy程式碼。雖然沒有相應的計算來證明， Python版本是遠優於 FFTPACK源，這個你可以從這裡瀏覽到。
那麼 FFTPACK是怎麼獲得這個最後一點的加速的呢？也許它只是一個詳細的記錄簿， FFTPACK花了大量時間來保證任何的子計算能夠被複用。我們這裡的numpy版本涉及到額外的記憶體的分配和複製，對於如Fortran的一些低階語言就能夠很容易的控制和最小化記憶體的使用。並且Cooley-Tukey演算法還能夠使其分成超過兩部分（正如我們這裡用到的Cooley-Tukey FFT基2演算法），而且，其它更為先進的FFT演算法或許也可以能夠得到應用，包括基於卷積的從根本上不同的方法（例如Bluestein的演算法和Rader的演算法）。結合以上的思路延伸和方法，就可使陣列大小即使不滿足2的冪，FFT也能快速執行。
我希望他們提供大量的基於FFT資料分析是怎樣進行的背景，儘管純Python函式在實際中並不適用。作為資料科學家，我們可以暫且引進能夠包含基本應用工具的黑盒子，是由我們有更多演算法思想的同事參入所構建，但是我堅定的相信：當我們對所應用的資料的底層演算法有更深的理解，我們也將會成為更優秀的從業者。