Python - 怎么做？在频率（FFT）或时间（filtfilt）中过滤掉多个频段

Question

我有一个时间序列信号，我想过滤掉不是我的基频倍数的任何东西。

我认为数字梳状滤波器是一个很好的起点，但它并没有在我所有的基频上给我一个 1.0 的系数，我不确定为什么。

from scipy import signal
from matplotlib import pyplot as plt

MySignal = [some time series data....]
BaseFreq = 90.0
fs = 36000
F0 = BaseFreq
F1 = BaseFreq * 2
Q0 = 50.0

# Basic Notch Filter for filtering out even harmonics
b,a = signal.iircomb(F0*2, Q0, ftype='notch', fs=fs)
freq, h = signal.freqz(b,a, fs=fs)
response = abs(h)

# Basic Peak Filter for
d,c = signal.iircomb(F0, Q0, ftype='peak', fs=fs)
freq2 ,h2 = signal.freqz(d,c, fs=fs)
response2 = abs(h2)

# Plots
fig, ax = plt.subplots(2, 1, figsize=(12,6),sharex=True)
ax[0].plot(freq, response)
ax[1].plot(freq2, response2)
plt.show()

我尝试这样做的另一种方法是在频域中使用 FFT，高、低和带通都很容易做到，但我似乎找不到一个很好的梳子来正确地剔除非基频。

我的另一个想法是在我的基频周围应用一个 Kaiser 窗口，但我不确定我什至可以在频域中如何使用它。

感谢和想法，即使它只是一个链接去做一堆阅读！

Answer 1

你的滤波器确实有 1.0 的增益峰值，但绘图分辨率太低，看不到它。默认情况下，scipy.signal.freqz 在 512 点对频率响应进行采样。这通常是不够的，尤其是解决响应中的窄槽口或峰值时。尝试增加到 10K 点，例如

signal.freqz(b,a, worN=10000, fs=fs)

.

还有另一个问题：当第一个滤波器根据需要对偶次谐波进行陷波时，峰值滤波器偏离预期响应的一半倍数。解决方法是使用“

pass_zero = True

”选项。如文档中所述，

pass_zero : bool, optional
如果为 False（默认值），则滤波器的陷波（空点）以频率

[0, w0, 2*w0, ...]

为中心，峰值以中点

[w0/2, 3*w0/2, 5*w0/2, ...]

为中心。如果为真，则峰值以

[0, w0, 2*w0, ...]

（通过零频率）为中心，反之亦然。

（旁注：

pass_zero

选项相对较新，在 1.9.0 版本中添加。如果您的安装缺少它，您可以简单地翻转

b

和

a

的最后一个系数上的符号。）

这是我通过这些更改得到的结果。我放大了前几个倍数。

代码：

# Copyright 2023 Google LLC.
# SPDX-License-Identifier: Apache-2.0

BaseFreq = 90.0
fs = 36000
F0 = BaseFreq
Q0 = 50.0

# Notch filter for removing the even harmonics.
b,a = signal.iircomb(F0*2, Q0, ftype='notch', fs=fs)
freq, h = signal.freqz(b,a, worN=10000, fs=fs)
response = abs(h)

# Peak filter for passing only multiples of F0.
d,c = signal.iircomb(F0, Q0, ftype='peak', fs=fs)
d[-1] *= -1  # Implement `pass_zero=True` pre-1.9.0.
c[-1] *= -1
freq2, h2 = signal.freqz(d,c, worN=10000, fs=fs)
response2 = abs(h2)

fig, ax = plt.subplots(2, 1, figsize=(12,6),sharex=True)
ax[0].plot(freq, response)
ax[0].set_xlim(0, 280)
ax[0].set_ylabel('Gain')
ax[1].plot(freq2, response2)
ax[1].set_xlim(0, 280)
ax[1].set_xlabel('Frequency (Hz)')
ax[1].set_ylabel('Gain')
plt.show()