.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/zapline/plot_06_ctf_meg.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_zapline_plot_06_ctf_meg.py: ZapLine on MEG-like data with many co-equal noise components. ============================================================= This example demonstrates the auto-mode fix from Issue #34. With high-channel-count MEG recordings (e.g., CTF 275-magnetometer systems), line noise couples into many sensors through different lead-field paths, producing several co-equal strong DSS components. The classic ``iterative_outlier_removal`` selector underdetects in this regime because no single eigenvalue stands out from the others; the new ``auto_select_components_robust`` falls back to a log-space knee detector that catches the cluster boundary. The example uses synthetic data because the issue does not provide a public dataset, but the eigenvalue pattern reproduced here mirrors the one reported in https://github.com/mne-tools/mne-denoise/issues/34. Authors: Sina Esmaeili (sina.esmaeili@umontreal.ca) .. GENERATED FROM PYTHON SOURCE LINES 22-24 Imports ------- .. GENERATED FROM PYTHON SOURCE LINES 24-36 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy import signal from mne_denoise.dss.utils.selection import ( auto_select_components_robust, detect_eigenvalue_knee, iterative_outlier_removal, ) from mne_denoise.viz import plot_component_score_curve, plot_psd_comparison from mne_denoise.zapline import ZapLine .. GENERATED FROM PYTHON SOURCE LINES 37-42 Simulate a MEG-like recording ----------------------------- Three harmonics (50, 100, 150 Hz) couple into each of 64 channels with an independent phase and amplitude. This is a deliberately simple stand-in for the spatially heterogeneous line-noise pickup seen in real MEG. .. GENERATED FROM PYTHON SOURCE LINES 42-65 .. code-block:: Python rng = np.random.default_rng(34) sfreq = 400 n_channels = 64 n_times = sfreq * 60 # 60 s times = np.arange(n_times) / sfreq line = np.zeros((n_channels, n_times)) for freq in (50.0, 100.0, 150.0): phases = rng.uniform(0, 2 * np.pi, n_channels) amps = rng.normal(1.0, 0.3, n_channels) * 5.0 line += amps[:, None] * np.sin(2 * np.pi * freq * times[None, :] + phases[:, None]) # 1/f-ish brain background background = rng.normal(0, 1, (n_channels, n_times)) sos = signal.butter(4, [1, 80], btype="bandpass", fs=sfreq, output="sos") background = signal.sosfiltfilt(sos, background, axis=-1) # Add a small 10 Hz alpha rhythm that we want to preserve alpha = 0.3 * np.sin(2 * np.pi * 10.0 * times) background += alpha[None, :] data = background + line .. GENERATED FROM PYTHON SOURCE LINES 66-68 Fit ZapLine in auto mode ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 68-73 .. code-block:: Python zap = ZapLine(sfreq=sfreq, line_freq=50.0, n_remove="auto", n_harmonics=3) zap.fit(data) cleaned = zap.transform(data) .. GENERATED FROM PYTHON SOURCE LINES 74-81 Inspect what each selection strategy returns -------------------------------------------- Before the fix, ``ZapLine`` called ``iterative_outlier_removal`` directly. On the user's CTF data it returned 0 even though seven DSS eigenvalues clearly stood above the noise floor. The new ``auto_select_components_robust`` additionally runs ``detect_eigenvalue_knee`` and uses the larger of the two counts. .. GENERATED FROM PYTHON SOURCE LINES 81-99 .. code-block:: Python eigs = zap.eigenvalues_ n_outlier = iterative_outlier_removal(eigs, sigma=zap.threshold) n_knee = detect_eigenvalue_knee( eigs, rel_floor=zap.knee_rel_floor, min_ratio=zap.knee_min_ratio ) n_robust = auto_select_components_robust( eigs, sigma=zap.threshold, knee_rel_floor=zap.knee_rel_floor, knee_min_ratio=zap.knee_min_ratio, ) print(f"iterative_outlier_removal -> {n_outlier}") print(f"detect_eigenvalue_knee -> {n_knee}") print(f"auto_select_components_robust -> {n_robust}") print(f"ZapLine n_removed_ -> {zap.n_removed_}") .. rst-class:: sphx-glr-script-out .. code-block:: none iterative_outlier_removal -> 6 detect_eigenvalue_knee -> 6 auto_select_components_robust -> 6 ZapLine n_removed_ -> 6 .. GENERATED FROM PYTHON SOURCE LINES 100-104 Visualize the eigenvalue spectrum --------------------------------- The knee marks the boundary between the line-noise cluster and the background floor. .. GENERATED FROM PYTHON SOURCE LINES 104-111 .. code-block:: Python fig, ax = plt.subplots(figsize=(8, 4)) plot_component_score_curve(zap, ax=ax, show=False) ax.set_title("ZapLine DSS eigenvalues (auto-selected components highlighted)") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/zapline/images/sphx_glr_plot_06_ctf_meg_001.png :alt: ZapLine DSS eigenvalues (auto-selected components highlighted) :srcset: /auto_examples/zapline/images/sphx_glr_plot_06_ctf_meg_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 112-114 Confirm the cleaning worked --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 114-130 .. code-block:: Python fig, ax = plt.subplots(figsize=(8, 4)) plot_psd_comparison( data, cleaned, sfreq=sfreq, line_freq=50.0, fmin=1, fmax=180, ax=ax, show=False, ) ax.set_title("PSD before vs. after ZapLine (MEG-like, auto-mode)") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/zapline/images/sphx_glr_plot_06_ctf_meg_002.png :alt: PSD before vs. after ZapLine (MEG-like, auto-mode) :srcset: /auto_examples/zapline/images/sphx_glr_plot_06_ctf_meg_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 131-143 Workarounds for the rank-loss warning ------------------------------------- Real CTF recordings can trip the rank-loss warning in :func:`mne_denoise.dss.linear.compute_dss` when channels with very different physical units are fit together, or when the baseline covariance spectrum spans many decades. Three options: 1. Pass ``normalize_input=True`` to :class:`mne_denoise.dss.linear.DSS` when fitting DSS directly (ZapLine sets this to ``False`` by design). 2. Lower ``reg`` if the relative cutoff discards too many PCA dimensions. 3. Fit homogeneous channel types separately, e.g. use ``raw.copy().pick("mag")`` for CTF magnetometers instead of passing all channels at once. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.493 seconds) .. _sphx_glr_download_auto_examples_zapline_plot_06_ctf_meg.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_06_ctf_meg.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_06_ctf_meg.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_06_ctf_meg.zip `