numpy中“漏积分器”的矢量化

jjhzyzn0 于 2022-12-23 发布在其他

关注(0)|答案(1)|浏览(237)

我需要一个泄漏积分器--IIR滤波器--来实现：

y[i] = x[i] + y[i-1] * leakiness

下面的代码是有效的，但是，我的x向量很长，而且这是在一个内部循环中，所以我的问题是：

为了提高效率，有没有办法在numpy中将其矢量化？
如果不是numpy，使用scipy.信号滤波器算法之一是否有利？

迭代代码如下所示。state只是前一个y[i-1]的值，它在连续调用中被结转：

import numpy as np

def leaky_integrator(x, state, leakiness):
    y = np.zeros(len(x), dtype=np.float32)
    for i in range(len(x)):
        if i == 0:
            y[i] = x[i] + state * leakiness
        else:
            y[i] = x[i] + y[i-1] * leakiness
    return y, y[-1]

>>> leakiness = 0.5
>>> a1 = [1, 0, 0, 0]
>>> state = 0
>>> print("a1=", a1, "state=", state)
a1= [1, 0, 0, 0] state= 0
>>> a2, state = leaky_integrator(a1, state, leakiness)
>>> print("a2=", a2, "state=", state)
a2= [1.    0.5   0.25  0.125] state= 0.125
>>> a3, state = leaky_integrator(a2, state, leakiness)
>>> print("a3=", a3, "state=", state)
a3= [1.0625    1.03125   0.765625  0.5078125] state= 0.5078125

numpy

来源：https://stackoverflow.com/questions/74876955/vectorizing-a-leaky-integrator-in-numpy

1条答案

按热度按时间

aiazj4mn1#

我可以看到两个选择：

最简单的（也是建议的解决方案）是扩展依赖列表并使用numba。
使用矩阵乘法，从矩阵运算的Angular 重新思考问题。

事实上如果x=[a, b, c]，s=state，l=leakiness那么

y = [a + s*l, b + (a + s*l)*l, c + (b + (a + s*l)*l)*l]
  = [a + s*l, b + a*l + s*l**2, c + b*l + a*l**2 + s*l**3]
  = [[1, 0, 0], [l, 1, 0], [l**2, l, 1]] @ x + s * [l, l**2, l**3]

然而，您可能需要生成一个大小为x.size**2的矩阵，并且即使对于较小的大小也可能会出现内存不足错误（例如，对于1M大小的数组，它会导致大约7 TiB，我认为这是不可行的）。
回到numba实现，向已经实现的函数添加一个@jit(nopython=True)装饰器就足够了。
在我的机器上用一个大小为1M的随机数组执行此操作：

%timeit leaky_integrator(a1, s, l)
2.07 s ± 99.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit leaky_integrator_jitted(a1, s, l)
7.66 ms ± 22.2 µs per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

（此处为numba性能备注）

from numba import jit

@jit(nopython=True)
def leaky_integrator_jitted(x, state, leakiness):
    y = np.zeros(len(x), dtype=np.float32)
    for i in range(len(x)):
        if i == 0:
            y[i] = x[i] + state * leakiness
        else:
            y[i] = x[i] + y[i-1] * leakiness
    return y, y[-1]

赞(0）回复(0）举报 2022-12-23

我来回答

numpy中“漏积分器”的矢量化

1条答案

相关问题

热门标签

最新问答