SciPy lfilter 用于任意过滤器，初始条件沿 N 维数组的任何轴应用

Question

这个pythonnumpy/Scipy问题有一个答案在运算axis!=0或当初始条件向量的长度zi 保持为 1，即对于一阶滤波器。一旦高阶滤波器被定义为对分配给 signal 的第一个轴的输入 axis 进行操作，它就会失败，其中 signal 是 2D。

例如：

import numpy as np
from scipy.signal import (lfilter, lfilter_zi, butter)

def apply_filter(B, A, signal, axis=-1):
   # apply filter, setting proper initial state (doesn't assume rest)
   filtered, zf = lfilter(B, A, signal, 
             zi= lfilter_zi(B, A) * np.take(signal, [0], axis=axis), axis=axis)
   return filtered

B, A = butter(1, 0.5)
x = np.random.randn(12, 50)
apply_filter(B, A, x, axis=1)    # works without error
apply_filter(B, A, x, axis=0)    # works without error


B, A = butter(2, 0.5)
x = np.random.randn(12, 50)
apply_filter(B, A, x, axis=1)    # works without error
apply_filter(B, A, x, axis=0)    # raises error

加注

ValueError: operands could not be broadcast together with shapes (2,) (1,50)

如何使解决方案更通用以适用于 axis=0 的任意滤波器长度？

Answer 1

当*这两个数组时，您的错误来自数组的形状。简单示例：

>>> np.random.rand(2) * np.random.rand(1,10) 
----> 1 np.random.rand(2) * np.random.rand(1,10)

ValueError: operands could not be broadcast together with shapes (2,) (1,10) 


>>> np.random.rand(2)[:,None] * np.random.rand(1,10) 
array([[0.05905608, 0.30028617, 0.12495555, 0.28012041, 0.15031258,
        0.05166653, 0.2035891 , 0.01499304, 0.31749996, 0.3146938 ],
       [0.06860488, 0.34883958, 0.14515967, 0.3254132 , 0.17461669,
        0.06002052, 0.23650752, 0.01741727, 0.36883668, 0.36557679]])

为了在函数中处理不同的轴，您可以编写如下函数：

def apply_filter(B, A, signal, axis=-1):
    tmp = lfilter_zi(B, A) if axis==1 else lfilter_zi(B, A)[:,None]
    filtered, zf = lfilter(B, A, signal, zi = tmp * np.take(signal, [0], axis=axis), axis=axis)
    return filtered

SciPy lfilter 用于任意过滤器，初始条件沿 N 维数组的任何轴应用

SciPy lfilter for arbitrary filter, with initial conditions applied along any axis of N-D array

python

numpy

filter

scipy