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

Question

关于在python/scipy.signal中制作低通滤波器，我有几个问题。如有任何答案，我将不胜感激

1. 我想了解^{}是如何工作的。它是否将滤波器系数与数据值相乘，这样对于data[500]，就可以了

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

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

1. 我也不明白当它开始过滤时，抽头的数量是如何影响的。对于不同的抽头数，我看到它从数据上的不同值开始。我假设它只有在有了必要的系数后才会开始。在我的示例中，我有683个来自scipy.signal.firwin的系数，但是过滤器从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)

1. 当fs=1，cutoff=fs/30时，我得到了一个带有firwin的低通滤波器，它延迟了很多，如上图所示。我能做些什么来改善延迟

2. 更改采样率将如何影响过滤器

3. 我在网上找到了两种近似抽头数的方法：

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

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是输入信号

barg是M+1分子（前馈）滤波器系数的数组，a是N+1分母（反馈）滤波器系数的数组。按照惯例，a[0]=1（否则可以对过滤器进行规范化，使其成为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的卷积

2. 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的y[n]和x[n]为零<；0，这意味着它计算第一个输出样本为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的抽头数的一半

3. What can I do to improve the delay?

尝试将轻敲683的数量减少到更小的数量。或者，如果不要求筛选为causal，请尝试scipy.ndimage.convolve1d，这会移动计算以使筛选器居中：

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

4. How would changing the sampling rate affect the filter?

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

5. 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可能会达到正确的标准，但不必依赖于这一标准来实现目标数量的波纹或过渡宽度。因此，一个可能的总体策略可能是这样一个迭代过程：

1. 使用kaiserord获得抽头数的初始估计值
2. 设计具有如此多抽头的滤波器
3. 使用scipy.signal.freqz获得滤波器的频率响应
4. 评估响应幅度与期望响应的接近程度。即，在通带上，从1计算最大绝对差，在阻带上，从0计算最大绝对差
5. 如果响应不够接近，增加抽头数并返回步骤2。否则，看看你是否可以通过减少轻敲次数来逃脱惩罚