state in Filter

[David van Leeuwen]

Hello João, 

I am looking how I can replace my own implementation of an IIR filter with that in DSP.  Apart from differences in API that can be streamlined, I wonder how you keep the state of the filter in TFFilter.  

In my own implementation I go through some effort to store the filter state in the filter itself (and not as external data), so after

using SignalProcessing

x = randn(10)

f1 = Filter([1, -0.97])

a1 = x | f1

a2 = x | f1

a1 and a2 are a little different, because the filter state of f1 is different in the second call.  How do I do the same in DSP?

using DSP

f2 = TFFilter([1., -0.97], [1])

b1 = filt(f2, x)

b2 = filt(f2, x)

filt() can have a state as input, but how is the sate recovered after the first call?

Cheers, 

---david

[João Felipe Santos]

Hi David,

you can use filt!, which will mutate the variable passed as initial state. In our implementation, filters do not store state internally.

Cheers

--
João Felipe Santos

--
You received this message because you are subscribed to the Google Groups "julia-dsp" group.
To unsubscribe from this group and stop receiving emails from it, send an email to julia-dsp+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Simon Kornblith]

filt! makes a copy of the initial state before it does anything, so the si input isn't actually overwritten with the final state. But I agree that we need some way to get the final state.

We could change filt! to do this, but it's a little complicated because filt! can also do filtering along the first dimension of a matrix (or other nd array), and then it accepts either a matrix of initial states or a single initial state to use for all columns. In the latter case overwriting the initial state in place doesn't make so much sense.

We could also implement a streaming filtering API in DSP.jl instead of/in addition to hacking this onto filt!, although I think the basic representation of the filters should be separate from the state.

Simon

[David van Leeuwen]

Hi, 

On Thursday, October 9, 2014 2:10:16 AM UTC+2, Simon Kornblith wrote:

filt! makes a copy of the initial state before it does anything, so the si input isn't actually overwritten with the final state. But I agree that we need some way to get the final state.

yeah.  The idea for including the state in the filter instance is that 

y = filt(f, x)

should be exactly the same as

y = [filt(f, xx) for xx in x]

 

in the 1D case.  This saves memory for potentially very large arrays, and can also be applied online. 

We could change filt! to do this, but it's a little complicated because filt! can also do filtering along the first dimension of a matrix (or other nd array), and then it accepts either a matrix of initial states or a single initial state to use for all columns. In the latter case overwriting the initial state in place doesn't make so much sense.

In my own implementation, where time runs along the 1s dimension (just like in the case of DSP if I am not mistaken), I do this by copying the filter (without state) for all columns, but this is (probably) only limited to 2D.  I see I use mapslices(), so maybe this works in higher dimensions as well. 

 

We could also implement a streaming filtering API in DSP.jl instead of/in addition to hacking this onto filt!, although I think the basic representation of the filters should be separate from the state.

In analogy to an analog filter, the state is part of the filter.  Maybe you can define a new type that includes the state.

---david

[Jay Kickliter]

David,

This is exactly what I've been working on with FIR filters:

julia> myfilt = FIRFilter( h, 3//17 )
FIRFilter{FIRRational}(FIRRational(3x3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 1.0  1.0  1.0,3//17,3,3,0,1,1),[0.0,0.0],2)

julia> y1 = filt( myfilt, x[1:5] )
1-element Array{Float64,1}:
 1.0

julia> y2 = filt( myfilt, x[6:23] )
4-element Array{Float64,1}:
  6.0
 12.0
 18.0
 23.0

The main filter object, FIRFilter, holds the state and a FIRKernel object which can be FIRStandard, FIRDecimator, FIRInterpolator...

I don't use IIR filters that often, but I was thinking about doing something similar for them.

[Simon Kornblith]

Including the state in the filter is desirable for that specific case, but here are also some reasons not to:

• ZPK representation is not directly implementable without conversion, so ZPKFilters can't have a state.
• Transfer functions can be implemented in multiple ways with different states and numbers of states (although right now we only implement transposed direct form II).
• Applying the same filter to multiple different signals becomes more complicated.

Maybe the right way to go about this is to differentiate "filter representations" from "filter implementations." What DSP.jl has now are representations: they have no state and they need not be implementable directly. Filter implementations would carry both coefficients and state, or maybe even n states to handle streaming filtering of n channels.

Simon

[Simon Kornblith]

There's now a stateful IIR filter (DF2TFilter) in DSP.jl. Once we merge JayKickliter's Multirate work (#101), there will be stateful FIR filters as well.

Simon