2

is there a way to have an alternative implementation of lfilter within scipy? I want to use the cusignal library and lfilter is not supported at the moment.

Here's my my current code that I want to speed up:

from scipy import signal
import numpy as np

data = np.random.rand(192,334)
a = [1,-1.086740193996892,0.649914553946275,-0.124948974636730]
b = [0.054778173164082,0.164334519492245,0.164334519492245,0.054778173164082]

x[range(0, len(x)),:] = signal.lfilter(b, a, x[range(0, len(x)),:])

Is there a way I can use numpy's convolve function or scipy's fftconvolve or firfilter to perform this operation? Ultimately, I want to perform the code snippet above faster than it's current version.

Any ideas or thoughts would be appreciated!

DumbCoder21
  • 113
  • 2
  • 7

1 Answers1

4

Even though your filter is in principle an infinite impulse response (IIR) filter, the impulse response for this particular filter decays very fast. You can compute the impulse response by running an impulse through it with lfilter like lfilter(b, a, [1] + [0]*99). Here is what I get:

impulse response plot1

As you can see, the taps are nearly zero above sample 20 or so. So you can take the first 20 samples of the impulse response to make an accurate truncated FIR approximation. From there, you can apply that FIR approximation with any FIR filtering function, like np.convolve, scipy.signal.convolve, or scipy.signal.fftconvolve.

Another thought: With any of these filtering functions, you could try casting all the args to np.float32. They might internally switch to a 32-bit float implementation that is faster than the 64-bit float implementation.

Pascal Getreuer
  • 2,906
  • 1
  • 5
  • 14
  • Thank you for your post! So taking what you said, how would I use scipy.signal.convolve or fftconvolve? Here's what I have so far: `impulse = signal.lfilter(b, a, [1] + [0]*99) impulse_20 = impulse[:20] conv = signal.convolve(data, impulse_20)` But I don't think that's entirely correct. – DumbCoder21 Oct 29 '20 at 15:09
  • [scipy.signal.convolve(x, y)](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve.html) returns (by default) the "full" convolution of length `len(x) + len(y) - 1`. To make it consistent with lfilter, trim off the excess elements from the right. Also, from your post it looks like `data` is a 2D array of shape (192, 334)? You'll need to either filter it row by row or work out how these functions operate on multidimensional inputs. Try `signal.convolve(impulse_20[np.newaxis, :], data)[:, :-19]`. – Pascal Getreuer Oct 29 '20 at 17:20
  • Yeah, the data is a shape of (192,334). So I've tried this but it's not matching the output from the original lfilter version: `conv = np.zeros((334,192)) conv[range(0, len(data)),:] = signal.convolve(data[range(0, len(data)),:], impulse_20[np.newaxis, :])[:, :-19]` – DumbCoder21 Oct 29 '20 at 17:53
  • I believe the problem there is your dimensions got swapped: the result of lfilter has the same shape as the signal, (192, 334), not (334, 192). Try this: `conv = signal.convolve(impulse_20[np.newaxis, :], data)[:, :-19]` and compare to `y = signal.lfilter(b, a, data)` – Pascal Getreuer Oct 30 '20 at 03:35