Add tutorial

tutorials/preprocessing/20_rejecting_bad_data.py

@@ -23,6 +23,8 @@
 
 import mne
 
+import numpy as np
+
 sample_data_folder = mne.datasets.sample.data_path()
 sample_data_raw_file = os.path.join(
     sample_data_folder, "MEG", "sample", "sample_audvis_filt-0-40_raw.fif"
@@ -203,8 +205,8 @@
 # %%
 # .. _`tut-reject-epochs-section`:
 #
-# Rejecting Epochs based on channel amplitude
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+# Rejecting Epochs based on peak-to-peak channel amplitude
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 #
 # Besides "bad" annotations, the :class:`mne.Epochs` class constructor has
 # another means of rejecting epochs, based on signal amplitude thresholds for
@@ -326,6 +328,99 @@
 epochs.drop_bad(reject=stronger_reject_criteria)
 print(epochs.drop_log)
 
+# %%
+# .. _`tut-reject-epochs-func-section`:
+#
+# Rejecting Epochs using callables (functions)
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+# Sometimes it is useful to reject epochs based criteria other than
+# peak-to-peak amplitudes. For example, we might want to reject epochs
+# based on the maximum or minimum amplitude of a channel.
+# In this case, the :class:`mne.Epochs` class constructor also accepts
+# callables (functions) in the ``reject`` and ``flat`` parameters. This
+# allows us to define functions to reject epochs based on our desired criteria.
+#
+# Let's begin by generating Epoch data with large artifacts in one eeg channel
+# in order to demonstrate the versatility of this approach.
+
+raw.crop(0, 5)
+raw.del_proj()
+chans = raw.info['ch_names'][-5:-1]
+raw.pick(chans)
+data = raw.get_data()
+
+new_data = data
+new_data[0, 180:200] *= 1e3
+new_data[0, 460:580] += 1e-3
+edit_raw = mne.io.RawArray(new_data, raw.info)
+
+# Create fixed length epochs of 1 second
+events = mne.make_fixed_length_events(edit_raw, id=1, duration=1.0, start=0)
+epochs = mne.Epochs(edit_raw, events, tmin=0, tmax=1, baseline=None)
+epochs.plot(scalings=dict(eeg=50e-5))
+
+# %%
+# As you can see, we have two large artifacts in the first channel. One large
+# spike in amplitude and one large increase in amplitude.
+
+# Let's try to reject the epoch containing the spike in amplitude based on the
+# maximum amplitude of the first channel.
+
+epochs = mne.Epochs(
+    edit_raw, 
+    events,
+    tmin=0,
+    tmax=1,
+    baseline=None,
+    reject=dict(eeg=lambda x: True if (np.max(x, axis=1) > 1e-2).any() else False),
+    preload=True
+)
+epochs.plot(scalings=dict(eeg=50e-5))
+
+# %%
+# Here, the epoch containing the spike in amplitude was rejected for having a
+# maximum amplitude greater than 1e-2 Volts. Notice the use of the ``any()``
+# function to check if any of the channels exceeded the threshold. We could
+# have also used the ``all()`` function to check if all channels exceeded the
+# threshold.
+
+# Next, let's try to reject the epoch containing the increase in amplitude
+# using the median.
+
+epochs = mne.Epochs(
+    edit_raw,
+    events,
+    tmin=0,
+    tmax=1,
+    baseline=None,
+    reject=dict(eeg=lambda x: True if (np.median(x, axis=1) > 1e-4).any() else False),
+    preload=True
+)
+epochs.plot(scalings=dict(eeg=50e-5))
+
+# %%
+# Finally, let's try to reject both epochs using a combination of the maximum
+# and median. We'll define a custom function and use boolean operators to
+# combine the two criteria.
+
+
+def reject_criteria(x):
+    max_condition = np.max(x, axis=1) > 1e-2
+    median_condition = np.median(x, axis=1) > 1e-4
+    return True if max_condition.any() or median_condition.any() else False
+
+
+epochs = mne.Epochs(
+    edit_raw,
+    events,
+    tmin=0,
+    tmax=1,
+    baseline=None,
+    reject=dict(eeg=reject_criteria),
+    preload=True
+)
+epochs.plot(events=True)
+
 # %%
 # Note that a complementary Python module, the `autoreject package`_, uses
 # machine learning to find optimal rejection criteria, and is designed to