[FIXED] vectorizing a "leaky integrator" in numpy

Issue

I need a leaky integrator -- an IIR filter -- that implements:

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

The following code works. However, my x vectors are long and this is in an inner loop. So my questions:

• For efficiency, is there a way to vectorize this in numpy?
• If not numpy, would it be advantageous to use one of the scipy.signal filter algorithms?

The iterative code follows. state is simply the value of the previous y[i-1] that gets carried forward over successive calls:

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

Solution

The current answer did not address this:

If not numpy, would it be advantageous to use one of the scipy.signal filter algorithms?

Yes, scipy.signal is well suited for IIR filtering through lfilter(), which is adequately vectorized. It takes input state and returns an output, output_state tuple as well, so for this particular case it would be:

from scipy.signal import lfilter

def leaky_integrator(x, state, leakiness):
    return lfilter([1], [1, -leakiness], x, zi=state)

Except that lfilter() expresses state (both when accepting and returning it) in a different way than your function. In this case, where your function would accept/return value, this function has [value*leakiness] (i.e. value is inside a 1-length array, and scaled by leakiness). See lfilter_zi() for more info.

Answered By - Alba Mendez