理解快速傅里葉變換（FFT）算法
快速傅里葉變換（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：
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函數，我們來對結果進行仔細檢查，
x = np.random.random(1024) np.allclose(DFT_slow(x), np.fft.fft(x))
輸出：
True
現在為了驗證我們的算法有多慢，對比下兩者的執行時間
%timeit DFT_slow(x) %timeit np.fft.fft(x)
輸出：
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
但我們不是到這步就停下來，只要我們的小傅里葉變換是偶倍數，就可以再作分治，直到分解出來的子問題小到無法通過分治提高效率，接近極限時，這個遞歸是 O(n logn) 級的。
這個遞歸算法能在python里快速實現，當子問題被分解到合適大小時，再用回原本那種“慢方法”。
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])</pre>
現在我們做個快速的檢查，看結果是否正確：
x = np.random.random(1024)
np.allclose(FFT(x), np.fft.fft(x))
True
然后與“慢方法”的運行時間對比下：
%timeit DFT_slow(x)
%timeit FFT(x)
%timeit np.fft.fft(x)
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。
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()</pre>
x = np.random.random(1024)
np.allclose(FFT_vectorized(x), np.fft.fft(x))
True
因為我們的算法效率更大幅地提升了，所以來做個更大的測試（不包括DFT_slow）
x = np.random.random(1024 * 16)
%timeit FFT(x)
%timeit FFT_vectorized(x)
%timeit np.fft.fft(x)
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函數在實際中并不適用。作為數據科學家，我們可以暫且引進能夠包含基本應用工具的黑盒子，是由我們有更多算法思想的同事參入所構建，但是我堅定的相信：當我們對所應用的數據的底層算法有更深的理解，我們也將會成為更優秀的從業者。
原文鏈接： jakevdp 翻譯： 伯樂在線 - Forever Shining
譯文鏈接： http://blog.jobbole.com/58246/
本文由用戶 jopen 自行上傳分享，僅供網友學習交流。所有權歸原作者，若您的權利被侵害，請聯系管理員。
轉載本站原創文章，請注明出處，并保留原始鏈接、圖片水印。
本站是一個以用戶分享為主的開源技術平臺，歡迎各類分享！