用Python实现DFT并绘制功率谱

Zarah ·
更新时间:2024-09-21
· 509 次阅读

知识背景:傅里叶变换可以分为连续傅里叶变化和离散傅里叶变换,分别是FT,FS,DTFT,DTFS。其中DTFT是我们常说的离散时间傅里叶变换,但这种变换并不一定能够由计算机进行处理,因为对于非周期信号来说其谱一般是连续谱,这样就无法由计算机完成了。所以DFT就出现了,我们知道DTFT是以2π为周期的,我们一般只需要取其主值(可以看作取一个完整的周期)进行分析,如果对DTFT在0到2π内均匀尽行采样的话得到的结果就一定是离散的,如果我们的采样是遵循一定规则的那就可以用采样后的谱完整的恢复原信号。

完整的过程是:原信号遵循采样定理采样得到离散时间信号序列,假设其序列长度为n,如果我们的DFT的点数为K,只要K>=n,那么我们就可以完全重构其频谱(也就是要恢复其DTFT的谱),再由恢复的谱就可以恢复出原信号。再到后来对DFT的计算方式进行改进就有了FFT,是按照DFT具有周期性对称性把多点的DFT分为点数最少的DFT来进行计算。

FFT的算法实现比较复杂,我们一般都是调用第三方库中已经实现的通用的FFT,但是DFT这种不关心运算效率的算法实现起来就很简单,无非就是两层循环。
在这里插入图片描述
在这里插入图片描述
外层循环就是对K的遍历,内层循环就是对n的遍历,相乘,求和。
比如我们要对在这里插入图片描述
这样一个信号做DFT,f1=0.1 f2=0.3,实现过程如下:

from math import * import numpy as np import matplotlib.pyplot as plt #懒得导包了,直接全抓进来 def signal(n): return (sin(0.2*pi*n+pi/3)+10*sin(0.6*pi*n+pi/4)) #生成WN项 def wn_k(k,n,N): return complex(cos(2*pi*n*k/N),sin(-2*pi*n*k/N)) amplitude=[]#准备一个空列表 power_spectrum=[] sums=0 #64点DFT,X(0)到X(63) for k in range(0,64): for n in range(1,257): #n的取值为从1到256 sums=sums+signal(n)*wn_k(k,n,64) amplitude.append(sums) sums=0 print(amplitude,len(amplitude)) for i in range(0,64): power_spectrum.append(amplitude[i]**2) plt.suptitle("xn=sin(2*pi*f1*n+pi/3)+10*sin(2*pi*f2*n+pi/4)") plt.subplot(2,1,1) plt.plot(np.abs(amplitude)) plt.title("amplitude_spectrum") plt.subplot(2,1,2) plt.plot(np.abs(power_spectrum)) plt.title("power_spectrum") plt.show()

其中那一段for循环的结果就是把变换结果全部放入了一个列表,列表里存储的是复数。之所以要把结果放入列表是因为matplotlib的绘图函数的入口参数是列表类型。我们要绘制幅度谱的话就取列表里每个复数值的模值,绘制功率谱的话就把每个值平方一下。为了验证我们所求结果是否正确,我们可以调用numpy库的FFT来对照一下
库函数版本:

import numpy as np import matplotlib.pyplot as plt from math import * power_spectrum=[] x=np.linspace(1,256,256) y=np.sin(0.2*np.pi*x+np.pi/3)+10*np.sin(0.6*np.pi*x+np.pi/4) #plt.subplot(3,1,1) #plt.plot(y) #plt.title("input xn") trans=np.fft.fft(y,n=64) plt.suptitle("xn=sin(2*pi*f1*n+pi/3)+10*sin(2*pi*f2*n+pi/4)") plt.subplot(2,1,1) plt.plot(np.abs(trans)) plt.title("amplitude by numpy") for x in range(0,64): power_spectrum.append(trans[x]**2) plt.subplot(2,1,2) plt.plot(np.abs(power_spectrum)) plt.title("power_spectrum by numpy") plt.show()

对比结果如图:结果一致
在这里插入图片描述
请注意:实际上是离散谱,但绘图时自动把点与点之间用线段给连接起来了


作者:naruhina



用python 功率谱 功率 dft Python

需要 登录 后方可回复, 如果你还没有账号请 注册新账号