.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/bss_cca/plot_02_bss_cca_properties.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_02_bss_cca_properties.py: BSS-CCA algorithm properties ============================ A deterministic, synthetic demonstration of what BSS-CCA does and where it breaks. This is **not** a replication of any published experiment — every signal here is simulated, and the numbers depend on the simulation. For the paper-oriented experiments see ``scripts/replicate_declercq_paper.py``. Four properties are shown: 1. how the error depends on how many components you remove; 2. that keeping every component is exactly a no-op, including on rank-deficient data; 3. that the correlation spectrum of realistic broadband data is compressed, so a high fixed threshold is not a safe operating point; 4. that near-Nyquist energy inverts the correlation ordering unless the data is band-limited first. .. GENERATED FROM PYTHON SOURCE LINES 20-52 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from mne_denoise.bss_cca import compute_bss_cca from mne_denoise.viz import themed_figure sfreq = 250.0 n_channels, n_times = 21, int(250.0 * 10) rng = np.random.default_rng(11) times = np.arange(n_times) / sfreq brain = 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), ] ) muscle = rng.standard_normal((3, n_times)) clean = rng.standard_normal((n_channels, 3)) @ brain observed = clean + 0.8 * (rng.standard_normal((n_channels, 3)) @ muscle) clean_centered = clean - clean.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) ) .. GENERATED FROM PYTHON SOURCE LINES 53-57 1. Error as a function of how many components are removed --------------------------------------------------------- Three sources generated the muscle activity, and the error bottoms out at three. Removing more starts eating the neural signal. .. GENERATED FROM PYTHON SOURCE LINES 57-77 .. code-block:: Python # The sweep can only go as far as the number of components the data supports. _, base_info = compute_bss_cca(observed, n_remove=0) counts = np.arange(0, base_info["input_rank"] + 1) errors = [ relative_error(compute_bss_cca(observed, n_remove=int(k), preserve_mean=False)[0]) for k in counts ] for count, error in zip(counts, errors, strict=True): print(f"n_remove={count}: relative error {error:.4f}") fig, ax = themed_figure(figsize=(6.0, 3.6)) ax.plot(counts, errors, marker="o") ax.axvline(3, linestyle="--", linewidth=1, label="true number of muscle sources") ax.set_xlabel("Components removed") ax.set_ylabel("Relative error") ax.set_title("Error versus components removed (synthetic)") ax.legend() fig.tight_layout() .. image-sg:: /auto_examples/bss_cca/images/sphx_glr_plot_02_bss_cca_properties_001.png :alt: Error versus components removed (synthetic) :srcset: /auto_examples/bss_cca/images/sphx_glr_plot_02_bss_cca_properties_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none n_remove=0: relative error 1.2272 n_remove=1: relative error 1.0298 n_remove=2: relative error 0.8167 n_remove=3: relative error 0.0348 n_remove=4: relative error 0.5863 n_remove=5: relative error 0.7967 n_remove=6: relative error 1.0000 .. GENERATED FROM PYTHON SOURCE LINES 78-84 2. Keeping everything is exactly the identity --------------------------------------------- This has to hold for rank-deficient data too. Average referencing, channel interpolation, and duplicated channels all reduce the rank, and an implementation that back-projects by inverting the canonical filters silently destroys whole channels in exactly these cases. .. GENERATED FROM PYTHON SOURCE LINES 84-107 .. code-block:: Python # Start from a full-rank recording so each manipulation removes a known amount # of rank, and the reported ranks differ. full_rank = rng.standard_normal((12, n_times)) average_referenced = full_rank - full_rank.mean(axis=0, keepdims=True) duplicated = full_rank.copy() duplicated[-1] = duplicated[0] flat = full_rank.copy() flat[5] = 3.0 for label, data in ( ("full rank", full_rank), ("average reference", average_referenced), ("duplicated channel", duplicated), ("flat channel", flat), ): identity, info = compute_bss_cca(data, n_remove=0, preserve_mean=True) residual = np.abs(identity - data).max() print( f"{label:20s} rank={info['input_rank']:3d} of {data.shape[0]} " f"max deviation {residual:.2e}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none full rank rank= 12 of 12 max deviation 6.22e-15 average reference rank= 11 of 12 max deviation 4.55e-15 duplicated channel rank= 11 of 12 max deviation 4.33e-15 flat channel rank= 11 of 12 max deviation 5.55e-15 .. GENERATED FROM PYTHON SOURCE LINES 108-113 3. Realistic data has a compressed correlation spectrum ------------------------------------------------------- The narrow-band sources above separate cleanly. Broadband neural background does not: every correlation lands in a narrow range, and a high fixed threshold would reject the entire recording. .. GENERATED FROM PYTHON SOURCE LINES 113-149 .. code-block:: Python freqs = np.fft.rfftfreq(n_times, 1.0 / sfreq) def band_limited(n_sources, fmin, fmax): """Random sources restricted to a frequency band.""" out = np.zeros((n_sources, n_times)) for index in range(n_sources): spectrum = rng.standard_normal(freqs.size) + 1j * rng.standard_normal( freqs.size ) spectrum[(freqs < fmin) | (freqs > fmax)] = 0.0 out[index] = np.fft.irfft(spectrum, n=n_times) return out / np.linalg.norm(out) broadband_brain = band_limited(6, 1.0, 45.0) broadband_muscle = band_limited(3, 30.0, sfreq / 2) realistic = rng.standard_normal((n_channels, 6)) @ broadband_brain realistic = realistic + 1.2 * (rng.standard_normal((n_channels, 3)) @ broadband_muscle) _, narrow_info = compute_bss_cca(observed, n_remove=3) _, wide_info = compute_bss_cca(realistic, n_remove=3) print(f"narrow-band correlations: {np.round(narrow_info['correlations'][:6], 3)}") print(f"broadband correlations: {np.round(wide_info['correlations'][:6], 3)}") print(f"broadband maximum: {wide_info['correlations'].max():.3f}") fig, ax = themed_figure(figsize=(6.0, 3.6)) ax.plot(narrow_info["correlations"], marker="o", label="narrow-band sources") ax.plot(wide_info["correlations"], marker="s", label="broadband sources") ax.set_xlabel("Component") ax.set_ylabel("Canonical correlation") ax.set_title("Correlation spectrum depends on the signal regime (synthetic)") ax.legend() fig.tight_layout() .. image-sg:: /auto_examples/bss_cca/images/sphx_glr_plot_02_bss_cca_properties_002.png :alt: Correlation spectrum depends on the signal regime (synthetic) :srcset: /auto_examples/bss_cca/images/sphx_glr_plot_02_bss_cca_properties_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none narrow-band correlations: [0.999 0.969 0.864 0.078 0.023 0.017] broadband correlations: [0.817 0.809 0.805 0.796 0.783 0.77 ] broadband maximum: 0.817 .. GENERATED FROM PYTHON SOURCE LINES 150-156 4. Near-Nyquist energy inverts the ordering ------------------------------------------- Canonical correlations come from singular values and are never negative, so a component oscillating at the Nyquist frequency is anti-correlated at lag 1 yet receives a correlation near 1. Band-limit before decomposing, as both source papers do, and check ``autocorrelations`` for negative entries. .. GENERATED FROM PYTHON SOURCE LINES 156-189 .. code-block:: Python aliased_sources = np.zeros((6, n_times)) aliased_sources[0] = (-1.0) ** np.arange(n_times) aliased_sources[1:] = rng.standard_normal((5, n_times)) aliased = rng.standard_normal((6, 6)) @ aliased_sources _, aliased_info = compute_bss_cca(aliased, n_remove=1) worst = int(np.argmin(aliased_info["autocorrelations"])) print(f"correlations: {np.round(aliased_info['correlations'], 3)}") print(f"autocorrelations:{np.round(aliased_info['autocorrelations'], 3)}") print( f"component {worst} ranks at correlation " f"{aliased_info['correlations'][worst]:.3f} but is anti-correlated " f"({aliased_info['autocorrelations'][worst]:.3f})" ) fig, ax = themed_figure(figsize=(6.0, 3.6)) index = np.arange(aliased_info["correlations"].size) ax.bar(index - 0.2, aliased_info["correlations"], width=0.4, label="canonical rho") ax.bar( index + 0.2, aliased_info["autocorrelations"], width=0.4, label="signed autocorrelation", ) ax.axhline(0.0, linewidth=1) ax.set_xlabel("Component") ax.set_ylabel("Value") ax.set_title("Unfiltered data: rho hides the sign (synthetic)") ax.legend() fig.tight_layout() plt.show() .. image-sg:: /auto_examples/bss_cca/images/sphx_glr_plot_02_bss_cca_properties_003.png :alt: Unfiltered data: rho hides the sign (synthetic) :srcset: /auto_examples/bss_cca/images/sphx_glr_plot_02_bss_cca_properties_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none correlations: [1. 0.065 0.053 0.044 0.017 0.006] autocorrelations:[-1. -0.031 -0.004 0.039 -0.002 0.001] component 0 ranks at correlation 1.000 but is anti-correlated (-1.000) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.607 seconds) .. _sphx_glr_download_auto_examples_bss_cca_plot_02_bss_cca_properties.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_02_bss_cca_properties.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_02_bss_cca_properties.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_02_bss_cca_properties.zip `