基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,认为我们的迭代算法正确。

仿真

打开APP阅读更多精彩内容
声明:本文内容及配图由入驻作者撰写或者入驻合作网站授权转载。文章观点仅代表作者本人,不代表电子发烧友网立场。文章及其配图仅供工程师学习之用,如有内容侵权或者其他违规问题,请联系本站处理。 举报投诉

全部0条评论

快来发表一下你的评论吧 !

×
20
完善资料,
赚取积分