JuliaDSP / DSP.jl

Filter design, periodograms, window functions, and other digital signal processing functionality
https://docs.juliadsp.org/stable/contents/
Other
380 stars 110 forks source link

Stateful filtering fails with multi-dim array #372

Open hergum opened 4 years ago

hergum commented 4 years ago

Using more than one dimension in the input array makes filt fail for stateful filter.

s = rand(10,3) df = digitalfilter(Lowpass(0.25), Butterworth(4)) f = DF2TFilter(df) filt(f,s)

ERROR: MethodError: no method matching filt(::DF2TFilter{SecondOrderSections{Float64,Float64},Array{Float64,2}}, ::Array{Float64,2})

One-dimensional input works as expected. Two-dimensional inputs works along the first dimension, as expected, for stateless filtering filt(df,s).

Any ideas on how to make this work according to the doc?

galenlynch commented 4 years ago

You're right that the docs suggest that this should work. However, that doc string was written for normal filter coefficients, and not stateful filters. It's not really clear to me what the expected behavior would be when a stateful filter is applied to a matrix. For stateless filtering each column is treated as a separate signal, where the filter state is reset at the start of each column of the signal, and the filter state is discarded at the end of each column. If you were to do the same with a stateful filter, it's not clear what the filter state should be either at the start of each column, or after all the columns have been filtered (i.e. which column of the signal should the state reflect).

Depending on what you are expecting the filter state to do between columns of your signal, you could write an explicit loop accordingly. The only thing that really makes sense to me, given your decision to use a stateful filter, is that each column is a continuation of the same signal. In this case, there is no reason to use a matrix instead of a vector. If you instead want to maintain separate filter state for each column, then you should make separate filter objects for each column.

Could you explain your use case a little bit more? Are you processing very long signals in batches, or are you doing some on-line filtering of signals? Are you filtering multiple channels of data that should be treated separately? What is the reason that you want to use stateful filters in the first place?

hergum commented 4 years ago

I would like the stateful filter to keep track of the state of the filter for each column.

My use-case is I get long 2D signals in batches. Each row is independent, and all rows should be processed the same way. The rows are continuous in time, and must be processed as if I had the full signal available, so I need to preserve the filter state. Its for a medical device, but the use case might as well be wanting to low-pass filter the audio track of 15 radio stations. In matlab I'm used to doing things like:

>> [b,a]=butter(4,0.25);
>> [s,state]=filter(b,a,rand(15,800),[],2); % [] means no state for the first batch
>> size(state)

ans =

     4    15
>> [s,state]=filter(b,a,rand(15,800),state,2);  % use the filter state from the end of last batch to initialize the next

The documentation suggested I could do this very conveniently in Julia with a stateful filter. I'll do it in a loop instead.

galenlynch commented 4 years ago

I see, that makes sense. Unfortunately the design of stateful filters in this package requires that they pre-allocate the state associated with the filter. This makes them rigid, and incapable of handling either changing the type of input that the can handle (your last issue), or the number of inputs (this issue). It seems like the Matlab approach of simply returning the state along with the filtered result is simple, and also more flexible. The down side of doing something similar in Julia is that functions can't change their number of outputs as easily as in Matlab, meaning that if we were to always return the state, then everyone using filt would have to handle the state output. This might still be worth the trade off, both because it would simplify the code, and eliminate both of the valid issues that you've raised.

That being said, I've always found filt treating multi-dimensional arrays as many signals, instead of only operating on vectors, to be a holdover from Matlab, and haven't been sure it's a great idea in the first place.

On Mon, Aug 17, 2020, 4:49 PM hergum notifications@github.com wrote:

I would like the stateful filter to keep track of the state of the filter for each column.

My use-case is I get long 2D signals in batches. Each row is independent, and all rows should be processed the same way. The rows are continuous in time, and must be processed as if I had the full signal available, so I need to preserve the filter state. Its for a medical device, but the use case might as well be wanting to low-pass filter the audio track of 15 radio stations. In matlab I'm used to doing things like:

[b,a]=butter(4,0.25); [s,state]=filter(b,a,rand(15,800),[],2); % [] means no state for the first batch size(state)

ans =

 4    15

[s,state]=filter(b,a,rand(15,800),state,2); % use the filter state from the end of last batch to initialize the next

The documentation suggested I could do this very conveniently in Julia with a stateful filter. I'll do it in a loop instead.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaDSP/DSP.jl/issues/372#issuecomment-675106967, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABUDZTYOQV6I53TKM6T3NZ3SBGJULANCNFSM4QBN2EKQ .

martinholters commented 14 hours ago

Adding this to the milestone for 0.8 so that we don't forget to consider any breaking changes (like adding the state as a second return value) required for fixing this and whether we can make them land in time for 0.8. If we cannot work out what the solution should be, we can always postpone to 0.9 (or any 0.8.x if we figure out a non-breaking way).