如何使用 scipy 和 lfilter 进行实时过滤?

Posted

技术标签:

【中文标题】如何使用 scipy 和 lfilter 进行实时过滤?【英文标题】:How to real-time filter with scipy and lfilter? 【发布时间】:2017-03-21 21:05:50 【问题描述】:

免责声明:我在 DSP 方面可能没有我应该做的那么好,因此我遇到的问题比我应该让这段代码工作更多。

我需要在传入信号发生时对其进行过滤。我试图让这段代码工作,但到目前为止我还没有做到。 引用scipy.signal.lfilter doc

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
from lib import fnlib

samples = 100
x = np.linspace(0, 7, samples)
y = [] # Unfiltered output
y_filt1 = [] # Real-time filtered

nyq = 0.5 * samples
f1_norm = 0.1 / nyq
f2_norm = 2 / nyq
b, a = scipy.signal.butter(2, [f1_norm, f2_norm], 'band', analog=False)
zi = scipy.signal.lfilter_zi(b,a)
zi = zi*(np.sin(0) + 0.1*np.sin(15*0))

这会将zi 初始设置为zi*y[0 ],在本例中为0。我从lfilter 文档中的示例代码中得到它,但我不确定这是否正确。

然后到了我不确定如何处理少数初始样本的地步。 系数ab 在这里是len(a) = 5。 当lfilter 将输入值从现在到 n-4 时,我是用零填充它,还是需要等到 5 个样本过去并将它们作为一个块,然后在相同的下一步中连续采样怎么办?

for i in range(0, len(a)-1): # Append 0 as initial values, wrong?
    y.append(0)

step = 0
for i in xrange(0, samples): #x:
    tmp = np.sin(x[i]) + 0.1*np.sin(15*x[i])
    y.append(tmp)

    # What to do with the inital filterings until len(y) ==  len(a) ?

    if (step> len(a)):
        y_filt, zi = scipy.signal.lfilter(b, a, y[-len(a):], axis=-1, zi=zi)
        y_filt1.append(y_filt[4])

print(len(y))
y = y[4:]
print(len(y))
y_filt2 = scipy.signal.lfilter(b, a, y) # Offline filtered

plt.plot(x, y, x, y_filt1, x, y_filt2)
plt.show()

【问题讨论】:

你找到解决办法了吗? 是的,我在一年前的某个时候解决了这个问题。我可以看看我是否有时间重新挖掘它,因为我忘记了很多细节。 输出的差异与块大小无关。通过循环在每次迭代中提供多少样本并不重要。在实时版本(第二个代码块)中,过滤器的状态未初始化。在离线版本(第一个代码块)中,它使用lfilter_zi 进行初始化。 【参考方案1】:

我想我有同样的问题,并在https://github.com/scipy/scipy/issues/5116找到了解决方案:

from scipy import zeros, signal, random

def filter_sbs():
    data = random.random(2000)
    b = signal.firwin(150, 0.004)
    z = signal.lfilter_zi(b, 1) * data[0]
    result = zeros(data.size)
    for i, x in enumerate(data):
        result[i], z = signal.lfilter(b, 1, [x], zi=z)
    return result

if __name__ == '__main__':
    result = filter_sbs()

这个想法是在对lfilter 的每个后续调用中传递过滤器状态z。对于前几个样本,过滤器可能会给出奇怪的结果,但后来(取决于过滤器长度)它开始表现正常。

【讨论】:

我认为这是不对的。输出应该与一次性过滤整个数据集完全相同:result = signal.lfilter(b,1,data)。前几个值不应该有“奇怪的”。这是问题的核心。 @dmedine :感谢您的评论!答案中的代码给出的结果与signal.lfilter(b, 1, data, zi=z) 完全相同。 zi 是一个选择问题,但它应该确保results[0] == data[0](请参阅lfilter_zi)。没有,或者有z=zeros(b.size-1)results[0] 将接近于 0。而且我仍然认为我们总是会得到“奇怪”的第一个值,至于 result 的第一个值,过滤器只能查看第一个值data,因此没有足够的样本来正确过滤。 说'过滤器......没有足够的样本来正确过滤'是绝对不正确的。根据定义,DFII 过滤器(lfilter 就是这样)需要精确的n>0 输入样本来计算正确的输出。如果其状态未定义,它将失败,但这是一个单独的问题。我不明白你说的“奇怪”是什么意思。 LTI 系统根据输入和内部状态计算输出。矩阵的线性组合从来没有任何“奇怪”的地方。这只是数学,没有人可以反驳(不确定性原理无法承受)。 另外,zi 并不是真正的“选择”问题。必须有理由为zi 选择值。如果你的系统是因果的,你必须用 0 初始化过滤器,因为你看不到未来。在你的循环中,你用它之前的状态初始化你的过滤器,你知道它是因为它已经发生了。我不明白为什么你会在提供输入之前选择初始状态为半余弦。 哦,我现在明白了。 lfilter_zi 根据滤波器系数计算初始稳态响应。是的,这是初始化过滤器的一个很好的选择,但是 OP 仅对离线版本执行此操作,因此输出存在差异。【参考方案2】:

问题不在于您如何缓冲输入。问题在于,在“离线”版本中,过滤器的状态是使用lfilter_zi 初始化的,它计算 LTI 的内部状态,以便在新样本到达输入时输出已经处于稳定状态。在“实时”版本中,您可以跳过此操作,以便过滤器的初始状态为 0。您可以将两个版本都初始化为使用 lfilter_zi,或者将两者都初始化为 0。然后,无论您有多少样本一次过滤。

注意,如果您初始化为 0,过滤器将在达到稳定状态之前“响”一段时间。在 FIR 滤波器的情况下,有一个用于确定该时间的解析解决方案。对于许多 IIR 滤波器,没有。

以下是正确的。为简单起见,我初始化为 0 并一次输入样本。但是,任何非零块大小都会产生等效的输出。

from scipy import signal, random
from numpy import zeros

def filter_sbs(data, b):
    z = zeros(b.size-1)
    result = zeros(data.size)
    for i, x in enumerate(data):
        result[i], z = signal.lfilter(b, 1, [x], zi=z)
    return result
    
def filter(data, b):
    result = signal.lfilter(b,1,data)
    return result

if __name__ == '__main__':
    data = random.random(20000)
    b = signal.firwin(150, 0.004)
    result1 = filter_sbs(data, b)
    result2 = filter(data, b)
    print(result1 - result2)

输出:

[ 0.00000000e+00  0.00000000e+00  0.00000000e+00 ... -5.55111512e-17
  0.00000000e+00  1.66533454e-16]

【讨论】:

见上面的评论。

以上是关于如何使用 scipy 和 lfilter 进行实时过滤?的主要内容,如果未能解决你的问题,请参考以下文章

如何在python中对信号应用过滤器

使用 JAX 和 SciPy 对不正确积分进行微分

如何使用 Scipy 模拟饱和度和阈值?

使用空间和时间变量在 python(scipy) 中进行聚类

使用 Scipy 进行线性回归曲线拟合 - 不知道出了啥问题

python中如何安装SciPY模块