如何在 scipy.signal 中制作低通滤波器？

Question

我有几个关于在 python/scipy.signal 中制作低通滤波器的问题。非常感谢任何答案。

我正在尝试了解 scipy.signal.lfilter 的工作原理。它是否采用滤波器系数并将它们乘以数据值，以便 data[500]，它会做

for b in range(0,len(coeff)):
    filtered = filtered + data[500-b]*coeff[b]

这样做和它正在做的有什么区别？

我也不明白点击次数对开始过滤有何影响。对于不同数量的点击，我看到它从数据的不同值开始。我假设它只有在具有必要的系数后才会启动。在我的示例中，我有来自 scipy.signal.firwin 的 683 个系数，但滤波器从 300-400 开始，正如您在图像 Firwin lowpass filter 中看到的（滤波器为蓝色；正弦波为红色；x 从 0 开始-1000)

from scipy import signal

a = signal.firwin(683, cutoff = 1/30, window = "hamming")
t = signal.lfilter(a,1,sig)

fs=1，cutoff = fs/30，我得到一个带有 firwin 的低通滤波器，如上图所示，延迟很多。我可以做些什么来改善延迟？
改变采样率对滤波器有何影响？
我在网上找到了两种估算点击次数的方法：

import math

(2/3) * math.log10(1/(10*ripple*attenuation)) * fs/transition_width
((-10*math.log10(ripple*attenuation) - 13)/(14.6)) * fs/transition_width

哪个是更好的近似值？

任何澄清将不胜感激。

Answer 1

I'm trying to understand how scipy.signal.lfilter works.

scipy.signal.lfilter(b, a, x)实现"infinite impulse response" (IIR), aka "recursive", filtering，其中b和a代表IIR滤波器，x是输入信号。

b arg 是 M+1 分子（前馈）滤波器系数的数组，a 是 N+1 分母（反馈）滤波器系数的数组。按照惯例，a[0] = 1（否则可以对过滤器进行归一化以使其如此），所以我假设 a[0] = 1。第 n 个输出样本 y[n] 计算为

  y[n] = b[0] * x[n] + b[1] * x[n-1] + ... + b[M] * x[n-M]
                     - a[1] * y[n-1] - ... - a[N] * y[n-N].

IIR 滤波的特殊之处在于 y[n] 的公式取决于先前的 N 输出值 y[n-1], ..., y[n-N]；该公式是递归的。因此，为了开始这个过程，通常通过假设 y[n] 和 x[n] 对于 n < 0 为零来对过滤器进行“零初始化”。这就是 scipy.signal.lfilter 默认情况下所做的。

您还可以使用 scipy.signal.lfilter 通过设置 a = [1] 来应用“有限脉冲响应”(FIR)，就像您在问题 2 中所做的那样。这样就没有递归反馈项了过滤公式，所以过滤就变成了 b 和 x.

的简单卷积

I also don't understand how the number of taps affects when it starts filtering. For different number of taps I see it starting at different values on the data. I assumed it'll start only after it has the necessary coefficients. In my example I have 683 coefficients from scipy.signal.firwin, but the filter starts at 300-400, as you can see in the image Firwin lowpass filter (filter is in blue; sinewave is in red; x goes from 0-1000)

scipy.signal.lfilter 立即开始过滤。正如我上面提到的，它（默认情况下）假设 n < 0 的 y[n] 和 x[n] 为零，这意味着它将第一个输出样本计算为 y[0] 计算为

  y[0] = b[0] * x[0].

但是，根据您的过滤器，b[0] 可能接近于零，这可以解释为什么一开始似乎什么也没有发生。

检查滤波器行为的一个好方法是计算其“脉冲响应”，即查看将单位脉冲 [1, 0, 0, 0, ...] 作为输入传递所产生的输出：

plot(scipy.signal.lfilter(b, a, [1] + [0] * 800))

这是我得到的 b = firwin(683, cutoff = 1/30, window = "hamming")、a = [1]:

从这个图中我们可以看出几件事：脉冲响应起初很小，然后上升并振荡，其峰值在样本索引 341 处，然后再次对称衰减到零。该滤波器的延迟为 341 = 683 // 2，即您在设计滤波器时指定给 firwin 的抽头数的一半。

What can I do to improve the delay?

尝试将点击次数 683 减少到更小的值。或者，如果您不需要过滤 causal, try scipy.ndimage.convolve1d，这会移动计算以使过滤器居中：

scipy.ndimage.convolve1d(sig, firwin_filter, mode='constant')

How would changing the sampling rate affect the filter?

对于大多数滤波器设计，只要截止频率小于采样率的 1/4，精确的采样率几乎没有影响。或者换句话说，除非截止频率适度接近奈奎斯特频率，否则这通常不是问题。

I've found two methods online to approximate the number of taps.

我不熟悉这些公式。请注意，实现目标特性所需的抽头数取决于特定的设计方法，因此请注意这些公式假设的上下文。

在 scipy.signal 内，我建议将 kaiserord 与 firwin 或 firwin2 一起用于目标纹波量和过渡宽度。 Here is their example，其中 65 是以 dB 为单位的阻带纹波，width 是以 Hz 为单位的过渡宽度：

Use kaiserord to determine the length of the filter and the parameter for the Kaiser window.
>>> numtaps, beta = kaiserord(65, width/(0.5*fs))
>>> numtaps
167
>>> beta
6.20426
Use firwin to create the FIR filter.
>>> taps = firwin(numtaps, cutoff, window=('kaiser', beta),
                  scale=False, nyq=0.5*fs)

编辑： 对于其他设计，kaiserord 可能会进入正确的范围，但不要依赖于它一定会达到目标数量的波纹或过渡宽度.所以一个可能的通用策略可能是这样的迭代过程：

使用kaiserord获得点击次数的初步估计。
设计具有那么多抽头的滤波器。
使用scipy.signal.freqz得到滤波器的频率响应。
评估响应幅度与所需响应的接近程度。即，在通带上，计算与 1 的最大绝对差，在阻带上，计算与 0 的最大绝对差。
如果响应不够接近，请增加点击次数并return到第 2 步。否则，看看您是否可以通过减少点击次数来解决问题。

如何在 scipy.signal 中制作低通滤波器？

How to make a low pass filter in scipy.signal?

python

signal-processing

scipy

lowpass-filter