Stream data to digital filter / live IIR filter

[s-celles]

Hello,

Following discussion started at https://discourse.julialang.org/t/scipy-signal-iirfilter-equivalent-in-dsp-jl-and-more/110783/4 I wonder if DSP.jl provides structure of digital filter which can be used for streaming input data ie processing live data and not a vector of precalculated values ?
ie is there some kind of filter which preserve its state each time a new input is coming.

I did this Pluto.jl notebook https://gist.github.com/scls19fr/0ae16d92ca39d3eb9c42cc0fc618c723 for experimenting this kind of problem (both for filtering a vector of input signal and for displaying streamed signal and output of filter thanks to PlutoUI Clock).

I end up to convert idea mentioned https://www.samproell.io/posts/yarppg/yarppg-live-digital-filter/ with Julia.

OnlineStatsBase.jl API looks interesting for this kind of work (this is API which is used by OnlineStats.jl

So I did

abstract type AbstractLiveFilter{T} <: OnlineStat{T} end

mutable struct LiveDigitalFilter{T} <: AbstractLiveFilter{T}
    value::T
    n::Int

    b::Vector
    a::Vector

    _xs::CircularBuffer
    _ys::CircularBuffer

    function LiveDigitalFilter{T}(pr::PolynomialRatio) where {T}
        b = reverse(pr.b.coeffs)
        a = reverse(pr.a.coeffs)
        _xs = CircularBuffer{T}(length(b))
        fill!(_xs, zero(T))
        _ys = CircularBuffer{T}(length(a) - 1)
        fill!(_ys, zero(T))
        new{T}(0, 0, b, a, _xs, _ys)
    end
end

function reset!(f::LiveDigitalFilter)
    f.n = 0
    empty!(f._xs)
    empty!(f._ys)
    fill!(f._xs, zero(eltype(f.b)))
    fill!(f._ys, zero(eltype(f.a)))	
end

function OnlineStatsBase._fit!(f::LiveDigitalFilter, x)
    f.n += 1
    pushfirst!(f._xs, x)
    y = sum(f.b .* f._xs) - sum(f.a[2:end] .* f._ys)
    y = y / f.a[1]
    pushfirst!(f._ys, y)		
    f.value = y
end

Notebook run on my side at around 10fps.

Any opinion?

[mbaz]

Options I see: use FIRFilter, which keeps its state between calls; or use a version of filt! that allows setting the filter state si.

The documentation here is quite unhelpful, unfortunately: for example, the two definitions of filt here are ambiguous.

[wheeheee]

Yes, in addition to what @mbaz mentioned, I believe filt with DF2TFilter (or filt!, plus a 1-element out buffer) does what _fit! does, although it should be much better to buffer x into a vector and filt! that.

julia> p = PolynomialRatio([1:3;], [2:10;]);

julia> a = LiveDigitalFilter{Float64}(p)
LiveDigitalFilter: n=0 | value=0.0

julia> _fit!(a, 10), _fit!(a, 23), _fit!(a, -1)
(5.0, 14.0, 6.5)

julia> _fit!(a, 2), _fit!(a, 3), _fit!(a, 4)
(-15.75, -37.375, 19.8125)

julia> sf = DF2TFilter(p);

julia> filt(sf, [10,23,-1])
3-element Vector{Float64}:
  5.0
 14.0
  6.5

julia> filt(sf, [2,3,4])
3-element Vector{Float64}:
 -15.75
 -37.375
  19.8125

[wheeheee]

The documentation here is quite unhelpful, unfortunately: for example, the two definitions of filt here are ambiguous.

I think type annotations within the docstring function definitions would be helpful?
Also, at least for DF2TFilter, there's this: https://docs.juliadsp.org/dev/filters/#DSP.filt

[s-celles]

Thanks @wheeheee and @mbaz

@wheeheee I'd do

p = PolynomialRatio([1:3;], [2:10;])
ldf1 = LiveDigitalFilter{Float64}(p)
fit!(ldf1, Float64[10, 23, -1])
fit!(ldf1, Float64[2, 3, 4])
println(value(ldf1))  # or println(ldf1)

What is odd with filt! or filt is that you have to create a 1 element array to process just one value. Isn't there a function which only takes one value?

The main advantage I see with working with OnlineStatsBase is that our filter can easily integrates with all online algorithm statistics that https://joshday.github.io/OnlineStats.jl/ provides

Maybe something like that

mutable struct LiveDigitalFilter2{T} <: AbstractLiveFilter{T}
    value::T
    n::Int
    coeff

    function LiveDigitalFilter2{T}(coeff) where {T}
        new{T}(0, 0, coeff)
    end
end

function OnlineStatsBase._fit!(f::LiveDigitalFilter2, x)
    f.n += 1
    f.value = filt(f.coeff, [x])[end]
end

(not sure if coeff is the appropriate name field... not sure what type should be used also)

and we can use it like

ldf2 = LiveDigitalFilter2{Float64}(DF2TFilter(p))
fit!(ldf2, Float64[10, 23, -1])
fit!(ldf2, Float64[2, 3, 4])
println(value(ldf2))  # or println(ldf2)

or with filt! calls inside fit!

mutable struct LiveDigitalFilter2{T} <: AbstractLiveFilter{T}
    value::T
    n::Int
    
    coeff
    out::Vector

    function LiveDigitalFilter2{T}(coeff) where {T}
        new{T}(0, 0, coeff, T[0])
    end
end

function OnlineStatsBase._fit!(f::LiveDigitalFilter2, x)
    f.n += 1
    filt!(f.out, f.coeff, [x])
    f.value = f.out[end]
end

This might be useful. If user wants to preserve older output values, b can be wrapped into StatLag

You are right. Type annotions would be helpful.

Moreover I think doc would be much better with some examples, some plots, some notebooks... it will help beginners to better understand how to start with this library.

[mbaz]

The documentation here is quite unhelpful, unfortunately: for example, the two definitions of filt here are ambiguous.

I think type annotations within the docstring function definitions would be helpful? Also, at least for DF2TFilter, there's this: https://docs.juliadsp.org/dev/filters/#DSP.filt

Yes, type annotations is what is needed.