DOC: improved some parts of the filtering and resampling tutorial

[robertoostenveld]

Reference issue

This is my review of #6872

Remaining issue

There is one section that needs further attention: it writes "and the data resampled from 600 Hz to 200 Hz gets automatically lowpass filtered at 100 Hz"

Should it not automatically be filtered at 67 Hz?

The subsequent figure also does not show the effect of the anti-aliasing filter, which suggests that the anti-aliasing filter was indeed at 100 Hz or higher (or absent). That is not sufficient, it should be at 1/3*200=67.

[cbrnr]

The theoretical cutoff frequency should be 100Hz - which is 1/3 times the original Nyquist frequency, which is 600/2=300Hz. Therefore, the cutoff frequency should be 1/3*300Hz=100Hz.

However, the practical cutoff frequency should probably be chosen lower than 100Hz due to non-ideal filter characteristics. Is that what you meant by setting it to 67Hz?

[robertoostenveld]

When (re)sampling, the original sampling frequency does not matter for the aliasing filter, the target sampling frequency does: the cutoff frequency of an anti aliasing filter should be as a rule-of-thumb be 1/3 or 1/4 of the desired sampling rate. Fir digitization at a certain frequency, you should ensure that the to-be-digitized data does not contain frequencies above the Nyquist frequency. If the filter (as in reality) does not have an infinitely steep drop-off, the anti-aliasing filter should therefore be below the Nyquist frequency. This is correctly explained in the text, which states " best practice is to low-pass filter the Raw data at or below 1/3 of the desired sample rate".

Resampling is not different than sampling. If the original data were analog (i.e. not yet sampled), its frequency content in principle goes to infinite. In that case you would also not relate the anti-aliasing filter to infinite, or 1/2 times infinite.

In this case the original data sampled at 600Hz had an anti aliasing filter at 172, which is a factor 172/600=1/3.5. When decimating this data, you could also decide to keep the filter at 1/3.5 of the sampling frequency, which means that the new anti-aliasing filter ends up at 57Hz.

[cbrnr]

Like I said, it's a rule of thumb. The only theoretical number you can state is that the anti-aliasing filter must be at the target Nyquist frequency, which is 100Hz in this example.

[cbrnr]

So what I mean is that what the resample function does is the theoretical bare minimum (100Hz). If you want to change this, the function would need to be adapted.

[robertoostenveld]

We discussed yesterday (in person with @agramfort and @drammock) that as rule of thumb the filter should be at 1/3 of the new sampling frequency and that its should be correspondingly documented. If the documentation is not consistent with the code (as it now seems), then either one should be fixed.

[cbrnr]

Yes, I agree. Maybe there could be an option where users can supply the desired cutoff frequency (as a multiplier in the interval (0, 0.5] where 0.5 should remain the default for backward compatibility). I just wanted to avoid confusion between theoretical limits and practical rule-of-thumb values.

[larsoner]

If the documentation is not consistent with the code (as it now seems), then either one should be fixed.

The documentation probably needs to be clarified here. I think in principle things are consistent -- there is as a difference between what we know when we construct and apply the antialiasing filter and then decimate (i.e., resample) versus what we know when users apply their own lowpass then decimate.

• When users decimate, we don't know their filter characteristics, so we use the 1/3 rule of thumb to check.
• When we do a resampling operation, we actually currently use a frequency-domain method, as opposed to e.g., a polyphase (MRG: Add polyphase resampling #5136) method. The AA filter here is essentially a brick-wall frequency-domain filter. This operation is probably undesirable compared to polyphase filtering for other reasons (circularity of the effective convolution for example), but aliasing might not be the problem with it.

We could resurrect #5136, but even in that case, rather than following the 1/3 rule of thumb when constructing our own antialiasing filter, it seems like it would make more sense to intentionally construct the FIR filter using known pass-/transition-/stop-band characteristics. Then we'd know how it behaves, and wouldn't have to follow the rule of thumb because the antialiasing will depend just on the stop/transition band (and thus effective FIR order) choices.

So I don't mind having a potential mismatch between what we do in .resample and what we check for in epochs.decimate. Expert users who have constructed a very sharp filter, for example, can do verbose='error' to ignore our -- in their case -- too conservative warning.

[cbrnr]

When users decimate, we don't know their filter characteristics, so we use the 1/3 rule of thumb to check

Are you saying that we currently do this already? I couldn't find this in the resample source, the only thing we're doing is min(lowpass, sfreq / 2.).

https://github.com/mne-tools/mne-python/blob/maint/0.19/mne/io/base.py#L1315-L1317

[larsoner]

We changed it to 1/3 here

https://github.com/mne-tools/mne-python/pull/6867/files

[cbrnr]

But that's just in Epochs.decim - the example uses raw.resample, where this is still sfreq / 2.

[larsoner]

But that's just in Epochs.decim - the example uses raw.resample, where this is still sfreq / 2.

In resample we don't need to sanity check because we know the filter we apply (brick wall), in decimate we look at info['lowpass'] and don't know characteristics so we use the 1/3 rule

[larsoner]

Travis says

tutorials/preprocessing/plot_30_filtering_resampling.py:171:80: E501 line too long (82 > 79 characters)

[cbrnr]

But @robertoostenveld referred to this section, which uses raw.resample and a cutoff frequency of 1/2 * fs, and he suggested that this should probably be 1/3 * fs.

[larsoner]

Yes this is what my text above is about.

• In our own .resample methods we create and apply the AA filter so the "rule of thumb" should not apply IMO (we know characteristics so don't need a rule of thumb).
• In .decimate methods we do not apply any AA filter ourselves so we use the rule of thumb to try to prevent users from making mistakes.

[cbrnr]

Alright, I didn't realize that you were responding to @robertoostenveld 😄.

[codecov]

Codecov Report

Merging #6882 into master will increase coverage by 0.01%.
The diff coverage is n/a.

@@            Coverage Diff             @@
##           master    #6882      +/-   ##
==========================================
+ Coverage   89.64%   89.65%   +0.01%     
==========================================
  Files         426      426              
  Lines       76048    76092      +44     
  Branches    12409    12415       +6     
==========================================
+ Hits        68170    68221      +51     
  Misses       5101     5101              
+ Partials     2777     2770       -7

[larsoner]

@drammock let me know when I should look

[drammock]

@larsoner @robertoostenveld revisions ready for review: https://15887-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/preprocessing/plot_30_filtering_resampling.html#resampling

[larsoner]

LGTM +1 for merge

[cbrnr]

Thanks @robertoostenveld and @drammock!