Add notch filter (#702)

Author: hoechenberger

docs/source/settings/preprocessing/filter.md

@@ -32,3 +32,6 @@ of tips! 😇
         - h_freq
         - l_trans_bandwidth
         - h_trans_bandwidth
+        - notch_freq
+        - notch_trans_bandwidth
+        - notch_widths

docs/source/v1.1.md.inc

@@ -1,8 +1,8 @@
 ## v1.1.0-dev (unreleased)
 
-[//]: # (### :new: New features & enhancements)
+### :new: New features & enhancements
 
-[//]: # (- Whatever (#000 by @whoever) )
+- Add support for notch filtering (#702 by @hoechenberger and @larsoner)
 
 [//]: # (### :warning: Behavior changes)
 
@@ -12,4 +12,3 @@
 
 [//]: # (### :bug: Bug fixes)
 
-- Nothing yet

mne_bids_pipeline/_config.py

@@ -683,13 +683,32 @@
 l_freq: Optional[float] = None
 """
 The low-frequency cut-off in the highpass filtering step.
-Keep it None if no highpass filtering should be applied.
+Keep it `None` if no highpass filtering should be applied.
 """
 
 h_freq: Optional[float] = 40.0
 """
 The high-frequency cut-off in the lowpass filtering step.
-Keep it None if no lowpass filtering should be applied.
+Keep it `None` if no lowpass filtering should be applied.
+"""
+
+notch_freq: Optional[Union[float, Iterable[float]]] = None
+"""
+Notch filter frequency. More than one frequency can be supplied, e.g. to remove
+harmonics. Keep it `None` if no notch filter should be applied.
+
+Note: Note
+    The notch filter will be applied before high- and lowpass filtering.
+
+???+ example "Example"
+    Remove line noise at 50 Hz:
+    ```python
+    notch_freq = 50
+    ```
+    Remove line noise at 50 Hz and its (sub-)harmonics
+    ```python
+    notch_freq = [25, 50, 100, 150]
+    ```
 """
 
 l_trans_bandwidth: Union[float, Literal["auto"]] = "auto"
@@ -706,6 +725,16 @@
 parameters.
 """
 
+notch_trans_bandwidth: float = 1.0
+"""
+Specifies the transition bandwidth of the notch filter. The default is `1.`.
+"""
+
+notch_widths: Optional[Union[float, Iterable[float]]] = None
+"""
+Specifies the width of each stop band. `None` uses the MNE default.
+"""
+
 raw_resample_sfreq: Optional[float] = None
 """
 Specifies at which sampling frequency the data should be resampled.

mne_bids_pipeline/steps/preprocessing/_03_frequency_filter.py

@@ -16,7 +16,7 @@
 
 import numpy as np
 from types import SimpleNamespace
-from typing import Optional, Union, Literal
+from typing import Optional, Union, Literal, Iterable
 
 import mne
 
@@ -58,7 +58,36 @@ def get_input_fnames_frequency_filter(
     )
 
 
-def filter(
+def notch_filter(
+    raw: mne.io.BaseRaw,
+    subject: str,
+    session: Optional[str],
+    run: str,
+    freqs: Optional[Union[float, Iterable[float]]],
+    trans_bandwidth: Union[float, Literal["auto"]],
+    notch_widths: Optional[Union[float, Iterable[float]]],
+    data_type: Literal["experimental", "empty-room", "resting-state"],
+) -> None:
+    """Filter data channels (MEG and EEG)."""
+    if freqs is None:
+        msg = f"Not applying notch filter to {data_type} data."
+    else:
+        msg = f"Notch filtering {data_type} data at {freqs} Hz."
+
+    logger.info(**gen_log_kwargs(message=msg))
+
+    if freqs is None:
+        return
+
+    raw.notch_filter(
+        freqs=freqs,
+        trans_bandwidth=trans_bandwidth,
+        notch_widths=notch_widths,
+        n_jobs=1,
+    )
+
+
+def bandpass_filter(
     raw: mne.io.BaseRaw,
     subject: str,
     session: Optional[str],
@@ -89,10 +118,6 @@ def filter(
         h_freq=h_freq,
         l_trans_bandwidth=l_trans_bandwidth,
         h_trans_bandwidth=h_trans_bandwidth,
-        filter_length="auto",
-        phase="zero",
-        fir_window="hamming",
-        fir_design="firwin",
         n_jobs=1,
     )
 
@@ -152,7 +177,17 @@ def filter_data(
         split=None,
     )
     raw.load_data()
-    filter(
+    notch_filter(
+        raw=raw,
+        subject=subject,
+        session=session,
+        run=run,
+        freqs=cfg.notch_freq,
+        trans_bandwidth=cfg.notch_trans_bandwidth,
+        notch_widths=cfg.notch_widths,
+        data_type="experimental",
+    )
+    bandpass_filter(
         raw=raw,
         subject=subject,
         session=session,
@@ -179,10 +214,10 @@ def filter_data(
         split_size=cfg._raw_split_size,
     )
     _update_for_splits(out_files, in_key)
+    fmax = 1.5 * cfg.h_freq if cfg.h_freq is not None else np.inf
     if exec_params.interactive:
         # Plot raw data and power spectral density.
         raw.plot(n_channels=50, butterfly=True)
-        fmax = 1.5 * cfg.h_freq if cfg.h_freq is not None else np.inf
         raw.plot_psd(fmax=fmax)
 
     del raw
@@ -227,7 +262,17 @@ def filter_data(
         )
 
         raw_noise.load_data()
-        filter(
+        notch_filter(
+            raw=raw_noise,
+            subject=subject,
+            session=session,
+            run=task,
+            freqs=cfg.notch_freq,
+            trans_bandwidth=cfg.notch_trans_bandwidth,
+            notch_widths=cfg.notch_widths,
+            data_type=data_type,
+        )
+        bandpass_filter(
             raw=raw_noise,
             subject=subject,
             session=session,
@@ -257,7 +302,6 @@ def filter_data(
         if exec_params.interactive:
             # Plot raw data and power spectral density.
             raw_noise.plot(n_channels=50, butterfly=True)
-            fmax = 1.5 * cfg.h_freq if cfg.h_freq is not None else np.inf
             raw_noise.plot_psd(fmax=fmax)
 
     assert len(in_files) == 0, in_files.keys()
@@ -300,8 +344,11 @@ def get_config(
         deriv_root=config.deriv_root,
         l_freq=config.l_freq,
         h_freq=config.h_freq,
+        notch_freq=config.notch_freq,
         l_trans_bandwidth=config.l_trans_bandwidth,
         h_trans_bandwidth=config.h_trans_bandwidth,
+        notch_trans_bandwidth=config.notch_trans_bandwidth,
+        notch_widths=config.notch_widths,
         raw_resample_sfreq=config.raw_resample_sfreq,
         crop_runs=config.crop_runs,
         rename_events=config.rename_events,

mne_bids_pipeline/tests/configs/config_ERP_CORE.py

@@ -57,7 +57,8 @@
 eog_channels = ["HEOG", "VEOG"]
 
 l_freq = 0.1
-h_freq = 40
+h_freq = None
+notch_freq = 60
 
 decode = True
 decoding_time_generalization = True
@@ -148,6 +149,7 @@
         "theta": [4, 7],
         "alpha": [8, 12],
         "beta": [13, 20, 30],
+        "gamma": [50, 63],
     }
     decoding_csp_times = [-0.2, 0.0, 0.2, 0.4]
 elif task == "LRP":

mne_bids_pipeline/tests/test_documented.py

@@ -50,8 +50,9 @@ def test_options_documented():
                     if val not in allowed_duplicates:
                         assert val not in in_doc, "Duplicate documentation"
                     in_doc.add(val)
-    assert in_doc.difference(in_config) == set(), "Extra values in doc"
-    assert in_config.difference(in_doc) == set(), "Values missing from doc"
+    what = "docs/source/settings doc"
+    assert in_doc.difference(in_config) == set(), f"Extra values in {what}"
+    assert in_config.difference(in_doc) == set(), f"Values missing from {what}"
 
 
 def test_datasets_in_doc():