为什么这个数字滤波器的效果如此差?

我构建了一个梳状倒滤波器,对既有信号滤波,理论上50Hz的整周期处,应该出现明显的FFT频谱幅度抑制,为什么实际滤波的效果这么差?理论上滤波器在中心点处有超过10倍的滤波效果。为什么没有表现出来?

img

'''
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幅度全部干掉。然后转回时域。。。。
我可以人工进行逐点滤波。。。。
我居然想不到这个。

不过还是想知道为什么那个滤波器不行。
是不是它的频率响应曲线,画错了?它实际上没有那么陡峭的谷底?

不知道你这个问题是否已经解决, 如果还没有解决的话:

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^