Memory efficient filtering with scipy.signal in Python

Question

I have several EEG channels that need filtering, and one channel may reach a recording length of 12 hours or potentially even more. This results in billions of data points stored which, when loaded, is represented as a NumPy array (by default as 16 bit long integer). In order not to load all data into memory at once, I make use of Python's mmap. This solution works until I call scipy.sosfiltfilt(). As I observed, this method creates copies of the entire signal which gets eventually filtered and is returned. The problem arises with really large arrays when making copies becomes simply unfeasible in terms of memory usage.

I have done some research, however, I have not come across a solution. Therefore, my question is whether there is any solution to filter such lengthy signals efficiently. My goal is to achieve this by for example, using an in place updating approach instead of creating copies without having to reimplement the method. Or does scipy.signal support mmap in any way? Or should I just create smaller chunks of the data? Or is there any other feasible way to do this?

Answer 1

The way to do this is to break the signal up into chunks and process each one at a time.

I asked in a comment if you could accept one-dimensional filtering, but I guess you can do bidirectional filtering, too; you just have to do it yourself.

For one-directional filtering:

• Break up the signal into chunks that can fit in memory (Is there a convenience function for this? Maybe np.nditer?)
• zi is initially all zeros
• For each chunk:
  • Run it through scipy.signal.sosfilt, with zi=zi
  • Set zi = zf (from the output of the filter)
  • Write the chunk of output filtered data to your new array

When you've assembled all the output chunks, it will be as if you had done one continuous sosfilt in one direction.

For bi-directional filtering:

• Do the whole thing again, but in the reverse direction.

For this to be exactly equal to sosfiltfilt, you probably need to keep some additional tail past the end of your signal, and then use that as your initial conditions when going backwards? Look at the code for sosfiltfilt or filtfilt and confirm how they work.

If you are ok with FIR filtering, then that's easier: You can make a zero-phase filter and just do that in one direction (using lfilter instead of sosfilt), with the same process as above, and then throw away the head of the output, so that it looks as if it were done with no offset.