Python中低通黄油滤波器的问题

Question

我想过滤（低通）我拥有的信号，当它不起作用时，我开始调查它为什么不起作用。 我进行了一些测试，我对巴特沃斯滤波器的行为感到有些惊讶。 我在this post中定义了它

def apply_filter(data, cutoff, fs, order=6, filter_type="low", analog=False):
    nyq = 0.5 * fs
    normalized_cutoff = cutoff / nyq
    b,a = butter(order, normalized_cutoff, btype=filter_type, analog=analog, output="ba")
    they = lfilter(b, a, data)
    return(they)

如果我采取 1000 个元素长的样本，像这样

x = np.linspace(0, 2*np.pi, 1000)
y = np.sin(x) + 0.3* np.sin(10*x)
sampling_frequency = 1/ (x[-1] * 1e-3)
sampling_frequency
>> 159.15494309189532
# because i have 1000 thousand points for a "time" going up to 2 pi

plt.plot(x, y, x, apply_filter(y, cutoff=1, fs= sampling_frequency)

我得到的

另一方面，如果我做完全相同的事情，但点数不同，比如 10000，我得到一个错误的结果，我不太明白为什么：

x = np.linspace(0, 2*np.pi, 10000)
y = np.sin(x) + 0.3* np.sin(10*x)
sampling_frequency = 1/ (x[-1] * 1e-4)
sampling_frequency
>> 1591.5494309189535
# because i have 10000 thousand points for a "time" going up to 2 pi

plt.plot(x, y, x, apply_filter(y, cutoff=1, fs= sampling_frequency)

这一次，我得到了

这显然是错误的。 有人可以解释为什么会这样吗？在 1000 点或更少的情况下，一切似乎都正常......

编辑：

我已经绘制了滤波器的频率响应，问题出现在这些图表上，虽然我也不知道为什么会这样。

sampling rate 
>> 159.1549430918953
b, a = butter(6, 1/(sampling_rate/2))
w, h = freqz(b, a, 8000)
plt.subplot(2,1,1)
plt.xlim(0, 15)
plt.plot(0.5*sampling_rate*w/np.pi, np.abs(h))

我得到

如果我这样做了

sampling_frequency *= 10
sampling_frequency
>> 1591.5494309189535
b, a = butter(6, 1/(sampling_rate/2))
w, h = freqz(b, a, 8000)
plt.subplot(2,1,1)
plt.xlim(0, 15)
plt.plot(0.5*sampling_rate*w/np.pi, np.abs(h))

然后我得到

我觉得函数 butterworth 可能由于某种原因在高点数上遇到了一些问题？

感谢您的帮助！

Answer 1

对于任何可能感兴趣的人，这实际上是与“ba”输出一起使用的黄油过滤器的“已知问题”。您应该改用“zpk”输出，请参阅this link。

您可以以一种相当简单的方式使用“zpk”输出，这与使用“ba”输出的操作非常相似。这似乎适用于多达 100 万个点，没有理由不进一步工作。

这是一个基本的例子：

point_number=1000000
# our "data"
x = np.linspace(0, 2*np.pi, point_number)
y = sin(x) + 0.3* sin(10*x)

# sampling frequency would be 1/ sampling_time
sampling_frequency = point_number/(2*np.pi)
# hence the nyquist frequency
nyq = sampling_frequency/2
# desired cutoff frequency, in Hertz
cutoff = 1
# normalized for the function butter
normalized_cutoff = cutoff/nyq

z,p,k = butter(6, normalized_cutoff, output="zpk")
lesos = zpk2sos(z, p, k)
# filtered data
y_filtered_sos = sosfilt(lesos, y)

# plot part
plt.plot(x, y_filtered_sos, label="filtered data")
plt.title("filtered data, with zpk")
plt.plot(x,y, label="data")
plt.legend()
plt.title("filtered data, with zpk")

给了