我构建了一个梳状倒滤波器,对既有信号滤波,理论上50Hz的整周期处,应该出现明显的FFT频谱幅度抑制,为什么实际滤波的效果这么差?理论上滤波器在中心点处有超过10倍的滤波效果。为什么没有表现出来?
'''
Demo comb filter effect.
but work terrible.
created by twicave @ Jun10,2023
V0.2.20230610
'''
#https://docs.scipy.org/doc/scipy/
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
# 定义信号,并加入白噪声。
fs = 1/(100e-6) # 采样频率
f_sigbase = 50;
f_sig = 50;
t = np.arange(0, 1/f_sigbase, 1/fs)
ppOfSig = 1000;
N = len(t);
print(N)
sig = np.random.normal(0, ppOfSig*0.50, len(t));
#追加一组与齿宽一致的方波信号:
f = 398 # 方波信号的频率
N = len(t) # 方波信号的采样点数
# 生成时间轴
Ts = 1 / (2 * f)
#t = np.linspace(0, (N-1)*Ts, N)
# 生成方波信号
duty = 0.5 # 方波信号的占空比,此处设置为50%
square_wave = ppOfSig * signal.square(2*np.pi*f*t, duty=duty)
sig = sig + square_wave;
sig_filtered = t;
# 设计并应用滤波器。
# 指定参数和采样频率
Q = 10 # 品质因数
# 计算滤波器的系数
numtaps = f_sigbase/fs; # 滤波器阶数
b, a = signal.iircomb(f_sigbase, Q, ftype='notch', fs=fs);
plt.scatter(b,a);
plt.show();
#滤波
sig_filtered = signal.lfilter(b,a, sig);
# 可视化信号、滤波器和滤波后的信号。
fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1, sharex=False, figsize=(6, 9))
ax1.scatter(t, sig, color='green', label="final signal");
ax1.plot(t, square_wave, color='red', label="square");
ax1.set_ylabel('Original Signal')
ax1.set_xlabel('Time (s)');
ax1.legend();
fft_sig = np.abs(np.fft.fft(sig)[:N//2])/N;
fft_sig[0] = 0;
fft_sig[1] = 0;
fft_sig[2] = 0;
ax2.scatter(np.arange(0, N//2)*fs/N, fft_sig);
ax2.set_ylabel('FFT of Original signal')
ax2.set_xlabel("Freq(Hz)");
ax3.plot(t, sig_filtered, color = 'green', label="filtered");
ax3.plot(t, sig, color='red', label="original signal");
ax3.legend();
ax3.set_ylabel('Filtered Signal')
ax3.set_xlabel('Time (s)')
fft_filted = np.abs(np.fft.fft(sig_filtered)[:N//2])/N;
fft_filted[0] = 0;
fft_filted[1] = 0;
fft_filted[2] = 0;
ax4.scatter(np.arange(0, N//2)*fs/N, fft_filted, color = 'green', label = "filtered")
ax4.scatter(np.arange(0, N//2)*fs/N, fft_sig, color = 'red', label = 'orignial signal')
ax4.legend();
ax4.set_ylabel('FFT of Filtered Signal')
ax4.set_xlabel("Freq(Hz)");
# 绘制频域响应
w, h = signal.freqz(b, a)
tomb_resp = 20*np.log10(abs(h))
tomb_freq = fs*w/(2*np.pi)
ax5.plot(tomb_freq[1:len(t)], tomb_resp[1:len(t)])
ax5.scatter(tomb_freq[1:len(t)], tomb_resp[1:len(t)])
ax5.set_title('Frequency response')
ax5.set_xlabel('Frequency (Hz)')
ax5.set_ylabel('Amplitude [dB]')
ax5.grid(True)
plt.ylim(-25,0);
plt.show()
发现自己接触不太熟悉的知识点时,脑子会糊涂。谱线就200点,滤波器不管用,我自己滤就完了。
把对应的周期倍数的FFT幅度全部干掉。然后转回时域。。。。
我可以人工进行逐点滤波。。。。
我居然想不到这个。
不过还是想知道为什么那个滤波器不行。
是不是它的频率响应曲线,画错了?它实际上没有那么陡峭的谷底?