.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/bss_cca/plot_01_bss_cca_basics.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_bss_cca_plot_01_bss_cca_basics.py: Basic reference-free BSS-CCA ============================ This example mixes narrow-band neural sources with broadband muscle sources, attenuates the muscle activity with BSS-CCA, and then runs the same workflow on an MNE Raw object. BSS-CCA needs no reference channel: it separates sources by how correlated they are with a delayed copy of the recording, and muscle activity — being close to temporally white — lands in the lowest-correlation components. .. GENERATED FROM PYTHON SOURCE LINES 13-39 .. code-block:: Python import mne import numpy as np from mne_denoise.bss_cca import BSSCCA, compute_bss_cca from mne_denoise.viz import plot_psd_comparison, plot_signal_overlay sfreq = 250.0 n_channels, n_times = 21, int(250.0 * 10) rng = np.random.default_rng(2006) times = np.arange(n_times) / sfreq # Narrow-band "brain" sources: delta, alpha, and low beta. brain_sources = np.vstack( [ np.sin(2 * np.pi * 2.0 * times), np.sin(2 * np.pi * 10.0 * times + 0.4), np.sin(2 * np.pi * 21.0 * times + 1.1), ] ) # Broadband "muscle" sources. muscle_sources = rng.standard_normal((3, n_times)) clean = rng.standard_normal((n_channels, 3)) @ brain_sources observed = clean + 0.8 * (rng.standard_normal((n_channels, 3)) @ muscle_sources) .. GENERATED FROM PYTHON SOURCE LINES 40-43 Removing the three lowest-correlation components recovers the neural signal. ``n_remove`` is the operating knob used in the original paper; sweeping it is how you choose an operating point. .. GENERATED FROM PYTHON SOURCE LINES 43-62 .. code-block:: Python cleaned, info = compute_bss_cca(observed, n_remove=3, sfreq=sfreq, preserve_mean=False) clean_centered = clean - clean.mean(axis=1, keepdims=True) observed_centered = observed - observed.mean(axis=1, keepdims=True) def relative_error(estimate): """Relative RMS error against the known clean signal.""" return float( np.linalg.norm(estimate - clean_centered) / np.linalg.norm(clean_centered) ) print(f"Canonical correlations: {np.round(info['correlations'], 3)}") print(f"Components removed: {info['n_removed']} of {info['input_rank']}") print(f"Relative error before: {relative_error(observed_centered):.3f}") print(f"Relative error after: {relative_error(cleaned):.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Canonical correlations: [0.999 0.969 0.864 0.06 0.027 0.017] Components removed: 3 of 6 Relative error before: 1.058 Relative error after: 0.047 .. GENERATED FROM PYTHON SOURCE LINES 63-68 The signed autocorrelation is the honest diagnostic. Canonical correlations cannot be negative, so a component carrying energy near the Nyquist frequency would rank high despite being anti-correlated at the lag. What matters is whether such a component is *retained*: a slightly negative value among the components you already discarded is just noise. .. GENERATED FROM PYTHON SOURCE LINES 68-74 .. code-block:: Python kept = info["kept_mask"] retained_aliased = kept & (info["autocorrelations"] < -0.1) print(f"Signed lag-1 autocorrelation: {np.round(info['autocorrelations'], 3)}") print(f"Retained but anti-correlated: {np.flatnonzero(retained_aliased).tolist()}") .. rst-class:: sphx-glr-script-out .. code-block:: none Signed lag-1 autocorrelation: [ 0.999 0.969 0.864 -0.054 0.022 -0.014] Retained but anti-correlated: [] .. GENERATED FROM PYTHON SOURCE LINES 75-76 One channel, before and after, against the known clean waveform. .. GENERATED FROM PYTHON SOURCE LINES 76-94 .. code-block:: Python plot_signal_overlay( observed, cleaned, times, pick=0, start=0.0, stop=2.0, scale_after=False, before_label="Observed", after_label="BSS-CCA", reference=clean_centered[0], reference_label="Clean signal", x_label="Time (s)", title="Reference-free BSS-CCA - channel 0", show=False, ) .. image-sg:: /auto_examples/bss_cca/images/sphx_glr_plot_01_bss_cca_basics_001.png :alt: Reference-free BSS-CCA - channel 0 :srcset: /auto_examples/bss_cca/images/sphx_glr_plot_01_bss_cca_basics_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none
.. GENERATED FROM PYTHON SOURCE LINES 95-96 Broadband power falls while the narrow-band neural peaks survive. .. GENERATED FROM PYTHON SOURCE LINES 96-106 .. code-block:: Python plot_psd_comparison( observed, cleaned, sfreq=sfreq, fmin=1.0, fmax=120.0, show=False, ) .. image-sg:: /auto_examples/bss_cca/images/sphx_glr_plot_01_bss_cca_basics_002.png :alt: PSD Comparison :srcset: /auto_examples/bss_cca/images/sphx_glr_plot_01_bss_cca_basics_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none
.. GENERATED FROM PYTHON SOURCE LINES 107-109 MNE containers use the same estimator and return a copy of the same object, with timing, annotations, bad channels, and unselected channels preserved. .. GENERATED FROM PYTHON SOURCE LINES 109-117 .. code-block:: Python info_mne = mne.create_info( [f"EEG {index:03d}" for index in range(n_channels)], sfreq, "eeg" ) raw = mne.io.RawArray(observed, info_mne, verbose=False) raw_clean = BSSCCA(n_remove=3).fit_transform(raw) print(type(raw_clean).__name__, raw_clean.get_data().shape) .. rst-class:: sphx-glr-script-out .. code-block:: none RawArray (21, 2500) .. GENERATED FROM PYTHON SOURCE LINES 118-121 Muscle artifacts are non-stationary, so a single operator across a long recording is often the wrong model. ``segment_len`` reproduces the contiguous 10-second block scheme used clinically. .. GENERATED FROM PYTHON SOURCE LINES 121-124 .. code-block:: Python raw_blocked = BSSCCA(segment_len=10.0, n_remove=3).fit_transform(raw) print(f"Block-wise output: {raw_blocked.get_data().shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Block-wise output: (21, 2500) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.498 seconds) .. _sphx_glr_download_auto_examples_bss_cca_plot_01_bss_cca_basics.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_bss_cca_basics.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_bss_cca_basics.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_bss_cca_basics.zip `