基2FFT的算法推导及python仿真

FFT的算法推导主要用到旋转因子的周期性、对称性和可约性：

基2FFT的频域抽取法，将x(n)按照n的自然顺序划分为前后两个部分：

所以当K为偶数时，前后两部分相加。当k为奇数时相减。将频域X(K)划分为奇偶两个序列，N点DFT就被分解为两个N/2点的DFT：

可以得到蝶形图如下：

进而可以得到基2FFT频域抽取代码的实现方法：

随后是数据倒换，如下图：

可以看到基2FFT频域抽取后的输出位置排序就是自然数二进制码按位倒读的值。

根据推导结果我们编写python实现代码：

首先根据FFT的点数计算需要迭代的次数，根据迭代次数例化一个loop_num+1*N的数组一共来存储输入及中间迭代的结果，同时将输入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的迭代，循环变量i一共来表征迭代的次数；循环变量p用来表征每次循环将将数据换分为几块；循环变量j用来进行蝶形运算。通过循环完成FFT的迭代及运算，代码如下：

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蝶形运算的结果进行输出倒序，定义rev2(k,N)递归函数达到按bit翻转的目的，最终输出FFT结果为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]

最后为了验证代码正确性，直接调用python的FFT库函数得到xf为库函数的结果，与fft_out相减并画图，观察误差。

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

输入1024点的任意复数：

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

波形如下：

运行python算法得到结果如下，图中蓝线是FFT计算的结果，橙线是FFT库函数计算结果与fft_out相减的差，差值为0，认为我们的迭代算法正确。