Open scls19fr opened 6 months ago
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.
Thanks @SamProell for this update. If you are interested in Julia you might like https://github.com/JuliaDSP/DSP.jl/issues/548
Maybe some classes for stateful version of digital filters should be implemented
sure, this was just a quick demonstration of how to use the filter state. Writing a class around this would be necessary.
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)
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
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.
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.
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.
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
also many thanks to you!
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
https://github.com/SamProell/yarppg/blob/613155c935e9b1f0659deaf977938bcb105f4613/yarppg/rppg/filters.py#L15
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?