Improve filtering (meta issue)

behinger commented 2 years ago

Widmann et al 2014 recomends not to use filtfilt, because cut-off frequency becomes dificult to estimate etc.

I also wouldn't use Butterworth but FIRWindow with a Hamming window; mirroring the firfilt plugin for eeglab from Andreas Widmann.

I'd be willing to re-implement that if this is of interest

Checklist

I (Rob) have hijacked the issue comment to list what's wanted and what's done

[x] Refactor filtering to use DSP.jl in a more intuitive way (expose filt and filtfilt)
[x] Refactor code so that EEG has a general filtering scheme, then sub types such as ASSR can have domain specific default parameters
[x] Rename XXXpass_filter to filte_XXXpass
[x] Add filterdelay function
[x] Add FIR filtering support
[ ] Ensure default EEG filtering is what @behinger suggests

rob-luke commented 2 years ago

Generally I think there are making too many absolute claims here (both in the existing code and your comment 😉). There are a broad variety of neuroscience applications where different analysis parameters are appropriate. For example, the use of Butterworth has been definetively demonstrated as appropriate for ASSR analysis. But it might be totally inappropriate for a different experimental paradigm as you suggest.

So I agree with you that we should overhaul the filtering system. But rather than changing the defaults to match EEGLab, we aim to make it more flexible (its currently too targeted to ASSRs, as this is my background). I think the dispatch system will make it easy to apply appropriate defaults for different data types. Let's throw around some API ideas and see what we both settle on. As a first suggestion how about...

We add a new function filter to compliment the existing highpass_filter, lowpass_filter, bandpass_filter.
We make the low level filter functions take a DSP filter object rather than setting defaults. We also ensure all functions that operate on NeuroImaging types (such as SSR or GeneralEEG) can also accept a DSP filter object. This provides maximum flexibility to the user, you can always specify your own filtering precisely.
The existing highpass_filter, lowpass_filter, bandpass_filter simply create appropriate DSP filter objects and pass this to the new filter function
We move the creation of the default filter generation to the specific data type. So the SSR data type would keep the current butterworth, but we change the GeneralEEG type to default to something else (I am happy with an FIR with Hamming window + delay compensation, I think this is what you suggest).
We add a filtfilt=True option to the filter functions that enables or disables filtfilt. If filtfilt=False then we will need to compensate for the delay of the filter.

So to an end user this would look something like...

# Maximum flexibility option
# 1. design the filter
using DSP
responsetype = Bandpass(300, 700; fs=48000)
designmethod = Butterworth(14)
zpg = digitalfilter(responsetype, designmethod)
fobj = DSP.Filters.DF2TFilter(zpg)
# 2. apply to EEG dat
s = read_EEG(data_path)
s = filter(s, fobj)  # Would apply a 14th order bandpass Butterworth to the GeneralEEG type using filtfilt

# Sane defaults for SSR type
s = read_SSR(data_path)
s = highpass_filter(s)  # Would apply the current filter settings to the SSR type
# But you could also apply the filter object from above the SSR if you wished...
s = filter(s, fobj, filtfilt=False) 

# Sane defaults for GeneralEEG type
s = read_EEG(data_path)
s = highpass_filter(s, l_freq=3u"Hz", filtfilt=False)  # Would apply zero phase fir hamming filter with lower pass-band edge of 3 Hz

Would something like the above work for you?

recomends not to use filtfilt, because cut-off frequency becomes dificult to estimate etc.

Acausal filters are commonly used in a huge variety of neuroscience and has great benefit. I agree the output needs to be interpreted carefully, but the same can be said for a single pass filter with no phase compensation. Every step of processing has implications to interpretation. I also prefer the discussion in https://www.cell.com/neuron/fulltext/S0896-6273(19)30174-6 Generally I think filtfilt is appropriate for certain applications, but I agree not all. So lets just add an argument so the user can use it or not.

mirroring the firfilt plugin for eeglab from Andreas Widmann.

I have never used eeglab and would prefer to base decisions on practical considerations rather than because another package does it. I am happy to make the change you suggest, but let's base decisions on the implications to analysis rather than what other libraries do.

I'd be willing to re-implement that if this is of interest

Absolutely of interest! And thank you for the offer to implement!

But I think I have opened a bigger can of worms that you originally suggested 😉. I greatly appreciate your offer of help here, but I am aware that the above suggestions may be more than you originally offered to implement. Are you willing to try and implement something like I suggest above (or whatever API we finally decide on together)? If this is too much to ask you, then we could split it up, I could do the refactoring and then you could set the new defaults. Just let me know your preference.

behinger commented 2 years ago

Hi! your proposal is good. Thanks for the extended answer!

I think I failed to see the scope of this package - of course in fMRI you'd have quite different filters to EEG. I mostly know about filters in EEG for ERPs and I think Widmann 2014 is currently the definite guide for that. But imho MNE use butterworth as default

Regarding filtfilt, the trouble is that your filterorder is not the order you specify. And most people would give the filterorder and/or cutoff, which is not the actual cutoff if you filter twice. But reading your answer, you seem to be aware of it. There is some discussion in Widmann 2014 on this.

Exactly as you write I dont want causal filters (for now), they o.c. have their own problems, but the linear phase filter + delay compensation

Defaults go to lowpass_filter, which internal build the filter and call filter. Quick q, do you think filter_lowpass could be the more discoverable term?

I think the implementation is not that complex, I will give it a shot!

rob-luke commented 2 years ago

And most people would give the filterorder and/or cutoff, which is not the actual cutoff if you filter twice.

Your concern seems to be with people reporting their analysis incorrectly in journal articles. I think that's beyond the scope here (that's why we have supervisors and reviewers etc), I would welcome a function that describes the analysis in plain English if you think that would help.

do you think filter_lowpass could be the more discoverable term

Agree, feel free to change them all to filter_xxxx

I think the implementation is not that complex, I will give it a shot!

Wonderful thanks. I will assign you to the issue then.

Not sure if this is documented anywhere. But the functions that operate on basic types like arrays or matrices or DSP objects etc live in the base level of src (https://github.com/rob-luke/Neuroimaging.jl/blob/main/src/preprocessing/filtering.jl), whereas the type specific wrapping code would live in the type directory (https://github.com/rob-luke/Neuroimaging.jl/blob/main/src/types/SSR/Preprocessing.jl).

rob-luke commented 2 years ago

Could you also please add a brief documentation page that describes an overview of the approach here? Filtering is pretty fundamental question for everyone who will use this package, so I think even a few sentences would be a great benefit. Feel free to add references etc. Somewhere around here: https://github.com/rob-luke/Neuroimaging.jl/tree/main/docs/src

rob-luke / Neuroimaging.jl

Improve filtering (meta issue) #164

Checklist