CODEWITHSUNDEEP
pythonnumpyscipy
dani_l
dani_l

Reputation: 159

Efficient way to implement simple filter with varying coeffients in Python/Numpy

I am looking for an efficient way to implement a simple filter with one coefficient that is time-varying and specified by a vector with the same length as the input signal.

The following is a simple implementation of the desired behavior:

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)

Is there a way to do this more efficiently with numpy or scipy?

Upvotes: 2

Views: 200

Answers (1)

Paul Panzer
Paul Panzer

Reputation: 53099

You can trade in the overhead of the loop for a couple of additional ops:

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)))

On my machine the vectorised version is several times faster. It uses a "closed form" solution of your recurrence equation.

Edit: For very long signal / weight (100,000 samples, say) this method doesn't work because of overflow. In that regime you can still save a bit (more than 50% on my machine) using the following trick, which has the added bonus that you needn't solve the recurrence formula, only invert it.

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)

This trick uses the fact that your recurrence is formally similar to a Gauss elimination on a matrix with only one nonvanishing subdiagonal. It piggybacks on a library function that specialises in doing precisely that.

Actually, quite proud of this one.

Upvotes: 4

Related Questions