Added implementation of zapline for power noise removal (#1032)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric Larson <[email protected]>

.circleci/config.yml

@@ -779,6 +779,7 @@ jobs:
 
   test_eeg_matchingpennies:
     <<: *imageconfig
+    resource_class: large  # memory for zapline
     steps:
       - attach_workspace:
           at: ~/

mne_bids_pipeline/_config.py

@@ -859,6 +859,17 @@
 Specifies the width of each stop band. `None` uses the MNE default.
 """
 
+zapline_fline: float | None = None
+"""
+Specifies frequency to remove using Zapline filtering. If None, zapline will not
+be used.
+"""
+
+zapline_iter: bool = False
+"""
+Specifies if the iterative version of the Zapline algorithm should be run.
+"""
+
 # ### Resampling
 #
 # If you have acquired data with a very high sampling frequency (e.g. 2 kHz)

mne_bids_pipeline/steps/preprocessing/_04_frequency_filter.py

@@ -20,6 +20,7 @@
 
 import mne
 import numpy as np
+from meegkit import dss
 from mne.io.pick import _picks_to_idx
 from mne.preprocessing import EOGRegression
 
@@ -64,6 +65,29 @@ def get_input_fnames_frequency_filter(
     )
 
 
+def zapline(
+    raw: mne.io.BaseRaw,
+    subject: str,
+    session: str | None,
+    run: str,
+    task: str | None,
+    fline: float | None,
+    iter_: bool,
+) -> None:
+    """Use Zapline to remove line frequencies."""
+    if fline is None:
+        return
+
+    msg = f"Zapline filtering data at with {fline=} Hz."
+    logger.info(**gen_log_kwargs(message=msg))
+    sfreq = raw.info["sfreq"]
+    picks = mne.pick_types(raw.info, meg=True, eeg=True)
+    data = raw.get_data(picks).T  # transpose to (n_samples, n_channels)
+    func = dss.dss_line_iter if iter_ else dss.dss_line
+    out, _ = func(data, fline, sfreq)
+    raw._data[picks] = out.T
+
+
 def notch_filter(
     raw: mne.io.BaseRaw,
     subject: str,
@@ -218,6 +242,15 @@ def filter_data(
         picks = np.unique(np.r_[picks_regress, picks_artifact, picks_data])
 
     raw.load_data()
+    zapline(
+        raw=raw,
+        subject=subject,
+        session=session,
+        run=run,
+        task=task,
+        fline=cfg.zapline_fline,
+        iter_=cfg.zapline_iter,
+    )
     notch_filter(
         raw=raw,
         subject=subject,
@@ -301,6 +334,8 @@ def get_config(
         l_freq=config.l_freq,
         h_freq=config.h_freq,
         notch_freq=config.notch_freq,
+        zapline_fline=config.zapline_fline,
+        zapline_iter=config.zapline_iter,
         l_trans_bandwidth=config.l_trans_bandwidth,
         h_trans_bandwidth=config.h_trans_bandwidth,
         notch_trans_bandwidth=config.notch_trans_bandwidth,

mne_bids_pipeline/tests/configs/config_eeg_matchingpennies.py

@@ -13,3 +13,8 @@
 decode = True
 
 interpolate_bads_grand_average = False
+
+l_freq = None
+h_freq = 100
+zapline_fline = 50
+zapline_iter = False

pyproject.toml

@@ -49,6 +49,7 @@ dependencies = [
   "mne[hdf5] >=1.7",
   "mne-bids[full]",
   "filelock",
+  "meegkit"
 ]
 dynamic = ["version"]