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

yvt65v4c 于 2023-08-05 发布在其他

关注(0)|答案(2)|浏览(160)

免责声明 *：我可能不像我应该的那样擅长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文档中的示例代码中得到了它，但我不确定这是否正确。
然后就到了我不确定如何处理最初的几个样本的地步。系数a和b在这里是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()

型

scipy

来源：https://stackoverflow.com/questions/40483518/how-to-real-time-filter-with-scipy-and-lfilter

2条答案

按热度按时间

v64noz0r1#

问题不在于如何缓冲输入。问题在于，在“离线”版本中，滤波器的状态是使用lfilter_zi初始化的，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]

型

赞(0）回复(0）举报 2023-08-05

piok6c0g2#

我想我也遇到了同样的问题，并在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()

字符串
其思想是在每次后续调用中将过滤器状态z传递给lfilter。对于最初的几个样本，过滤器可能会给予奇怪的结果，但稍后（取决于过滤器长度）它开始正常运行。

赞(0）回复(0）举报 2023-08-05

我来回答

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

2条答案

相关问题

热门标签

最新问答