快速傅里葉變換（FFT）算法【詳解】

快速傅里葉變換（Fast Fourier Transform）是信號處理與數據分析領域裏最重要的算法之一。我打開一本老舊的算法書，欣賞了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 import numpy as np
2 def DFT_slow(x):
3     """Compute the discrete Fourier Transform of the 1D array x"""
4     x = np.asarray(x, dtype=float)
5     N = x.shape[0]
6     n = np.arange(N)
7     k = n.reshape((N, 1))
8     M = np.exp(-2j * np.pi * k * n / N)
9     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<N ， 而 n 的範圍是 0≤n<M≡N/2 ， 從上面的對稱性特色來看，咱們只需對每一個子問題做一半的計算量。O(N^2) 變成了 O(M^2) 。

但咱們不是到這步就停下來，只要咱們的小傅里葉變換是偶倍數，就能夠再做分治，直到分解出來的子問題小到沒法經過分治提升效率，接近極限時，這個遞歸是 O(n logn) 級的。

這個遞歸算法能在python裏快速實現，當子問題被分解到合適大小時，再用回本來那種「慢方法」。

 1 def FFT(x):
 2     """A recursive implementation of the 1D Cooley-Tukey FFT"""
 3     x = np.asarray(x, dtype=float)
 4     N = x.shape[0]
 5  
 6     if N % 2 > 0:
 7         raise ValueError("size of x must be a power of 2")
 8     elif N <= 32:  # this cutoff should be optimized
 9         return DFT_slow(x)
10     else:
11         X_even = FFT(x[::2])
12         X_odd = FFT(x[1::2])
13         factor = np.exp(-2j * np.pi * np.arange(N) / N)
14         return np.concatenate([X_even + factor[:N / 2] * X_odd,
15                                X_even + factor[N / 2:] * X_odd])

如今咱們作個快速的檢查，看結果是否正確：

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

 1 def FFT_vectorized(x):
 2     """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
 3     x = np.asarray(x, dtype=float)
 4     N = x.shape[0]
 5  
 6     if np.log2(N) % 1 > 0:
 7         raise ValueError("size of x must be a power of 2")
 8  
 9     # N_min here is equivalent to the stopping condition above,
10     # and should be a power of 2
11     N_min = min(N, 32)
12  
13     # Perform an O[N^2] DFT on all length-N_min sub-problems at once
14     n = np.arange(N_min)
15     k = n[:, None]
16     M = np.exp(-2j * np.pi * n * k / N_min)
17     X = np.dot(M, x.reshape((N_min, -1)))
18  
19     # build-up each level of the recursive calculation all at once
20     while X.shape[0] < N:
21         X_even = X[:, :X.shape[1] / 2]
22         X_odd = X[:, X.shape[1] / 2:]
23         factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
24                         / X.shape[0])[:, None]
25         X = np.vstack([X_even + factor * X_odd,
26                        X_even - factor * X_odd])
27  
28     return X.ravel()

 

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也能快速執行。