0.16.1
  • Install
  • Documentation
  • API
  • Examples
  • Contribute
  • Site
    • Development
    • v0.16 (stable)
    • v0.15
    • v0.14
    • v0.13
    • v0.12
    • v0.11

    Logo

    • Explore event-related dynamics for specific frequency bands
      • References

    Note

    Click here to download the full example code

    Explore event-related dynamics for specific frequency bands¶

    The objective is to show you how to explore spectrally localized effects. For this purpose we adapt the method described in [1] and use it on the somato dataset. The idea is to track the band-limited temporal evolution of spatial patterns by using the Global Field Power (GFP).

    We first bandpass filter the signals and then apply a Hilbert transform. To reveal oscillatory activity the evoked response is then subtracted from every single trial. Finally, we rectify the signals prior to averaging across trials by taking the magniude of the Hilbert. Then the GFP is computed as described in [2], using the sum of the squares but without normalization by the rank. Baselining is subsequently applied to make the GFPs comparable between frequencies. The procedure is then repeated for each frequency band of interest and all GFPs are visualized. To estimate uncertainty, non-parametric confidence intervals are computed as described in [3] across channels.

    The advantage of this method over summarizing the Space x Time x Frequency output of a Morlet Wavelet in frequency bands is relative speed and, more importantly, the clear-cut comparability of the spectral decomposition (the same type of filter is used across all bands).

    References¶

    [1]Hari R. and Salmelin R. Human cortical oscillations: a neuromagnetic view through the skull (1997). Trends in Neuroscience 20 (1), pp. 44-49.
    [2]Engemann D. and Gramfort A. (2015) Automated model selection in covariance estimation and spatial whitening of MEG and EEG signals, vol. 108, 328-342, NeuroImage.
    [3]Efron B. and Hastie T. Computer Age Statistical Inference (2016). Cambrdige University Press, Chapter 11.2.
    # Authors: Denis A. Engemann <denis.engemann@gmail.com>
    #
    # License: BSD (3-clause)
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    import mne
    from mne.datasets import somato
    from mne.baseline import rescale
    from mne.stats import _bootstrap_ci
    

    Set parameters

    data_path = somato.data_path()
    raw_fname = data_path + '/MEG/somato/sef_raw_sss.fif'
    
    # let's explore some frequency bands
    iter_freqs = [
        ('Theta', 4, 7),
        ('Alpha', 8, 12),
        ('Beta', 13, 25),
        ('Gamma', 30, 45)
    ]
    

    We create average power time courses for each frequency band

    # set epoching parameters
    event_id, tmin, tmax = 1, -1., 3.
    baseline = None
    
    # get the header to extract events
    raw = mne.io.read_raw_fif(raw_fname, preload=False)
    events = mne.find_events(raw, stim_channel='STI 014')
    
    frequency_map = list()
    
    for band, fmin, fmax in iter_freqs:
        # (re)load the data to save memory
        raw = mne.io.read_raw_fif(raw_fname, preload=True)
        raw.pick_types(meg='grad', eog=True)  # we just look at gradiometers
    
        # bandpass filter and compute Hilbert
        raw.filter(fmin, fmax, n_jobs=1,  # use more jobs to speed up.
                   l_trans_bandwidth=1,  # make sure filter params are the same
                   h_trans_bandwidth=1,  # in each band and skip "auto" option.
                   fir_design='firwin')
        raw.apply_hilbert(n_jobs=1, envelope=False)
    
        epochs = mne.Epochs(raw, events, event_id, tmin, tmax, baseline=baseline,
                            reject=dict(grad=4000e-13, eog=350e-6), preload=True)
        # remove evoked response and get analytic signal (envelope)
        epochs.subtract_evoked()  # for this we need to construct new epochs.
        epochs = mne.EpochsArray(
            data=np.abs(epochs.get_data()), info=epochs.info, tmin=epochs.tmin)
        # now average and move on
        frequency_map.append(((band, fmin, fmax), epochs.average()))
    

    Out:

    Opening raw data file /home/circleci/mne_data/MNE-somato-data/MEG/somato/sef_raw_sss.fif...
        Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
    Ready.
    Current compensation grade : 0
    111 events found
    Event IDs: [1]
    Opening raw data file /home/circleci/mne_data/MNE-somato-data/MEG/somato/sef_raw_sss.fif...
        Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
    Ready.
    Current compensation grade : 0
    Reading 0 ... 269399  =      0.000 ...   897.077 secs...
    Setting up band-pass filter from 4 - 7 Hz
    Filter length of 991 samples (3.300 sec) selected
    111 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    Loading data for 111 events and 1202 original time points ...
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
    3 bad epochs dropped
    Subtracting Evoked from Epochs
        The following channels are not included in the subtraction: EOG 061
    [done]
    108 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    0 bad epochs dropped
    Opening raw data file /home/circleci/mne_data/MNE-somato-data/MEG/somato/sef_raw_sss.fif...
        Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
    Ready.
    Current compensation grade : 0
    Reading 0 ... 269399  =      0.000 ...   897.077 secs...
    Setting up band-pass filter from 8 - 12 Hz
    Filter length of 991 samples (3.300 sec) selected
    111 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    Loading data for 111 events and 1202 original time points ...
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
    3 bad epochs dropped
    Subtracting Evoked from Epochs
        The following channels are not included in the subtraction: EOG 061
    [done]
    108 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    0 bad epochs dropped
    Opening raw data file /home/circleci/mne_data/MNE-somato-data/MEG/somato/sef_raw_sss.fif...
        Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
    Ready.
    Current compensation grade : 0
    Reading 0 ... 269399  =      0.000 ...   897.077 secs...
    Setting up band-pass filter from 13 - 25 Hz
    Filter length of 991 samples (3.300 sec) selected
    111 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    Loading data for 111 events and 1202 original time points ...
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
    3 bad epochs dropped
    Subtracting Evoked from Epochs
        The following channels are not included in the subtraction: EOG 061
    [done]
    108 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    0 bad epochs dropped
    Opening raw data file /home/circleci/mne_data/MNE-somato-data/MEG/somato/sef_raw_sss.fif...
        Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
    Ready.
    Current compensation grade : 0
    Reading 0 ... 269399  =      0.000 ...   897.077 secs...
    Setting up band-pass filter from 30 - 45 Hz
    Filter length of 991 samples (3.300 sec) selected
    111 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    Loading data for 111 events and 1202 original time points ...
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
        Rejecting  epoch based on EOG : ['EOG 061']
    3 bad epochs dropped
    Subtracting Evoked from Epochs
        The following channels are not included in the subtraction: EOG 061
    [done]
    108 matching events found
    No baseline correction applied
    Not setting metadata
    0 projection items activated
    0 bad epochs dropped
    

    Now we can compute the Global Field Power We can track the emergence of spatial patterns compared to baseline for each frequency band, with a bootstrapped confidence interval.

    We see dominant responses in the Alpha and Beta bands.

    fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True, sharey=True)
    colors = plt.get_cmap('winter_r')(np.linspace(0, 1, 4))
    for ((freq_name, fmin, fmax), average), color, ax in zip(
            frequency_map, colors, axes.ravel()[::-1]):
        times = average.times * 1e3
        gfp = np.sum(average.data ** 2, axis=0)
        gfp = mne.baseline.rescale(gfp, times, baseline=(None, 0))
        ax.plot(times, gfp, label=freq_name, color=color, linewidth=2.5)
        ax.axhline(0, linestyle='--', color='grey', linewidth=2)
        ci_low, ci_up = _bootstrap_ci(average.data, random_state=0,
                                      stat_fun=lambda x: np.sum(x ** 2, axis=0))
        ci_low = rescale(ci_low, average.times, baseline=(None, 0))
        ci_up = rescale(ci_up, average.times, baseline=(None, 0))
        ax.fill_between(times, gfp + ci_up, gfp - ci_low, color=color, alpha=0.3)
        ax.grid(True)
        ax.set_ylabel('GFP')
        ax.annotate('%s (%d-%dHz)' % (freq_name, fmin, fmax),
                    xy=(0.95, 0.8),
                    horizontalalignment='right',
                    xycoords='axes fraction')
        ax.set_xlim(-1000, 3000)
    
    axes.ravel()[-1].set_xlabel('Time [ms]')
    
    ../../_images/sphx_glr_plot_time_frequency_global_field_power_001.png

    Out:

    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    Applying baseline correction (mode: mean)
    

    Total running time of the script: ( 1 minutes 5.359 seconds)

    Download Python source code: plot_time_frequency_global_field_power.py
    Download Jupyter notebook: plot_time_frequency_global_field_power.ipynb

    Gallery generated by Sphinx-Gallery

    Institutions
    • GitHub
    • ·
    • Mailing list
    • ·
    • Gitter
    • ·
    • What's new
    • ·
    • Cite MNE
    • Back to top

    © Copyright 2012-2018, MNE Developers. Last updated on 2018-06-24.