1024基4FFT python实现

import os
from scipy.fftpack import fft,ifft
from math import log
from math import cos, sin,pi


def but_f(x1,x2,x3,x4):
    """
蝶形运算
    """
    x_1 = x1 + x2
    x_2 = x1 - x2 
    x_3 = x3 + x4
    x_4 = x3 - x4
    x_4 = x_4  * -1j
    y1 = x_1 + x_3
    y2 = x_2 + x_4
    y3 = x_1 - x_3
    y4 = x_2 - x_4 
    return y1,y2,y3,y4 

def w_nk(n):
    """
    定义选择因子
    """
    result =  cos((n/1024) * 2 * pi) - (sin((n/1024)* 2 * pi) * 1j)
    return result


def to_bin(value, num):#十进制数据，二进制位宽
	bin_chars = ""
	temp = value
	for i in range(num):
		bin_char = bin(temp % 2)[-1]
		temp = temp // 2
		bin_chars = bin_char + bin_chars
	return bin_chars.upper()

def r4_1024fft(n):
    """
    1024基4fft
    """
    stages = 5

    a = [0]*len(n)
    """
    逆序排列地址 0 - 1024  转换成 0 512 256 768 ...
    """
    for i in range(1024):
        res_i_1 = ""
        res_i = to_bin(i,10)[::-1]
        res_i = int(res_i,2)
        # res_i = list(res_i)
        #res_i[0:2],res_i[4:6] = res_i[4:6] ,res_i[0:2]
                   
        # for j in range(6):
        #     res_i_1 = res_i_1 + str(res_i[j])
        # res_i_1 = int(res_i_1,2)
        a[i] = n[res_i]        
        
    for stage in range(stages):#第一级无需成旋转因子， 第二级16个数据一组 第三级 64 个数据一组
        count_buf = 4 ** (stage) 
        
        groups = int(1024 / (4 ** (stage+1)))
        n =  4 ** stage
        m = int(log(4**(stage + 1),2))
        
        p = 0
        for i in range(groups):#256 ，64，16，4，1
            p = (4 ** (stage+1)) * i  

            for j in range(count_buf):#1 #4,16，64，256
                
                if(stage == 0):
                    wn1_bin =0
                    wn2_bin =0
                    wn3_bin = 0
                    n1 = 0
                    n2 = 0
                    n3 = 0
                    k1 = 0
                    k2 = 0
                    k3 = 0
                    wn1 = 0
                    wn2 = 0
                    wn3 = 0
                else:
                    wn1_bin = to_bin(p+n+j,10)[::-1]
                    wn2_bin = to_bin(p+n*2+j,10)[::-1]
                    wn3_bin = to_bin(p+n*3+j,10)[::-1]
                    n1 = wn1_bin[m-2:m]
                    n2 = wn2_bin[m-2:m]
                    n3 = wn3_bin[m-2:m]
                    k1 = wn1_bin[m-3::-1]
                    k2 = wn2_bin[m-3::-1]
                    k3 = wn3_bin[m-3::-1]
                    n1 = int(n1,2)
                    n2 = int(n2,2)
                    n3 = int(n3,2)
                    k1 = int(k1,2)
                    k2 = int(k2,2)
                    k3 = int(k3,2)

                    wn1 = n1*k1*groups
                    wn2 = n2*k2*groups
                    wn3 = n3*k3*groups




                a[p+n+j] = w_nk(wn1) * a[p+n+j]
                a[p+n*2+j] = w_nk(wn2)*a[p+n*2+j] 
                a[p+n*3+j] = w_nk(wn3)*a[p+n*3+j]            
                a[p+j],a[p+n+j],a[p+n*2+j],a[p+n*3+j] = but_f(a[p+j],a[p+n+j],a[p+n*2+j],a[p+n*3+j]) 
                
                        
    return a 
if __name__ == '__main__':
    a = [1]*1024
    for i in range(1024):
         a[i]=i
    c = r4_1024fft(a)
    print(c[:20])
    print(len(c))
    b = list(fft(a))
    print(b[:20])
   

直接运算与使用包的方法一致，每一个蝶形的初始地址产生由p+i+j

蝶形运算的旋转因子是根据地址得到的，具体方式需要转换成二进制数进行倒序处理具体做法是根据一篇论文：可在知网下载基于CORDIC算法的高性能FFT设计与实现--《南开大学》2009年硕士论文 (cnki.com.cn)

 论文中每一级的旋转因子的分母不同，(16,64,...) 但我的方法中选择因子的分母是固定1024，所以获取到选择因子地址后还需乘以一个倍数，即代码中的 groups

蝶形运算采用的方式是：

 现阶段是理论实现阶段，以此博客记录，接下来开始着手fpga实现，采用基4的方法可以加快速度（基2有10级，基4有5级）

结果：上半部分为自己方法结果，下半部分是fft包的结果