将matlab的fft移植到python（蝶形运算）

将matlab的fft移植到python（蝶形运算）

1. 缘起

matlab源代码由一名大佬编写而成，其中的原理是蝶形运算。（我不太懂只能移植过来用了，以此加深对fft的理解）。
该文章链接：Peter_831

%----------------函数调用方法：Myfft(xn,N)------------------
%----------------xn为要进行FFT计算的序列，N为FFT点数=2^n------------
function return_num=Myfft(xn,N)   

M=log2(N);          %计算级数
len=length(xn);     %序列xn的长度

if(len<N)
    xn=[xn,zeros(1,N-len)];  %如果xn序列长度不足N点FFT，则补零到长度为N
end

%--------------------倒序部分调用自带函数fliplr------------------------
binary_index=dec2bin((0:N-1),M);  %将输入的n=0,1,...N-1转换成二进制
reverse_binary_index=fliplr(binary_index);  %将二进制序列逆序
reverse_decimal_index=bin2dec(reverse_binary_index)+1; %逆序完将索引值转为十进制，记得加1，matlab数组索引从1开始到N结束
xn=xn(reverse_decimal_index);   %按照新的索引值重新排序，完成逆序输入

%--------------------分级进行蝶形运算---------------------------------
for L=1:M                               %DIT-FFT变换，M级蝶形变换
    B=2^(L-1);                          %两个输入数据相距B个点
    for J=0:B-1;                        %旋转因子处理
        P=2^(M-L)*J;                    %旋转因子上标
        W=exp(-1i*2*pi*P/N);            %对应旋转因子
        for k=(J+1):2^L:N;              %每级相同旋转因子间隔，数组索引从1开始到N结束，每次循环间隔2^L点          
            T=xn(k)+xn(k+B)*W;          %进行蝶形运算,先暂存，此时的xn(k)还要参与下一个xn(k+B)计算，不能更新存储值
            xn(k+B)=xn(k)-xn(k+B)*W;    %等相减部分计算完存储后
            xn(k)=T;                    %再存回
        end
    end
end
return_num=xn;

2. 部分原理

3. 代码移植

fft代码

def Myfft(xn,N):
    from math import log2
    from math import pi
    from cmath import exp
    import numpy as np
    M=int(log2(N));          #计算级数
    length=len(xn);     #序列xn的长度
    if length<N:
        xn=xn+[0]*(N-length)   #如果xn序列长度不足N点FFT，则补零到长度为N
    
    xn = np.array(xn)
    xn=xn.astype(np.complex)#特别注意astype改完数据类型要重新赋值
    
     #--------------------倒序操作---------------------------------
    tn=np.array([0]*N)
    tn=tn.astype(np.complex)
    for i in range(N):
        tn[i]=xn[i]
            
    for i in range(N):
         t = bin(i)           #%将输入的n=0,1,...N-1转换成二进制
         t=int(t[:2] + t[2:].zfill(M)[ : :-1], 2) #将二进制序列逆序
         xn[i]=tn[t]             #按照新的索引值重新排序，完成逆序输入              

         
 
    
    #--------------------分级进行蝶形运算---------------------------------

    for L in range(1,M+1):                               #DIT-FFT变换，M级蝶形变换
        B=2**(L-1)                         #两个输入数据相距B个点
        for J in range(B):                       #旋转因子处理
            P=2**(M-L)*J                    #旋转因子上标
            W=  np.array(exp(complex(0,(-2*pi*P/N))))            #对应旋转因子
            W=W.astype(np.complex)
            
            for k in range(J,N,2**L):             #每级相同旋转因子间隔，数组索引从1开始到N结束，每次循环间隔2^L点          
            #for k in range(J,2**(M-L)):   
                T=xn[k]+xn[k+B]*W         #进行蝶形运算,先暂存，此时的xn(k)还要参与下一个xn(k+B)计算，不能更新存储值
                xn[k+B]=xn[k]-xn[k+B]*W    #等相减部分计算完存储后
                xn[k]=T                    #再存回
                
    return xn

测试代码

from matplotlib.pyplot import stem
import numpy as np
from numpy.fft import *
from Myfft import *
import matplotlib.pyplot as plt
xn=[0,1,2,3,4,5,6,7]
y1=Myfft(xn,8)
y2=fft(xn,8)
plt.figure(1)
plt.subplot(211)
stem(y1)
plt.subplot(212)
stem(y2)

对比一下别人手打的代码。（我不配写出来这么优雅的代码）
数字图像处理的python实践(11)——傅里叶变换和快速傅里叶变换

import numpy as np
import math
import copy
 
def fft1(src, dst = None):
	'''
	src: list is better.One dimension.
	'''
	l = len(src)
	n = int(math.log(l,2))
 
	bfsize = np.zeros((l), dtype = "complex")
 
	for i in range(n + 1):
		if i == 0:
			for j in range(l):	
				bfsize[j] = src[Dec2Bin_Inverse2Dec(j, n)]
		else:
			tmp = copy.copy(bfsize)
			for j in range(l):
				pos = j%(pow(2,i))
				if pos < pow(2, i - 1):
					bfsize[j] = tmp[j] + tmp[j + pow(2, i - 1)] * np.exp(complex(0, -2*np.pi*pos/pow(2,i)))
					bfsize[j + pow(2, i - 1)] = tmp[j] - tmp[j + pow(2, i - 1)] * np.exp(complex(0, -2*np.pi*pos/(pow(2,i))))
	return bfsize
 
 
def Dec2Bin_Inverse2Dec(n, m):
	'''
	Especially for fft.To find position.
	'''
	b = bin(n)[2:]
	if len(b) != m:
		b = "0"*(m-len(b)) + b
	b = b[::-1]
	return int(b,2)
 
src = [1,2,3,4,5,6,7,8]
print(ifft1(fft1(src)))

4. 调试总结
1、Python 的数组下标起始值是从0 开始，而matlab是从1开始，Python 还有负号的下标。这也就意味着 Python是0到(n-1),而matlab是从1到n;在for循环的时候要特别注意。
2、特别注意^在Matlab 中是次方的意思，可是在Python 中的次方是** ，在Python中^的作用是异或。
3、Matlab 用’%’ 符号代表程序中的注释，Python用’#’ 注释。
4、在Matlab中，；号表示是否打印语句结果，可以用来调试。Python得用print()来显示结果。
5、Matlab的矩阵运算可以直接操作，而Python进行矩阵运算只能借助于特定的库Numpy。
6、Matlab对复数很友好(5+4i)，在Python中乖乖用complex吧。
7、Matlab对数组的索引是a(n) ,在Python数组索引是用中括号[… ] 表示，在元组是用小括号(… ) 表示。（这点很重要）
8、Numpy使用注意事项 ：特别注意astype改完数据类型要重新赋值才能确保数据类型更改完成。
9、Python的math库对复数运算兼容不好，使用cmath库进行替代。当然还可以使用np.sin()操作。（哈哈后来才知道）。

----------------------------------分割线----------------------------------------
补上另一份测试代码（有点小问题）（还是移植（抄））

from math import pi
from cmath import sin 

F=50;		#输入信号频率
N=32;		#N点FFT
T=np.array([0.000625,0.005,0.0046875,0.004,0.000625])		#抽样间隔
n=[i for i in range(64)]   #对连续时间信号x(t)采样64个点   #采样点数应足够大
n=np.array(n)
fig=plt.figure(1)
for i in range(5):

    if i<4:
        T1=T[i]
        N=32
    else :
        T1=T[i]
        N=64
    xn=np.sin(2*pi*F*n*T1)   #将x(t)抽样成离散x(nT)
   # print(xn)
    Xk=fft(xn,N)   #对抽样得到的xn进行FFT变换，此处采用自己编写的FFT函数，可换成系统自带函数Xk=fft(xn,N);
    #print(Xk)
    abs_Xk=abs(Xk)  #取模值
    xlabel_f=np.array([j for j in range(N)])/(N*T1)
    print(xlabel_f)
    plt.subplot(3,2,i+1)   
    plt.grid()  # 生成网格子
    stem(xlabel_f,abs_Xk[:N])  #由于实序列的dft共轭对称，即模值关于k=N/2偶对称，所以只需查看前N/2个点即可  
    plt.xlabel('频率/Hz')
    plt.ylabel('振幅')
    plt.title(['F=50'+'  '+'N='+str(N)+'  '+'T=',str(T1)])

fig.tight_layout()