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,认为我们的迭代算法正确。
全部0条评论
快来发表一下你的评论吧 !