I have seen code examples of how to use the scipy Butterworth filter to filter data stored in a numpy array.
I need to do something slightly different. I want to process the data in blocks.
As an example, if the total signal is 4096 samples long, I would like to process the first 1024 samples and output the result. Then the next 1024 samples, and so on.
But when I concatenate the output arrays I would like the result to be exactly as if I had processed all 4096 samples in one operation.
My guess is that I will need to carry the state of the filter from one block to the next, but any pointers on how to do this would be very helpful.
I have used 1024 as an example. I could make the buffers larger or smaller if that helps. But I can’t put all the data in a single buffer because that would limit the total data I could process to the available memory and the data might be bigger than that in some cases.
My guess is that I will need to carry the state of the filter from one block to the next, but any pointers on how to do this would be very helpful.
Yes, that’s what you need to do.
Here’s an example of how to do this with scipy.signal.sosfilt(). This example is based on the one in the SciPy docs.
Imagine you have a signal composed of a 10 Hz signal and a 20 Hz signal. You want to remove the 10 Hz signal and leave the 20 Hz signal alone.
Normally, you would do that in one call, like this:
filtered = signal.sosfilt(sos, sig)
(where sig is the signal you want to filter, and sos is the second-order section filter found using scipy.signal.butter() with output='sos'. See the example in the SciPy docs for the full code.)
But you might want to do this without having the whole array in memory. Here is an example which does this in multiple steps, keeping the filter state between invocations.
sig_split = np.split(sig, 20)
zi = np.zeros((sos.shape[0], 2))
ret = []
for block in sig_split:
filt, zi = signal.sosfilt(sos, block, zi=zi)
ret.append(filt)
filtered = np.concatenate(ret)
You can see that you can do the same thing as the previous call, and that if you remove the zi argument passed to sosfilt, the output is no longer correct. (It looks like this if you don’t pass the zi argument.)
See also the sosfilt() documentation.
Most other functions in SciPy which evaluate a filter on a signal also accept a zi parameter, and you can do something similar with them.