.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/csc/plot_cross_frequency_coupling.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_csc_plot_cross_frequency_coupling.py: ================================================================== Extracting cross-frequency coupling waveforms from rodent LFP data ================================================================== This example illustrates how to learn univariate atoms on a univariate time-serie. The data is a single LFP channel recorded on a rodent's striatum [1]_. Interestingly in this time-serie, the high frequency oscillations around 80 Hz are modulated in amplitude by the low-frequency oscillation around 3 Hz, a phenomenon known as cross-frequency coupling (CFC). The convolutional sparse coding (CSC) model is able to learn the prototypical waveforms of the signal, on which we can clearly see the CFC. .. [1] G. Dallérac, M. Graupner, J. Knippenberg, R. C. R. Martinez, T. F. Tavares, L. Tallot, N. El Massioui, A. Verschueren, S. Höhn, J.B. Bertolus, et al. Updating temporal expectancy of an aversive event engages striatal plasticity under amygdala control. Nature Communications, 8:13920, 2017 .. GENERATED FROM PYTHON SOURCE LINES 21-29 .. code-block:: Python # Authors: Tom Dupre La Tour # Mainak Jas # Umut Simsekli # Alexandre Gramfort # # License: BSD (3-clause) .. GENERATED FROM PYTHON SOURCE LINES 30-31 Let us first load the data sample. .. GENERATED FROM PYTHON SOURCE LINES 31-43 .. code-block:: Python import mne import numpy as np import matplotlib.pyplot as plt # sample frequency sfreq = 350. # We load the signal. It is an LFP channel recorded on a rodent's striatum. data = np.load('../rodent_striatum.npy') print(data.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (1, 630000) .. GENERATED FROM PYTHON SOURCE LINES 44-46 As the data contains severe artifacts between t=0 and t=100, we use a section not affected by artifacts. .. GENERATED FROM PYTHON SOURCE LINES 46-56 .. code-block:: Python data = data[:, 35000:] # We also remove the slow drift, which accounts for a lot of variance. data = mne.filter.filter_data(data, sfreq, 1, None) # To make the most of parallel computing, we split the data into trials. data = data.reshape(50, -1) data /= data.std() .. rst-class:: sphx-glr-script-out .. code-block:: none Setting up high-pass filter at 1 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal highpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Lower passband edge: 1.00 - Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz) - Filter length: 1155 samples (3.300 s) .. GENERATED FROM PYTHON SOURCE LINES 57-60 This sample contains CFC between 3 Hz and 80 Hz. This phenomenon can be described with a comodulogram, computed for instance with the `pactools `_ Python library. .. GENERATED FROM PYTHON SOURCE LINES 60-69 .. code-block:: Python from pactools import Comodulogram comod = Comodulogram(fs=sfreq, low_fq_range=np.arange(0.2, 10.2, 0.2), low_fq_width=2., method='duprelatour') comod.fit(data) comod.plot() plt.show() .. image-sg:: /auto_examples/csc/images/sphx_glr_plot_cross_frequency_coupling_001.png :alt: plot cross frequency coupling :srcset: /auto_examples/csc/images/sphx_glr_plot_cross_frequency_coupling_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [ ] 0% | 0.00 sec | comodulogram: DAR(10, 1) [ ] 2% | 0.57 sec | comodulogram: DAR(10, 1) [. ] 4% | 0.90 sec | comodulogram: DAR(10, 1) [.. ] 6% | 1.23 sec | comodulogram: DAR(10, 1) [... ] 8% | 1.56 sec | comodulogram: DAR(10, 1) [.... ] 10% | 1.89 sec | comodulogram: DAR(10, 1) [.... ] 12% | 2.24 sec | comodulogram: DAR(10, 1) [..... ] 14% | 2.57 sec | comodulogram: DAR(10, 1) [...... ] 16% | 2.91 sec | comodulogram: DAR(10, 1) [....... ] 18% | 3.25 sec | comodulogram: DAR(10, 1) [........ ] 20% | 3.56 sec | comodulogram: DAR(10, 1) [........ ] 22% | 3.90 sec | comodulogram: DAR(10, 1) [......... ] 24% | 4.25 sec | comodulogram: DAR(10, 1) [.......... ] 26% | 4.58 sec | comodulogram: DAR(10, 1) [........... ] 28% | 4.90 sec | comodulogram: DAR(10, 1) [............ ] 30% | 5.24 sec | comodulogram: DAR(10, 1) [............ ] 32% | 5.58 sec | comodulogram: DAR(10, 1) [............. ] 34% | 5.92 sec | comodulogram: DAR(10, 1) [.............. ] 36% | 6.24 sec | comodulogram: DAR(10, 1) [............... ] 38% | 6.55 sec | comodulogram: DAR(10, 1) [................ ] 40% | 6.89 sec | comodulogram: DAR(10, 1) [................ ] 42% | 7.23 sec | comodulogram: DAR(10, 1) [................. ] 44% | 7.54 sec | comodulogram: DAR(10, 1) [.................. ] 46% | 7.88 sec | comodulogram: DAR(10, 1) [................... ] 48% | 8.22 sec | comodulogram: DAR(10, 1) [.................... ] 50% | 8.52 sec | comodulogram: DAR(10, 1) [.................... ] 52% | 8.83 sec | comodulogram: DAR(10, 1) [..................... ] 54% | 9.14 sec | comodulogram: DAR(10, 1) [...................... ] 56% | 9.47 sec | comodulogram: DAR(10, 1) [....................... ] 58% | 9.78 sec | comodulogram: DAR(10, 1) [........................ ] 60% | 10.10 sec | comodulogram: DAR(10, 1) [........................ ] 62% | 10.40 sec | comodulogram: DAR(10, 1) [......................... ] 64% | 10.71 sec | comodulogram: DAR(10, 1) [.......................... ] 66% | 11.02 sec | comodulogram: DAR(10, 1) [........................... ] 68% | 11.32 sec | comodulogram: DAR(10, 1) [............................ ] 70% | 11.62 sec | comodulogram: DAR(10, 1) [............................ ] 72% | 11.92 sec | comodulogram: DAR(10, 1) [............................. ] 74% | 12.22 sec | comodulogram: DAR(10, 1) [.............................. ] 76% | 12.53 sec | comodulogram: DAR(10, 1) [............................... ] 78% | 12.83 sec | comodulogram: DAR(10, 1) [................................ ] 80% | 13.14 sec | comodulogram: DAR(10, 1) [................................ ] 82% | 13.44 sec | comodulogram: DAR(10, 1) [................................. ] 84% | 13.74 sec | comodulogram: DAR(10, 1) [.................................. ] 86% | 14.05 sec | comodulogram: DAR(10, 1) [................................... ] 88% | 14.36 sec | comodulogram: DAR(10, 1) [.................................... ] 90% | 14.66 sec | comodulogram: DAR(10, 1) [.................................... ] 92% | 14.96 sec | comodulogram: DAR(10, 1) [..................................... ] 94% | 15.24 sec | comodulogram: DAR(10, 1) [...................................... ] 96% | 15.54 sec | comodulogram: DAR(10, 1) [....................................... ] 98% | 15.85 sec | comodulogram: DAR(10, 1) [........................................] 100% | 16.16 sec | comodulogram: DAR(10, 1) [........................................] 100% | 16.16 sec | comodulogram: DAR(10, 1) .. GENERATED FROM PYTHON SOURCE LINES 70-71 We fit a CSC model on the data. .. GENERATED FROM PYTHON SOURCE LINES 71-88 .. code-block:: Python from alphacsc import learn_d_z params = dict( n_atoms=3, n_times_atom=int(sfreq * 1.0), # 1000. ms reg=5., n_iter=10, solver_z='l-bfgs', solver_z_kwargs=dict(factr=1e9), solver_d_kwargs=dict(factr=1e2), random_state=42, n_jobs=5, verbose=1) _, _, d_hat, z_hat, _ = learn_d_z(data, **params) .. rst-class:: sphx-glr-script-out .. code-block:: none V_0/10 /github/workspace/alphacsc/update_d.py:224: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) lhs_c_copy[0] += lambd /github/workspace/alphacsc/update_d.py:172: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) lambd_hats[k] = lambd_hat ......... .. GENERATED FROM PYTHON SOURCE LINES 89-91 Plot the temporal patterns. Interestingly, we obtain prototypical waveforms of the signal on which we can clearly see the CFC. .. GENERATED FROM PYTHON SOURCE LINES 91-109 .. code-block:: Python n_atoms, n_times_atom = d_hat.shape n_columns = min(6, n_atoms) n_rows = int(np.ceil(n_atoms // n_columns)) figsize = (4 * n_columns, 3 * n_rows) fig, axes = plt.subplots(n_rows, n_columns, figsize=figsize, sharey=True) axes = axes.ravel() for kk in range(n_atoms): ax = axes[kk] time = np.arange(n_times_atom) / sfreq ax.plot(time, d_hat[kk], color='C%d' % kk) ax.set_xlim(0, n_times_atom / sfreq) ax.set(xlabel='Time (sec)', title="Temporal pattern %d" % kk) ax.grid(True) fig.tight_layout() plt.show() .. image-sg:: /auto_examples/csc/images/sphx_glr_plot_cross_frequency_coupling_002.png :alt: Temporal pattern 0, Temporal pattern 1, Temporal pattern 2 :srcset: /auto_examples/csc/images/sphx_glr_plot_cross_frequency_coupling_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (3 minutes 29.722 seconds) .. _sphx_glr_download_auto_examples_csc_plot_cross_frequency_coupling.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_cross_frequency_coupling.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_cross_frequency_coupling.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_cross_frequency_coupling.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_