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?

  • $\begingroup$ Can you accept one-directional IIR filtering or FIR filtering? $\endgroup$ – endolith Mar 21 at 17:37

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.

  • $\begingroup$ (I'm not sure about np.nditer. "The nditer will try to provide chunks that are as large as possible to the inner loop." But I don't know what "as large as possible" means.) $\endgroup$ – endolith Mar 21 at 18:06
  • 1
    $\begingroup$ I would figure out how to read the file in chunks, and write it out in chucks. You should only need to keep two input and two output chunks in memory at a time. $\endgroup$ – TimWescott Mar 21 at 19:01
  • $\begingroup$ @TimWescott Yeah I assume both in and out arrays would be mmap $\endgroup$ – endolith Mar 21 at 22:56
  • $\begingroup$ I would create my own FIR filter as a generator object in which you pass in a file object (which is an iterator) and write the output stream to a file. This approach is both simple once constructed and can stream indefinitely up to the limit of the file storage space available. $\endgroup$ – Dan Boschen Mar 21 at 23:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.