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)

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

其实，挺为这个感到骄傲的。