python - 在 Python/Numpy 中实现具有不同系数的简单滤波器的有效方法

Question

我正在寻找一种有效的方法来实现一个简单的滤波器，该滤波器具有一个随时间变化的系数，并由与输入信号长度相同的向量指定。

以下是所需行为的简单实现:

def myfilter(signal, weights):
    output = np.empty_like(weights)
    val = signal[0]
    for i in range(len(signal)):
        val += weights[i]*(signal[i] - val)
        output[i] = val
    return output

weights = np.random.uniform(0, 0.1, (100,))
signal = np.linspace(1, 3, 100)
output = myfilter(signal, weights)

有没有办法使用 numpy 或 scipy 更有效地完成此操作？

Answer 1

您可以用循环的开销换取一些额外的操作:

import numpy as np

def myfilter(signal, weights):
    output = np.empty_like(weights)
    val = signal[0]
    for i in range(len(signal)):
        val += weights[i]*(signal[i] - val)
        output[i] = val
    return output

def vectorised(signal, weights):
    wp = np.r_[1, np.multiply.accumulate(1 - weights[1:])]
    sw = weights * signal
    sw[0] = signal[0]
    sws = np.add.accumulate(sw / wp)
    return wp * sws

weights = np.random.uniform(0, 0.1, (100,))
signal = np.linspace(1, 3, 100)
print(np.allclose(myfilter(signal, weights), vectorised(signal, weights)))

在我的机器上，矢量化版本快了好几倍。它使用递归方程的“封闭形式”解。

编辑:对于很长的信号/权重(例如 100,000 个样本)，此方法由于溢出而不起作用。在这种情况下，您仍然可以使用以下技巧节省一点(在我的机器上节省 50% 以上)，它的额外好处是您不需要求解递归公式，只需反转它。

from scipy import linalg

def solver(signal, weights):
    rw = 1 / weights[1:]
    v = np.r_[1, rw, 1-rw, 0]
    v.shape = 2, -1
    return linalg.solve_banded((1, 0), v, signal)

这个技巧利用了这样一个事实，即递归在形式上类似于只有一个非零下对角线的矩阵上的高斯消除。它依赖于专门从事这一工作的库函数。

事实上，我对此感到非常自豪。