.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/example_dss.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_example_dss.py: DSS example =========== Find the linear combinations of multichannel data that maximize repeatability over trials. This example uses synthetic data with a known repeated source, so we can check whether the first DSS component better matches the ground truth than the raw trial average. Uses meegkit.dss0(). References ---------- .. [1] de Cheveigne, A., & Parra, L. C. (2014). Joint decorrelation, a versatile tool for multichannel data analysis. NeuroImage, 98, 487-505. .. GENERATED FROM PYTHON SOURCE LINES 19-27 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from meegkit import dss from meegkit.utils import fold, rms, tscov, unfold rng = np.random.default_rng(5) .. GENERATED FROM PYTHON SOURCE LINES 28-30 Create simulated data ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 30-56 .. code-block:: Python # Data are time * channel * trials. n_samples = 100 * 3 n_chans = 30 n_trials = 100 noise_dim = 20 # dimensionality of noise # Source signal source = np.hstack(( np.zeros((n_samples // 3,)), np.sin(2 * np.pi * np.arange(n_samples // 3) / (n_samples / 3)).T, np.zeros((n_samples // 3,))))[np.newaxis].T s = source * rng.standard_normal((1, n_chans)) # 300 * 30 s = s[:, :, np.newaxis] s = np.tile(s, (1, 1, 100)) # Noise noise = np.dot( unfold(rng.standard_normal((n_samples, noise_dim, n_trials))), rng.standard_normal((noise_dim, n_chans))) noise = fold(noise, n_samples) # Mix signal and noise SNR = 0.1 data = noise / rms(noise.flatten()) + SNR * s / rms(s.flatten()) .. GENERATED FROM PYTHON SOURCE LINES 57-59 Apply DSS to clean them ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 59-90 .. code-block:: Python # Compute original and biased covariance matrices c0, _ = tscov(data) # In this case the biased covariance is simply the covariance of the mean over # trials c1, _ = tscov(np.mean(data, 2)) # Apply DSS [todss, _, pwr0, pwr1] = dss.dss0(c0, c1) z = fold(np.dot(unfold(data), todss), epoch_size=n_samples) # Find best components best_comp = np.mean(z[:, 0, :], -1) trial_avg_by_channel = np.mean(data, 2) raw_average = trial_avg_by_channel.mean(axis=1) # DSS components are sign-indeterminate. Flip to align with the known source # so the comparison is easier to read. comp_sign = np.sign(np.corrcoef(source[:, 0], best_comp)[0, 1]) if comp_sign == 0: comp_sign = 1.0 best_comp = comp_sign * best_comp q10, q90 = np.percentile(trial_avg_by_channel, [10, 90], axis=1) shown_chans = np.linspace(0, n_chans - 1, min(n_chans, 8), dtype=int) # Compare the recovered component to the known source waveform. source_corr_raw = np.corrcoef(source[:, 0], raw_average)[0, 1] source_corr_dss = np.corrcoef(source[:, 0], best_comp)[0, 1] .. GENERATED FROM PYTHON SOURCE LINES 91-97 Plot results ----------------------------------------------------------------------------- What to look for: - The raw trial average should still be visibly noisy. - The first DSS component should follow the known source more closely. - The reported correlation with the source should improve after DSS. .. GENERATED FROM PYTHON SOURCE LINES 97-128 .. code-block:: Python f, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 7), sharex=True, constrained_layout=True) time = np.arange(n_samples) ax1.plot(source, label="source") ax1.grid(alpha=.2) ax1.set_ylabel("Amplitude") ax1.set_title("Ground-truth source") channel_lines = ax2.plot(time, trial_avg_by_channel[:, shown_chans], color="0.7", lw=.8, alpha=.5) channel_lines[0].set_label(f"{len(shown_chans)} channels") for line in channel_lines[1:]: line.set_label("_nolegend_") ax2.plot(time, raw_average, color="C0", lw=2.2, label="channel mean") ax2.fill_between(time, q10, q90, color="C0", alpha=.18, linewidth=0, label="10-90% range") ax2.axhline(0, color="0.5", lw=.8, ls=":", zorder=0) ax2.grid(alpha=.2) ax2.set_ylabel("Amplitude") ax2.set_title(f"Trial-averaged observed data (corr = {source_corr_raw:.2f})") ax3.plot(time, best_comp, label="recovered") ax3.axhline(0, color="0.5", lw=.8, ls=":", zorder=0) ax3.grid(alpha=.2) ax3.set_ylabel("Amplitude") ax3.set_xlabel("Samples") ax3.set_title(f"Recovered DSS component (corr = {source_corr_dss:.2f})") ax1.legend(loc="upper right", fontsize="small") ax2.legend(loc="upper right", fontsize="small") ax3.legend(loc="upper right", fontsize="small") print(f"Correlation with source, raw trial average: {source_corr_raw:.3f}") print(f"Correlation with source, first DSS component: {source_corr_dss:.3f}") plt.show() .. image-sg:: /auto_examples/images/sphx_glr_example_dss_001.png :alt: Ground-truth source, Trial-averaged observed data (corr = -0.71), Recovered DSS component (corr = 1.00) :srcset: /auto_examples/images/sphx_glr_example_dss_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Correlation with source, raw trial average: -0.712 Correlation with source, first DSS component: 1.000 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.402 seconds) .. _sphx_glr_download_auto_examples_example_dss.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_dss.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_dss.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_dss.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_