Live digital filters as a separate package / library

[scls19fr]

Hello,

I read with interest your blog post https://www.samproell.io/posts/yarppg/yarppg-live-digital-filter/

I wonder if a separate package just for implementing live digital filters in Python won't make sense.

It could be based on list as you did in

yarppg/yarppg/rppg/filters.py

Line 15 in 613155c

self._xs = [0]*len(b)

or on deque as mentioned in your blog post.

A last idea could be to rely on circular buffer / ring buffer on top of Numpy such as https://github.com/eric-wieser/numpy_ringbuffer or https://github.com/Dennis-van-Gils/python-dvg-ringbuffer

Any opinion?

[SamProell]

Hi thank you for the feedback.
Sounds like an interesting idea. And I did not find any similar library in a quick search..
I have also learned more about this topic since I wrote the post and I would suggest a slightly less complicated solution using scipy's filter state (zi). One could do something like below, without the need to track past inputs and outputs, and solely relying on the scipy implementation:

zi = scipy.signal.lfilter_zi(b, a)
ys = []
for y in yraw:
    yout, zi = scipy.signal.lfilter(b, a, [y], zi=zi)
    ys.append(yout)

I intend to make an update to the post once I get around to it.

[scls19fr]

Thanks @SamProell for this update.
If you are interested in Julia you might like JuliaDSP/DSP.jl#548

[scls19fr]

Maybe some classes for stateful version of digital filters should be implemented

[SamProell]

sure, this was just a quick demonstration of how to use the filter state. Writing a class around this would be necessary.

[scls19fr]

I tried

class StateFulLiveLFilter(LiveFilter):
    def __init__(self, b, a):
        """Initialize live filter based on difference equation.

        Args:
            b (array-like): numerator coefficients obtained from scipy.
            a (array-like): denominator coefficients obtained from scipy.
        """
        self.b = b
        self.a = a
        self.zi = scipy.signal.lfilter_zi(self.b, self.a)

    def _process(self, y):
        """Filter incoming data with standard difference equations.
        """
        yout, self.zi = scipy.signal.lfilter(b, a, [y], zi=self.zi)
        return yout

but MAE is a bit higher

lfilter error: 1.8117e-15 (with previous version)
lfilter error: 0.039018 (with this version)

[scls19fr]

Same problem also occurs with

class StatefulLiveSosFilter(LiveFilter):
    """Live implementation of digital filter with second-order sections.
    """
    def __init__(self, sos):
        """Initialize live second-order sections filter.

        Args:
            sos (array-like): second-order sections obtained from scipy
                filter design (with output="sos").
        """
        self.sos = sos
        self.zi = scipy.signal.sosfilt_zi(sos)

    def _process(self, y):
        """Filter incoming data with cascaded second-order sections.
        """
        yout, self.zi = scipy.signal.sosfilt(self.sos, [y], zi=self.zi)
        return yout

sosfilter error was 0 and is now 0.039018

[SamProell]

In principal, your implementation(s) are good. I played around with it myself and the "problem" lies with setting the initial state.
You should be able to see, that the results do not match only in the beginning of the signal.
Apparently, to reproduce the behavior of lfilter(b, a, y), zi has to be initialized differently. The documentation for lfilter states

If zi is None or is not given then initial rest is assumed. See lfiltic for more information.

After some playing around, I found that zi = scipy.signal.lfiltic(b, a, [0]) does the trick. At least for lfilter, I cannot find something similar for sosfilt.

[scls19fr]

Ok so StateFulLiveLFilter should probably be

class StateFulLiveLFilter(LiveFilter):
    def __init__(self, b, a, yi=None, xi=None):
        """Initialize live filter based on difference equation.

        Args:
            b (array-like): numerator coefficients obtained from scipy.
            a (array-like): denominator coefficients obtained from scipy.
        """
        self.b = b
        self.a = a
        if yi is None:
            yi = [0.0]
        self.zi = scipy.signal.lfiltic(b, a, yi, xi)

    def _process(self, y):
        """Filter incoming data with standard difference equations.
        """
        yout, self.zi = scipy.signal.lfilter(b, a, [y], zi=self.zi)
        return yout

no idea about how to overcome this sosfilt error.

[SamProell]

looks good. I started working on a "minor" rework of the internal structure, including formatting. If it's ok with you, I might just replace the current implementation of DigitalFilter with this one. Or you can wait and then make a PR to get the full credit ;)

I'd say that having some deviation in the beginning is not that big of a deal for most applications. Focus is on online filtering, where some startup differences should be no problem. And also, you cannot get the same behavior with the original sosfilt anyway, i guess.

[scls19fr]

Yes we can live with this initialization problem... but maybe that's a lack on Scipy side
I don't really care about the "full credit" 😄 Thanks @SamProell

[SamProell]

also many thanks to you!