BUG/API with filtering

[dengemann]

@jdammers just demonstrated to me how to silently produce nans with our band pass filter.

https://gist.github.com/dengemann/7580547

Or old non-Python routines default to order 16 butterworth band pass filters (FFT implementation) ... If you go for lower orders, in certain application evil filter artifacts would occur.

I think we should

1) add a clear RuntimeError if NaNs are returned by scipy.signal.filtfilt

2) @jdammers suggested to port his IDL fft code, currently we don't support order type construction for fft filters.

WDYT?

[jdammers]

is anybody interest in converting this https://gist.github.com/jdammers/7580737 to Python (FFT filtering using butterworth) ? The code is available in IDL only at the moment.

[agramfort]

can you clarify where the problem is with python? the filter coefs are given from scipy than we rely on filtfilt in scipy.signal. So if there is a problem in scipy it should be fixed there and back port it.

my 2c

[dengemann]

@agramfort I think we don't have equivalent filters in scipy ... the point is that currently we ignore the order and ftype arguments if method == 'fft'. The filter coefs seem ok, albeit small, utlimately filtfilt fails, for our narrow bandpass use case already with order > 4. My point is that we should at least catch nans that come from filfilt + we should support higher orders via FFT.

[agramfort]

@agramfort I think we don't have equivalent filters in scipy ... the point is that currently we ignore the order and ftype arguments if method == 'fft'.

yes fft based filters as in MNE-C are different. It's the transition that only matters.

The filter coefs seem ok, albeit small, utlimately filtfilt fails, for our narrow bandpass use case already with order > 4. My point is that we should at least catch nans that come from filfilt + we should support higher orders via FFT.

does it means the FFTs used in filtfilt are too small?

[dengemann]

does it means the FFTs used in filtfilt are too small?

good question ....

Btw. it seems this butterworth recipe is working well at high orders, additional tests with filtfilt needed ... http://stackoverflow.com/questions/12093594/how-to-implement-band-pass-butterworth-filter-with-scipy-signal-butter

[mluessi]

@dengemann do you have an example use-case where the filter gives NaN's? I think we should try to fix it first and add different types of filters later.

[dengemann]

@mluessi see https://gist.github.com/dengemann/7580547

[mluessi]

@dengemann my assumption is that the transition bandwidth is too small and the filter was very high order and becomes unstable due to numerical accuracy. We should add a check for filter stability (poles within unit circle).

[dengemann]

was very high order and becomes unstable due to numerical accuracy. We should add a check for filter stability (poles within unit circle).

I don't know how 'very high' is actually to be construed, but filtering at order 8 should not be something too exciting. Larger trans widths do not seem to change the behavior. Check sounds good.

[jdammers]

@mluessi @agramfort @dengemann although we use order 16 only for phase analysis in our IDL filter implementation the FFT filter has no problems with filtering the (same) data using higher orders. However, our student wanted to filter MEG data using MNE and was playing with the filter order, with the result that the filter generated always NaNs when order was set larger than 4.

[larsoner]

The code chunk @dengemann posted will indeed just wrap to scipy since it's IIR filtering, and doesn't use any FFT-based filtering (which is only used for FIR filtering). Similarly, the trans widths are ignored when doing IIR filtering when the order is supplied -- they are, however, used to calculate the order if it is not supplied. So it is not surprising @dengemann that you changing the trans_bandwidth had no effect when having an order=8 in there. See:

https://github.com/mne-tools/mne-python/blob/master/mne/filter.py#L402

Like @agramfort said, if there is an issue when using IIR mode, it is effectively a problem with scipy (and by association, our use of it). The sanest approach to me would be to 1) add a check on our end that works, to fix this problem ASAP, and 2) port it over to scipy so that other people don't hit this issue. Basically I think it will boil down to filter stability, and how to establish it. And it's especially important since we do calculate the order, if it's not supplied, for IIR filters using the trans_bandwidth. So we'll have to figure out how to do the check... a post-hoc check for nan isn't pretty, but it could work in the short-term if we can't find a more reliable solution.

@rkmaddox knows more about this than me, but I think that such high-order filters using BA-filtering like scipy uses for flitfilt will sometimes (often?) not be stable. To get this working, we'd have to implement ZPK filtering, which should really also go in scipy since they already provide options for producing coefficients in ZPK instead of BA form. Maybe if somebody feels like doing a really difficult Sunday hack :) I might email scipy to make sure they don't have it implemented and I'm missing it, given that you can already create the appropriate coefficients.

Thinking about it, I'm curious what the motivation is for using high-order Butterworth? IIR butterworth should have a flat pass-band (assuming numerical issues don't crop up), but truncating the impulse response to obtain a FIR approximation should introduce ripple (albeit small)... in which case there could be other transforms that might work better or at least as well? Just throwing out ideas.

[larsoner]

FYI apparentnly there is some progess on this issue in scipy:

https://github.com/scipy/scipy/pull/3085

[dengemann]

@Eric89GXL I agree with you. It seems we're discussing separate issues. 1) deal with the scipy / iir issue. 2) implement a working FFT alternative (here/ scipy?). 3) an API question, related to many parameters effictively not passed / used and as @jdammers pointed out to me: no option to chose a filtertype for the FFT route. I hope this admittedly dense summary makes sense. I would classify your comment as a proposal on how to deal with 1)

[larsoner]

Yeah, it does just deal with 1). @agramfort @mluessi WDYT about adding a check on our end until scipy gets ZPK / second-order-sections filtering?

In terms of 2). My hesitation with including the Butterworth FIR approximation code is that, where do we stop? If we include that, it's weird not to also include the many other standard filter types that are available in e.g., scipy.signal. And if we do include all filter types because we have the code for it, then they actually seem to be a better fit for going into scipy. We stopped with what scipy included because that's a standard, tested (although apparently here not quite tested enough) package.

In terms of 3). We do have our own "filter design" in the sense that we have mandated that FIR filters are trapezoids in frequency, essentially. We could do things like allow people to use causal filters instead of forward-backward filters, or allow for different shapes. But again, this becomes a "where do we stop" problem. In terms of certain parameters not being passed, we can update our docstring, or change the API. The latter is a pain, but might be better long-term.

What do you think about exposing an API for using arbitrary frequency-domain kernels (H)? That would allow your lab to implement any FIR filters you want, effectively solving 2) and most of 3). An additional question is whether you want support for causal (only forward) filtering, or if you're happy with our current use of forward-backward filtering. It would make sense to me to expose both options potentially. We'd have to think about the API, but it would be doable.

[dengemann]

Yo yo. Hello folks,

some updates: https://github.com/scipy/scipy/issues/2980

the a, b representation of a butterworth kernel lends itself towards instability. But rescaling can help. Tested it and higher orders fly with our butterworth once you scale the data to Tesla, Volt, etc.

While @Eric89GXL is merging https://github.com/scipy/scipy/pull/3717 for scipy 0.15 I suggest to version check scipy and apply the rescaling trick if the filter crashes (coefs outside unit circle).

cc @jdsitt @agramfort @jdammers @fboers

[larsoner]

We have a check for poles outside the unit circle already, closing this one