Python信号处理小试牛刀——快速傅里叶变换(FFT)

输入：仿真一个理想的多频信号y，频率为3Hz、10Hz，然后在这个理想信号上添加一个白噪声，得到一个带有白噪声的多频信号y_noise；

处理过程：分别对两个信号进行快速傅里叶变换得到对应的频谱图；

输出：显示添加白噪声前后的时域图、频域图，并分析白噪声对信号时域图和频域图的影响。

Python：

import numpy as np
import scipy       #数据处理与分析
import matplotlib.pyplot as plt #结果显示
#-----信号仿真-----#
dt=0.01 #采样间隔
t=np.arange(0,3,dt) #采样时间序列
N=len(t) #采样点数
fs=1/dt #采样频率
n=np.arange(0,N,1)
f=fs/N*n #频谱图横轴频率序列
e=np.random.rand(1,len(t)) #白噪声
y=0.5*np.cos(2*np.pi*3*t+np.pi/3)+np.cos(2*np.pi*10*t+np.pi/3) #理想信号，y=0.5cos(2π*3*t+π/3)+cos(2π*10*t+π/3)
y_noise=np.array(y+e).ravel() #添加白噪声
#-----时域图显示-----#
plt.figure(1)
plt.rcParams['font.sans-serif'] = ['Songti SC']  # 配置以后可以显示中文
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
fig,ax=plt.subplots(2,1) #设置子图的数量
ax[0].plot(t,y,marker='o',markersize=0,color="b", linewidth=1, linestyle="solid", label="")
ax[0].set_title('理想信号',fontsize=15)#设置图的标题
ax[0].set_xlabel('时间(s)',fontsize=10)  #设置横轴名称以及字体大小
ax[0].set_ylabel('幅值',fontsize=10)   #设置纵轴
ax[0].legend(loc="best")#图例放到图中的最佳位置

ax[1].plot(t,y_noise,marker='o',markersize=0,color="b", linewidth=1, linestyle="solid", label="")
ax[1].set_title('理想信号+白噪声',fontsize=15)#设置图的标题
ax[1].set_xlabel('时间(s)',fontsize=10)  #设置横轴名称以及字体大小
ax[1].set_ylabel('幅值',fontsize=10)   #设置纵轴
ax[1].legend(loc="best")#图例放到图中的最佳位置

plt.tight_layout()  # 当有多个子图时，可以使用该语句保证各子图标题不会重叠
plt.show()
#-----时频转换-----#
y_fft=scipy.fftpack.fft(y) #理想信号FFT，得到虚数序列
Ay=abs(y_fft) #取模
Ayy=Ay*2/N #得到信号在频域对应的幅值
y_noise_fft=scipy.fftpack.fft(y_noise) #添加噪声后FFT，得到虚数序列
Ay_noise=abs(y_noise_fft) #取模
Ayy_noise=Ay_noise*2/N #得到信号在频域对应幅值
#-----频域图显示-----#
plt.figure(2)
plt.rcParams['font.sans-serif'] = ['Songti SC']  # 配置以后可以显示中文
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
fig,ax=plt.subplots(2,1) #设置子图的数量
ax[0].plot(f[0:int(N/2)],Ayy[0:int(N/2)],marker='o',markersize=0,color="b", linewidth=1, linestyle="solid", label="")
ax[0].set_title('理想信号频谱图',fontsize=15)#设置图的标题
ax[0].set_xlabel('频率(Hz)',fontsize=10)  #设置横轴名称以及字体大小
ax[0].set_ylabel('幅值',fontsize=10)   #设置纵轴
ax[0].legend(loc="best")#图例放到图中的最佳位置

ax[1].plot(f[0:int(N/2)],Ayy_noise[0:int(N/2)],marker='o',markersize=0,color="b", linewidth=1, linestyle="solid", label="")
ax[1].set_title('理想信号+白噪声频谱图',fontsize=15)#设置图的标题
ax[1].set_xlabel('频率(Hz)',fontsize=10)  #设置横轴名称以及字体大小
ax[1].set_ylabel('幅值',fontsize=10)   #设置纵轴
ax[1].legend(loc="best")#图例放到图中的最佳位置

plt.tight_layout()  # 当有多个子图时，可以使用该语句保证各子图标题不会重叠
plt.show()

运行结果：

① 时域图

 ② 频域图

 结果分析如下：

• 白噪声会让平滑的时域信号变得不平滑，有许多毛刺；
• 白噪声会使信号在频域上有个很大的直流分量；
• 去除白噪声的方法：高通滤波，也就是去除这个直流分量。

快速傅里叶变换(FFT)的原理讲解见如下链接：

傅里叶变换_Dfreedom.的博客-CSDN博客_csdn傅立叶变换