(1)原始信号包含50Hz/100Hz两个频率成分,命名为sig。
(2)带噪信号是在sig的基础上增加一个服从正态分布的随机噪声,命名为sig_noise
(3)对带噪信号sig_noise做FFT变换(sig_noise_fft),先求取FFT变换后的abs值(sig_noise_fft_abs),并绘图。
(4)滤波:对sig_noise_fft_abs数组抽样,抽取abs<0.3(噪声)样本的下标索引,将抽取的数组索引对应的FFT值赋0,去噪并绘图
(5)对信号做IFFT变换并绘图与原始信号对比。
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['Microsoft YaHei']
sample_rate = 1024
t = np.linspace(0, 1, sample_rate)
fre1 = 50
fre2 = 100
sig = np.cos(2*np.pi*fre1*t) + np.cos(2*np.pi*fre2*t) # 原始信号
sig_noise = sig + np.random.normal(0, 0.6, len(sig)) # 带噪信号
# 绘制原始信号
plt.figure(figsize=(20, 10))
plt.subplot(2, 2, 1)
plt.plot(t, sig, color='blue')
plt.xlim(0, 0.1)
plt.xlabel("time(s)")
plt.title("原始信号=blue,去噪后IFFT信号=red")
# 绘制带噪信号
plt.subplot(2, 2, 2)
plt.plot(t, sig_noise, color='green')
plt.xlim(0, 0.1)
plt.xlabel("time(s)")
plt.title("带噪信号")
# 计算并绘制带噪信号的FFT频域图
n_fft = 1024
fre = np.linspace(0, sample_rate/2, int(n_fft/2)+1)
sig_noise_fft = np.fft.rfft(sig_noise, n_fft)
sig_noise_fft_abs = np.abs(sig_noise_fft)*2/n_fft
plt.subplot(2, 2, 3)
plt.plot(fre, sig_noise_fft_abs)
plt.xlim(0, 513)
plt.xlabel("Fre")
plt.title("带噪信号FFT")
# 对带噪信号滤波并绘制FFT频域图
found_fre = np.where(sig_noise_fft_abs < 0.3) # 滤波:判定数组内abs<0.3的值为噪声,并抽取其数组索引
filter_sig_noise_fft = sig_noise_fft.copy()
filter_sig_noise_fft[found_fre] = 0 # 滤波:将抽取的数组索引对应的FFT值赋0,去噪
filter_sig_noise_fft_abs = sig_noise_fft_abs.copy()
filter_sig_noise_fft_abs[found_fre] = 0
plt.subplot(2, 2, 4)
plt.plot(fre, filter_sig_noise_fft_abs)
plt.xlim(0, 513)
plt.xlabel("Fre")
plt.title("滤波后信号FFT")
# 对滤波后的带噪信号进行IFFT并绘制时域图与原信号对比
filter_sig_noise_ifft = np.fft.irfft(filter_sig_noise_fft, 1024)
plt.subplot(2, 2, 1)
plt.plot(t, filter_sig_noise_ifft, color='red')
plt.xlim(0, 0.1)
plt.show()
代码运行结果:
思考:带噪信号FFT去噪后再IFFT的波形似乎和原始波形会有微小的差异?因为本文的滤除随机噪声仅针对频域滤除,但是有部分噪声信号依然叠加在了50Hz和100Hz的频率成分上,造成了这个差异。