演算法導論FFT程式碼

《演算法導論》是我目前為止看到的最好的有關於寫FFT的書，雖然和訊號處理一類的樹不太相同，演算法導論對於FFT的出發點更加純粹一點，那就是如果快速計算DFT，而並不深究DFT的物理意義，切入點也是多項式乘法。說實話，這個演算法很難，並不是很直觀，理解困難的很大原因是演算法構造的如此精妙，以至於讓人迷失了究竟演算法的目的在何處，當然，現在很多庫整合的FFT高效且穩定，不過深入學習一下FFT的原理我認為還是有必要的。

詳細的原理就不說了，直接上python實現的程式碼

python由於支援了複數運算，所以和書中的虛擬碼基本差不多，無須做更多修改，計算結果和numpy中自帶的FFT函式完全一致，當然，此演算法只能應用於輸入時2的冪次的時候。

其中第一個函式是遞迴版本的FFT，第二個則是迭代版本的FFT。

def recursiveFFT(a):
    n = len(a)
    if n == 1:
        return a
    wm = np.exp((2*np.pi*1j)/n) #旋轉因子
    #wm = np.cos(2*np.pi/n)+1j*np.sin(2*np.pi/n) # 尤拉展開後的旋轉因子，與上式相同
    w = 1
    a0 = a[[x for x in range(0,n) if x%2 == 0]] # 提取出下標為偶數的係數
    a1 = a[[x for x in range(0,n) if x%2 == 1]] # 提取出下標為奇數的係數
    y0 = recursiveFFT(a0) #　遞迴呼叫
    y1 = recursiveFFT(a1)
    y = np.empty(n,dtype=complex) #注意這裡型別必須宣告為complex型，否則運算時會自動捨去虛部
    for k in range(0,n/2):  # 歸一
        #t = w*y1[k]
        y[k] = y0[k]+(w*y1[k])
        y[k+n/2] = y0[k] - w*y1[k]
        w = w*wm
    return y

def rev(x,num):
    n = int(np.log2(num)) # 求出最大的二進位制位數
    ans = []
    while x!= 0:
        ans.append(x%2)
        x /= 2
    while len(ans)!= n:
        ans.append(0)   # 空位補零
    reval = 0
    for index,i in enumerate(ans):
        reval += np.power(2,n-index-1)*i  #二進位制轉十進位制
    return reval

def bitReserveCopy(a):
    n = len(a)
    newA = np.empty(n,dtype=complex) # 此處必須宣告型別為complex，否則下面的FFT運算會自動捨去虛部
    for k in range(0,n):
        newA[rev(k,n)] = a[k]  # 二進位制逆序排列
    return newA

def MyFFT(a):
    A = bitReserveCopy(a)
    n = len(a)
    for s in range(1,int(np.log2(n))+1): # 遍歷樹的所有層
        m = np.power(2,s) # m很關鍵，代表了接下來n次單位根中的n，也表示接下來遍歷新序列時的步長
        wm = np.exp((2*np.pi*1j)/m)
        for k in range(0,n,m):
            w = 1
            for j in range(0,m/2):
                t = w*A[k+j+m/2] # 蝶形運算
                u = A[k+j]
                A[k+j] = u+t
                A[k+j+m/2] = u-t
                w = w*wm
    return A


x = np.random.rand(8)
print x
my1 = recursiveFFT(x)
print my1
my2 = MyFFT(x)
print my2
y = np.fft.fft(x)
print y