filtering evoked and epochs

[jjnurminen]

Hi, there is no facility for filtering Epochs or Evoked, except by Savitzky-Golay polynomial fit. Apparently this is to encourage filtering at raw data level instead (to avoid edge artifacts). Couple of points/questions:

• Why Savitzky-Golay? As Ronald Schafer has pointed out [1], S-G is equivalent to convolution with a finite impulse response, but actually has properties inferior to standard FIR design methods.
• Why are edge artifacts so important? If significant ones occur with a given filtering method, they are limited in duration and in most cases easily ignored (IMO).

Filtering raw data is very slow in many cases, so it would be nice to have .filter() method to quickly play around with frequency domain properties of Evoked and Epochs. For example, the Elekta software allows this. Of course one has to be careful not to distort data with e.g. carelessly applied highpass filters, but the same applies equally to filtering Raw.

[1] https://ai.berkeley.edu/~ee123/sp15/docs/SGFilter.pdf

[agramfort]

the motivation for SG was that it tends to keep peak amplitudes better. It's quite important for evoked data with ERF/ERPs.

how would you avoid users making mistakes with edge artifacts?

can you mimick what Elekta software does. Doing what Elekta does seems like a good plan :)

[jjnurminen]

the motivation for SG was that it tends to keep peak amplitudes better.

But this is equivalent to having a maximally flat passband response. As shown in the paper I cited, an equally long Parks-McClellan FIR filter has even flatter passband and better stopband attenuation too. (Or you could use e.g. Butterworth IIR). So AFAIK, S-G does not have any advantage in this respect.

how would you avoid users making mistakes with edge artifacts?

Is the S-G implementation free of edge artifacts then? At least the scipy docs talk about different padding options. IMO, edge artifacts are a fact of life when doing filtering and users should just be aware of them. On the other hand, having tested both S-G and scipy.filtfilt on my evoked data, both produced minimal if any edge artifacts.

can you mimick what Elekta software does. Doing what Elekta does seems like a good plan :)

Well I'm not 100% sure about that :) However it's very convenient for playing with evoked data. For example, today I was playing with different highpass/lowpass corners for the HF-SEF data. In mne-python this is pretty slow and clumsy, as you need to reload, re-filter and re-epoch big chunks of raw data. Instead I ended up just applying filtering directly on evoked.data; this is pretty simple to do, but a ready-cooked method would be even better.

[jjnurminen]

To clarify my actual point: if the API is going to offer users a possibility to filter Evoked or Epochs (and I don't see why it shouldn't), a proper FIR or IIR filter would be preferable to savgol_filter.

[agramfort]

what you recommend in terms of code / API change?

[jjnurminen]

Well as for my 5 eurocents, I would just put a .filter() method in there and default to some reasonable FIR or IIR filter.

[agramfort]

I have no objection but I'd like to see experimental evidence that the defaults we provide for such short signals with ERF and ERPs are at least as good or better than the SG filter we provide currently.

can you do a quick demo / proof of concept?

[jjnurminen]

Ok, here's some sandboxing. I filtered HF-SEF evoked data with 5th order Butterworth IIR (forward-backward) and Savitzky-Golay. I specified 100 Hz as the lowpass corner for both. Filtered data below:

Savitzky-Golay:

Butterworth:

The edge effects are minimal (and similar for both filters). Getting the peak response amplitudes gives this:

original: 6.20207e-12 T/m S-G filtered: 5.33125e-12 T/m Butterworth filtered: 4.88476e-12 T/m

So indeed it looks like Butterworth is reducing the response amplitude a bit more. But why is this?

Plotting the filter frequency responses, we get the following:

The Butterworth lowpass corner, i.e. -3 dB point was specified at 100 Hz. (This becomes -6 dB in the freq response because of forward-backward filtering; the filter is applied twice.) However we see that S-G barely does anything at 100 Hz; its -6 dB point is around 187.5 Hz. (You can see this by eye in the filtered data: the S-G filtered data clearly has more high frequency content)

To design more comparable filters, we can put the Butterworth -6 dB point also at 187.5 Hz. Then the peak amplitudes become:

original: 6.20207e-12 T/m S-G filtered: 5.33125e-12 T/m Butterworth filtered: 5.38546e-12 T/m

Now the results are very similar. In other words, the "peak amplitude preserving" feature of S-G is due to it being a bad filter! Some of the energy of the response peak is in high frequencies; S-G does not reduce that very efficiently, and thus more of the peak amplitude is preserved. The same is pointed out by Schafer in the above article:

"S-G filters are often preferred (even revered in some circles) because, when they are appropriately designed to match the waveform of an oversampled signal corrupted by noise, they tend to preserve the width and height of peaks in the signal waveform. While such performance features are often explained in terms of matching fitted polynomial slopes to signal slopes or to the preservation of signal moments, the reason for this behavior is more obvious from the frequency-domain properties of the filters. Specifically, S-G filters have extremely flat passbands with modest attenuation in their stopbands..."

So my personal conclusion would be:

• it's better to use some appropriate IIR or FIR filter instead of Savitzky-Golay: you'll have more control of the filter properties and better performance

• loss of response amplitude is a natural consequence of suppressing high frequencies; don't filter your data more aggressively than is necessary

[agramfort]

I'll try to look soon. Someone else please pitch in.

[larsoner]

I'm okay with adding it in principle. But to do it right will take more than the code changes, which would be something like a half-dozen-line Mixin wrapper to filter_data, because we need to try to stop people from making mistakes (which is perhaps easier to do when filtering is done at the evoked stage rather than raw stage). Specifically I think we'd need to:

1. Add to the filtering background tutorial to talk about edge artifacts, how they can arise, and how they can be avoided or minimized. I think this is necessary both from a pedagogical standpoint, but more critically also to make sure we've explored all of the potential consequences so we can avoid letting people do something bad with what looks like a reasonable operation.

2. From explorations in (1) we should ensure that our warnings about signal vs filter lengths are correct.

3. We might want to expose a pad_type argument -- right now IIRC we use "odd" but people might want other choices. (This could be a subsequent PR depending on whether or not our padding creates problems in (1).)

4. Document the differences in processes (maybe in a tutorial?). Going raw->filter->epochs->evoked will be different from raw->epochs->evoked->filter in a number of potentially important ways:

   1. peak-to-peak trial rejection during epoching
   2. epochs decimation capabilities
   3. introducing differences in frequency content during covariance estimation (if Epochs are unfiltered) vs in evoked data (if Evoked ends up being filtered differently)

Point (iii) might not be a problem if it's assumed that the spatial pattern is conserved across frequency regions (we probably make this assumption already in other places), but we should point it out nonetheless.

Based on what you've shown / paper you've linked to @jjnurminen I'd actually be +1 for deprecating SG filtering after having this, ideally with a nice addendum to the tutorial showing an example of inferior stop-band performance.

[larsoner]

@jjnurminen if you're up for doing some or all of this that would be awesome. It will be time consuming, though, so I'm also happy to help or do it when I have time, but it might take me a while to get to it.

[jjnurminen]

Thanks for the comments @Eric89GXL I would be up for doing some of it, sure. I can start to look into the necessary things - it will take me a while too. (Just as an example, I have zero idea what a 'Mixin wrapper' is :)

[larsoner]

That would be this. Basically you'd want to take the code from mne/io/base.py for BaseRaw.filter and move it to FilterMixen.filter. The FilterMixin is inherited by Epochs and Evoked, and you'd want to make it also inherited by BaseRaw.

[larsoner]

@jjnurminen have you started on this? If not I can probably do it before 0.15