Live digital filters as a separate package / library about yarppg HOT 11 OPEN

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.

from yarppg.

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

from yarppg.

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

from yarppg.

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

from yarppg.

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)

from yarppg.

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

from yarppg.

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.

from yarppg.

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.

from yarppg.

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.

from yarppg.

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

from yarppg.

also many thanks to you!

from yarppg.