快速入門 FFT
本文适合不想完全了解 FFT 原理不過想實作 FFT 算法的人
FFT 簡介
FFT 也就是常說的快速傅裡葉變換，其實本質上和傅裡葉變換沒有什麼不同，也就是通過算法減少了一些乘法和加法的運算；DFT 的時間複雜度是 N² 級别的，而 FFT 的卻是對數級别的複雜度，進而在 FFT 在資料量大的時候運算時間并不會高太多，這就使得 FFT 具有了在工程上的使用價值
FFT 算法簡介
上面這張圖是來自維基的介紹，簡單來說的話，FFT 将資料按照索引位址按奇數和偶數分為兩半，這兩半資料還可以繼續抽取，直到最後分解為單獨的資料；然後利用了機關根的一些特性，減少了一些乘加步驟
但是順序輸入的資料，最後得到的結果是亂序的；為了得到順序的輸出，一般我們會将輸入資料提前進行整理；因為 FFT 的輸出索引位址會位反轉，可以利用這個特性提前整理好輸入資料，例如：輸入 [1 2 3 4 5 6 7 8] 就整理為 [1 5 3 7 2 6 4 8]
FFT 第一個版本實作
最開始我為了能弄懂 FFT ，我就通過手算來了解 FFT 的過程，也就是通過上面那張圖中的公式，一級級計算 FFT 的結果
最後根據手算的步驟實作了下面這個版本，其中機關根的複指數計算就直接手動計算完存表表了，在實作後發現隻能夠适用 4 個數的輸入
import cmath
table = [1+0j,2+0j,3+0j,4+0j]
rotate_table = [1+0j,0-1j,-1+0j,0+1j]
result = []
N = 4
k = 0
r = 0
tmp = 0
while k < len(table):
while r < len(table)/2:
tmp = tmp + table[2*r]*rotate_table[(2*r*k)%4] + table[2*r+1]*rotate_table[((2*r+1)*k)%4]
r = r + 1
result.append(tmp)
tmp = 0
k = k + 1
r = 0
print(result)
運作結果如下
PS D:\Files\dsp\python> python.exe .\my_fft.py
[(10+0j), (-2+2j), (-2+0j), (-2-2j)]
matlab fft 的結果如下
>> x = [1 2 3 4];X = fft(x)
X =
10 + 0i -2 + 2i -2 + 0i -2 - 2i
FFT 第二個版本的實作
最開始看到書上的 FFT 特點的時候，還不是太了解，是以上面的那個實作非常的簡陋
原址運算的特點使得 FFT 占用空間為 N，每個蝶形運算結束後，輸出節點的兩個值直接放到原輸入的兩個節點中；這個特點給了我對算法的基本思路，算法是從整理輸入資料順序開始的
不滿足 2 的 N 次方的情況直接不考慮了
通過觀察，可以看出每層蝶形運算的規律：如果将每層的蝶形運算分組的話，每層的蝶形運算的組數和每組的蝶形運算次數的乘積是一定的，為 N/2。具體可以參考這篇文章
基于這些特點就有了下面這個第二個版本的 FFT
import cmath
import math
def reverseBits(n, bits): # 整理資料子函數
return int(bin(n)[2:].rjust(int(bits), '0')[::-1], base=2)
table = [1,2,3,4,5,6,7,8] # 輸入資料
result = [] # 存儲結果
L = math.log(len(table), 2) # 蝶形運算的層數
N = len(table) # 輸入資料的長度
tmp = 0
tmp_i = 0
while tmp_i < N/2: # 整理資料過程
tmp = table[tmp_i]
table[tmp_i] = table[reverseBits(int(tmp_i), int(L))]
table[reverseBits(tmp_i, L)] = tmp
tmp_i = tmp_i + 1
print(table)
k = 0 # 一組裡的 k 次蝶形
m = 1
i = len(table)/2 # i 組蝶形
j = 0 # 第 j 組蝶形
while m <= L:
while j < i: # i組蝶形 4/2=2
while k < N/2/i: # k次蝶形運算
r = (j*(2**m) + k)*(2**(L-m))
tmp_1 = table[j*(2**m) + k]
tmp_2 = table[j*(2**m) + k+2**(m-1)] * cmath.exp(-1j*2*cmath.pi*r/N)
table[j*(2**m) + k] = tmp_1 + tmp_2
table[j*(2**m) + k+2**(m-1)] = tmp_1 - tmp_2
k = k + 1
j = j + 1
k = 0
i = i/2
j = 0
m = m+1
print(table)
運作結果如下：
PS D:\Files\dsp\python> python.exe .\my_fft_2.py
[1, 5, 3, 7, 2, 6, 4, 8]
[(36+1.3471114790620885e-14j), (-4.000000000000009+9.656854249492383j), (-4.000000000000006+3.9999999999999982j), (-4.000000000000005+1.656854249492377j), (-4-1.0042103753008295e-14j), (-3.9999999999999942-1.6568542494923832j), (-3.999999999999994-4.000000000000002j), (-3.999999999999991-9.656854249492376j)]
matlab 中的結果如下：
>> x = [1 2 3 4 5 6 7 8];X = fft(x)
X =
Columns 1 through 4:
36.0000 + 0.0000i -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i
Columns 5 through 8:
-4.0000 + 0.0000i -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i
和 matlab 中的對比還存在精度的問題，這個與 python 的 cmath 庫提供的指數運算有關，我覺得可以考慮像第一版的實作一樣，将指數運算優化掉，不過暫時不打算不去考慮這個問題了