feat: Adaptive DSS with segmentation, auto-selection, and smoothing

docs/changes/devel/feature.rst

@@ -0,0 +1 @@
+Add adaptive DSS with automatic segmentation (:class:`~mne_denoise.dss.utils.CovarianceSegmenter`, :class:`~mne_denoise.dss.utils.FixedWindowSegmenter`), component auto-selection (:func:`~mne_denoise.dss.utils.eigenvalue_ratio_selection`, :func:`~mne_denoise.dss.utils.max_gap_selection`), smoothing decomposition, L2-normalised patterns, and proportional Q-mode for the periodic denoiser, by `Scott Huberty`_.

examples/dss/plot_01_dss_fundamentals.py

@@ -202,7 +202,7 @@
 print("Reconstructing data from first component...")
 sources = dss_evoked.transform(epochs)
 # To reconstruct using only specific components, we zero out the others
-sources[:, 1:, :] = 0
+sources[1:, :, :] = 0
 epochs_denoised = dss_evoked.inverse_transform(sources)
 epochs_denoised = mne.EpochsArray(epochs_denoised, info)
 
@@ -270,7 +270,7 @@
 # We concatenate epochs for continuous reconstruction if desired, or keep as epochs
 # Here we keep as epochs to use plot_psd_comparison
 sources = dss_osc.transform(epochs)
-sources[:, 1:, :] = 0
+sources[1:, :, :] = 0
 epochs_osc = dss_osc.inverse_transform(sources)
 epochs_osc = mne.EpochsArray(epochs_osc, info)
 
@@ -424,7 +424,7 @@
 # Denoising Comparison
 print("Reconstructing M100 component...")
 sources = dss_m100.transform(epochs_real)
-sources[:, 1:, :] = 0
+sources[1:, :, :] = 0
 epochs_m100 = dss_m100.inverse_transform(sources)
 epochs_m100 = mne.EpochsArray(epochs_m100, epochs_real.info)
 

examples/dss/plot_04_spectral_dss.py

@@ -220,7 +220,7 @@
 
 # We use method='iir' to replicate a traditional Notch filter approach
 notch_bias_60 = LineNoiseBias(freq=60, sfreq=sfreq, method="iir", bandwidth=2)
-dss_notch = DSS(n_components=None, bias=notch_bias_60)
+dss_notch = DSS(n_components=3, bias=notch_bias_60)
 dss_notch.fit(raw_noisy)
 
 print(f"\nDSS Eigenvalues: {dss_notch.eigenvalues_[:3]}")

examples/zapline/plot_04_adaptive_mode.py

@@ -17,12 +17,12 @@
 from matplotlib.gridspec import GridSpec
 from scipy import signal
 
+from mne_denoise.dss.utils.segmentation import CovarianceSegmenter
 from mne_denoise.viz import plot_psd_comparison
 from mne_denoise.zapline.adaptive import (
     check_artifact_presence,
     find_fine_peak,
     find_noise_freqs,
-    segment_data,
 )
 
 # Suppress warnings for cleaner output
@@ -204,9 +204,12 @@ def generate_nonstationary_data(
 
 print("\n--- Step 2: Adaptive Segmentation ---")
 target_freq = detected_freqs[0] if detected_freqs else 50.0
-segments = segment_data(
-    data, sfreq, target_freq, min_chunk_len=PAPER_PARAMS["minChunkLength"]
+segmenter = CovarianceSegmenter(
+    sfreq=sfreq,
+    min_chunk_len=PAPER_PARAMS["minChunkLength"],
+    bandpass=(target_freq - 3, target_freq + 3),
 )
+segments = segmenter.segment(data)
 print(f"Number of segments: {len(segments)}")
 for i, (start, end) in enumerate(segments):
     print(

mne_denoise/dss/__init__.py

@@ -40,6 +40,13 @@
 
 # Utils (exposed for convenience if needed)
 from .utils import convergence, whitening
+from .utils.segmentation import CovarianceSegmenter, FixedWindowSegmenter
+from .utils.selection import (
+    auto_select_components,
+    eigenvalue_ratio_selection,
+    iterative_outlier_removal,
+    max_gap_selection,
+)
 
 # Variants (Modules)
 from .variants import narrowband, ssvep, tsr

mne_denoise/dss/denoisers/periodic.py

@@ -140,6 +140,12 @@ class CombFilterBias(LinearDenoiser):
         Default 3.
     q_factor : float
         Quality factor for each peak. Default 30.
+    q_mode : ``"fixed"`` | ``"proportional"``
+        How Q scales across harmonics. ``"fixed"`` uses the same
+        ``q_factor`` for every harmonic (bandwidth narrows as frequency
+        increases). ``"proportional"`` scales Q as ``q_factor * h`` for
+        the *h*-th harmonic, maintaining approximately constant absolute
+        bandwidth across all harmonics. Default ``"fixed"``.
     weights : array-like, optional
         Weights for each harmonic. If None, uses 1/harmonic_number
         weighting (decreasing importance of higher harmonics).
@@ -154,6 +160,10 @@ class CombFilterBias(LinearDenoiser):
     >>> bias = CombFilterBias(
     ...     fundamental_freq=12, sfreq=500, n_harmonics=4, weights=[1.0, 1.0, 1.0, 1.0]
     ... )
+    >>> # Adaptive Q for constant bandwidth across harmonics
+    >>> bias = CombFilterBias(
+    ...     fundamental_freq=50, sfreq=1000, n_harmonics=3, q_mode="proportional"
+    ... )
 
     See Also
     --------
@@ -177,13 +187,22 @@ def __init__(
         *,
         n_harmonics: int = 3,
         q_factor: float = 30.0,
+        q_mode: str = "fixed",
         weights: np.ndarray | None = None,
     ) -> None:
         self.fundamental_freq = fundamental_freq
         self.sfreq = sfreq
         self.n_harmonics = n_harmonics
         self.q_factor = q_factor
 
+        # Validate q_mode
+        allowed_q_modes = ("fixed", "proportional")
+        if q_mode not in allowed_q_modes:
+            raise ValueError(
+                f"q_mode must be one of {allowed_q_modes}, got {q_mode!r}"
+            )
+        self.q_mode = q_mode
+
         # Set up weights
         if weights is None:
             self.weights = np.array([1.0 / h for h in range(1, n_harmonics + 1)])
@@ -212,7 +231,11 @@ def _create_harmonic_filters(self) -> None:
             w0 = freq / nyq
             weight = self.weights[h - 1]
 
-            b, a = signal.iirpeak(w0, self.q_factor)
+            # Proportional Q scales linearly with harmonic number,
+            # maintaining constant absolute bandwidth across harmonics
+            q = self.q_factor * h if self.q_mode == "proportional" else self.q_factor
+
+            b, a = signal.iirpeak(w0, q)
             sos = signal.tf2sos(b, a)
             self._peak_filters.append((sos, weight))
 

mne_denoise/dss/denoisers/temporal.py

@@ -12,6 +12,10 @@
        Journal of Neuroscience Methods, 189(1), 113-120.
 .. [2] de Cheveigné, A. & Simon, J.Z. (2008). Denoising based on spatial filtering.
        Journal of Neuroscience Methods, 171(2), 331-339.
+.. [3] de Cheveigné, A. (2020). ZapLine: A simple and effective method to remove
+       power line artifacts. NeuroImage, 207, 116356. (Period-matched
+       smooth/residual decomposition: spatially clean only the residual branch
+       and add the smooth branch back.)
 """
 
 from __future__ import annotations
@@ -151,7 +155,12 @@ def _prediction_bias(self, data: np.ndarray) -> np.ndarray:
 class SmoothingBias(LinearDenoiser):
     """Unified temporal smoothing bias (Moving Average).
 
-    Uses a boxcar moving average filter to smooth the data."
+    Uses a boxcar moving average filter to smooth the data. When used to split
+    the signal into a smooth branch and a residual (``data - smooth``), fitting
+    DSS on the residual and adding the smooth branch back follows ZapLine's
+    period-matched decomposition (de Cheveigné, 2020 [3]_): with
+    ``window = round(sfreq / f_line)`` the smoother has zeros at ``f_line`` and
+    its harmonics, so the residual concentrates the narrowband artifact.
 
     Parameters
     ----------