将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;                    %再存回endend
end
return_num=xn;

2. 部分原理

3. 代码移植

fft代码

def Myfft(xn,N):from math import log2from math import pifrom cmath import expimport numpy as npM=int(log2(N));          #计算级数length=len(xn);     #序列xn的长度if length<N:xn=xn+[0]*(N-length)   #如果xn序列长度不足N点FFT，则补零到长度为Nxn = 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 copydef 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 bfsizedef Dec2Bin_Inverse2Dec(n, m):'''Especially for fft.To find position.'''b = bin(n)[2:]if len(b) != m:b = "0"*(m-len(b)) + bb = 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=32else :T1=T[i]N=64xn=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()

将matlab的fft移植到python（蝶形运算）相关推荐

matlab中xtem,快速傅里叶变换_蝶形运算_按频率抽取基2-fft算法_MATLAB代码
function y=MyFFT_FB(x,n) %MYFFT_TB:My Fast Fourier Transform Frequency Based %按频率抽取基2-fft算法 %input: ...
Matlab 计算 FFT 的方法及幅值问题
欢迎转载,但请一定要给出原文链接,标注出处,支持原创! 谢谢~ https://blog.csdn.net/qq_29225913/article/details/105467006 目录 1.Mat ...
C语言实现FFT和IFFT，并与MATLAB编写显示的结果相对比，进行验证（蝶形运算）
本次实验中在Microsoft Visual Studio 2010环境下编写,实现FFT和IFFT,并用MATLAB编写显示的结果,两者相对比,进行验证. #include "stdafx ...
【Java】快速傅里叶变换FFT的程序实现（时间抽取的基-2FFT、倒位计算、蝶形运算）
Java--快速傅里叶变换(FFT)的程序实现好久没来更新了,阿汪大三了. 这学期阿汪要学习两门课<数字信号处理>和<Java程序设计>,刚好前几天老师告诉我们不久后会有个实 ...
matlab函数fftshift,matlab中fft算法_matlab中fftshift函数_matlab中fft函数的用法(2)
plot([0 : PointNum/2 - 1], x1(1:PointNum/2)); grid on subplot(3,1,2); % [REX IMX] am = sqrt(abs(REX. ...
Matlab中fft作频谱横纵坐标
关于这个问题,在很早之前就分享过,也通过了解实现了算法,当时看的明白,想的明白,突然要用的时候,又开始疑问,不免有些纠结,与其每次使用的时候都查,浪费时间,还不如,一次搞定. 真心没把哪门没学好的课程 ...
python变量存为matlab,详解如何在python中读写和存储matlab的数据文件(*.mat)
背景在做deeplearning过程中,使用caffe的框架,一般使用matlab来处理图片(matlab处理图片相对简单,高效),用python来生成需要的lmdb文件以及做test产生结果.所以 ...
快速傅立叶变换(FFT)算法(原来这就是蝶形变换)
快速傅立叶变换(FFT)算法(原来这就是蝶形变换) 为了实现FFT的海面模拟,不得不先撸个FFT算法实现. 离散傅立叶变换(DFT) 学习FFT之前,首先要先了解什么是DFT,我们都知道傅立叶变换是将 ...
示波器数据用matlab进行fft,示波器CSV波形数据导入Matlab进行FFT分析.doc
示波器CSV波形数据导入Matlab进行FFT分析 1,将CSV文件拖到workspace窗口,弹出的Import Wizard窗口中,点选"Next",新窗口中选第二项" ...

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

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

将matlab的fft移植到python（蝶形运算）相关推荐

最新文章

热门文章