基2FFT的算法推導(dǎo)及python仿真

FFT的算法推導(dǎo)主要用到旋轉(zhuǎn)因子的周期性、對稱性和可約性：

基2FFT的頻域抽取法，將x(n)按照n的自然順序劃分為前后兩個部分：

所以當(dāng)K為偶數(shù)時，前后兩部分相加。當(dāng)k為奇數(shù)時相減。將頻域X(K)劃分為奇偶兩個序列，N點(diǎn)DFT就被分解為兩個N/2點(diǎn)的DFT：

可以得到蝶形圖如下：

進(jìn)而可以得到基2FFT頻域抽取代碼的實(shí)現(xiàn)方法：

隨后是數(shù)據(jù)倒換，如下圖：

可以看到基2FFT頻域抽取后的輸出位置排序就是自然數(shù)二進(jìn)制碼按位倒讀的值。

根據(jù)推導(dǎo)結(jié)果我們編寫python實(shí)現(xiàn)代碼：

首先根據(jù)FFT的點(diǎn)數(shù)計算需要迭代的次數(shù)，根據(jù)迭代次數(shù)例化一個loop_num+1*N的數(shù)組一共來存儲輸入及中間迭代的結(jié)果，同時將輸入X送入第一行作為輸入：

import numpy as np
import matplotlib.pyplot as plt


#頻域抽取的基2FFT
loop_num= int(np.log2(N))
data=np.zeros((loop_num+1,N),dtype=np.complex)
data[0]=x

隨后開始FFT的迭代，循環(huán)變量i一共來表征迭代的次數(shù)；循環(huán)變量p用來表征每次循環(huán)將將數(shù)據(jù)換分為幾塊；循環(huán)變量j用來進(jìn)行蝶形運(yùn)算。通過循環(huán)完成FFT的迭代及運(yùn)算，代碼如下：

for i in range(loop_num):
    k=i+1
    for p in range(2**i):
        for j in range(N//(2**k)):
            data[i+1][j           +p*(N//(2**i))] =  data[i,j+p*(N//(2**i))] + data[i,j+N//(2**k) +p*(N//(2**i))]
            data[i+1][j+N//(2**k) +p*(N//(2**i))] = (data[i,j+p*(N//(2**i))] - data[i,j+N//(2**k) +p*(N//(2**i))])*np.e**(-1j*2*j*np.pi*(2**i)/N)

最終將FFT蝶形運(yùn)算的結(jié)果進(jìn)行輸出倒序，定義rev2(k,N)遞歸函數(shù)達(dá)到按bit翻轉(zhuǎn)的目的，最終輸出FFT結(jié)果為fft_out：

def rev2(k,N):
    if (k==0):
        return (0)
    else:
        return(((rev2(k//2,N)//2)+(k%2)*(N//2)))


#輸出倒序
fft_out = np.ones_like(data[0,:])
for k  in range (N):
    fft_out[rev2(k,N)] = data[loop_num,k]

最后為了驗(yàn)證代碼正確性，直接調(diào)用python的FFT庫函數(shù)得到xf為庫函數(shù)的結(jié)果，與fft_out相減并畫圖，觀察誤差。

xf = np.fft.fft(x)
plt.plot(abs(xf))
plt.plot(abs(fft_out-xf))

輸入1024點(diǎn)的任意復(fù)數(shù)：

x = [int(np.round(np.sin(i)*1024))+int(np.round(np.cos(i)*1024))*1j for i in n]

波形如下：

運(yùn)行python算法得到結(jié)果如下，圖中藍(lán)線是FFT計算的結(jié)果，橙線是FFT庫函數(shù)計算結(jié)果與fft_out相減的差，差值為0，認(rèn)為我們的迭代算法正確。