使用 python 评估滤波器 matlab 函数

Question

根据这两个函数的描述，filter(b, a, x) matlab 函数和 lfilter(b, a, x, axis=-1, z_i=None) scipy.signal 应该给出相同的结果。

我举这个例子，结果完全不同：

import numpy as np
from scipy.signal import lfilter

bb = np.arange(0, 100, 1)

aa = 1

xx = np.tile(0.9, (100, 1))

yy = lfilter(b=bb, a=aa, x=xx)

print(yy[:10])

array([[0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.]])
print(yy.shape)
(100, 1)

# Same example on Matlab using filter
bb = 0:1:99
aa  =  1
xx =repmat(0.9, 100, 1)
dd = filter(bb, aa, xx)
dd(1:10)
ans =

         0
    0.9000
    2.7000
    5.4000
    9.0000
   13.5000
   18.9000
   25.2000
   32.4000
   40.5000
print(size(dd)) # (100 ,1)

Answer 1

MATLAB 和 NumPy 处理数组形状的方式不同。 NumPy 具有通用的 n 维数组。 scipy.signal.lfilter 接受一个 n 维数组，并沿输入数组的每个一维切片应用过滤器。哪个片是由axis参数决定的。默认情况下，lfilter 在 last 轴 (axis=-1) 上运行。您给 lfilter 一个形状为 (100, 1) 的数组。将 lfilter 和 axis=-1 应用于该输入应用过滤器 100 次，每行长度为 1 一次——当然不是您想要的！相反，您想沿 axis=0 应用过滤器（在这种二维情况下，意味着沿列应用 lfilter）。如果将 lfilter 的调用更改为

yy = lfilter(b=bb, a=aa, x=xx, axis=0)

返回的值将与 MATLAB 代码的值匹配。

从长远来看，我建议不要将自己局限于二维数组。在这种情况下，将 xx 创建为一维数组（例如 xx = np.full(100, fill_value=0.9)）更有意义。然后 lfilter 将按照您期望的方式对给定的一维数组进行操作，而无需指定 axis.

Answer 2

@warren 的回答很好。下面详细介绍函数的使用

filter(b, a, x, zi, dim) return zf 值，即使 zi 未作为输入给出。

在上面的例子中：

[y, zf] = filter(b=bb, a=aa, x=xx)

zf(1:10)

ans =

   1.0e+03 *

    4.4550
    4.4541
    4.4523
    4.4496
    4.4460
    4.4415
    4.4361
    4.4298
    4.4226
    4.4145

但是在 lfilter(a, b, x, zi=None, axis=-1) 中，要获得 zf，您必须提供 zi。为此，您可以使用 lfilter_zi(b, x) 为 lfilter 构造阶跃响应稳态的初始条件。

z_i = lfilter_zi(b=bb, a=aa)
y, zf = lfilter(b=bb, x=xx, a=aa, zi=z_i, axis=0)

zf[1:10]

[[4455. ]
 [4454.1]
 [4452.3]
 [4449.6]
 [4446. ]
 [4441.5]
 [4436.1]
 [4429.8]
 [4422.6]
 [4414.5]]

使用 python 评估滤波器 matlab 函数

Evaluate filter matlab function using python

python

matlab

scipy