Compute source power using DICS beamfomer

Compute a Dynamic Imaging of Coherent Sources (DICS) [1] filter from single-trial activity to estimate source power for two frequencies of interest.

References

[1]Gross et al. Dynamic imaging of coherent sources: Studying neural interactions in the human brain. PNAS (2001) vol. 98 (2) pp. 694-699
# Author: Roman Goj <roman.goj@gmail.com>
#         Denis Engemann <denis.engemann@gmail.com>
#
# License: BSD (3-clause)

import mne
from mne.datasets import sample
from mne.time_frequency import csd_epochs
from mne.beamformer import dics_source_power

print(__doc__)

data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
event_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'
fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
subjects_dir = data_path + '/subjects'

Read raw data

raw = mne.io.read_raw_fif(raw_fname)
raw.info['bads'] = ['MEG 2443']  # 1 bad MEG channel

# Set picks
picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False,
                       stim=False, exclude='bads')

# Read epochs
event_id, tmin, tmax = 1, -0.2, 0.5
events = mne.read_events(event_fname)
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=picks, baseline=(None, 0), preload=True,
                    reject=dict(grad=4000e-13, mag=4e-12))
evoked = epochs.average()

# Read forward operator
forward = mne.read_forward_solution(fname_fwd, surf_ori=True)

# Computing the data and noise cross-spectral density matrices
# The time-frequency window was chosen on the basis of spectrograms from
# example time_frequency/plot_time_frequency.py
# As fsum is False csd_epochs returns a list of CrossSpectralDensity
# instances than can then be passed to dics_source_power
data_csds = csd_epochs(epochs, mode='multitaper', tmin=0.04, tmax=0.15,
                       fmin=15, fmax=30, fsum=False)
noise_csds = csd_epochs(epochs, mode='multitaper', tmin=-0.11,
                        tmax=-0.001, fmin=15, fmax=30, fsum=False)

# Compute DICS spatial filter and estimate source power
stc = dics_source_power(epochs.info, forward, noise_csds, data_csds)

for i, csd in enumerate(data_csds):
    message = 'DICS source power at %0.1f Hz' % csd.frequencies[0]
    brain = stc.plot(surface='inflated', hemi='rh', subjects_dir=subjects_dir,
                     time_label=message, figure=i)
    brain.set_data_time_index(i)
    brain.show_view('lateral')
    # Uncomment line below to save images
    # brain.save_image('DICS_source_power_freq_%d.png' % csd.frequencies[0])
  • ../../_images/sphx_glr_plot_dics_source_power_001.png
  • ../../_images/sphx_glr_plot_dics_source_power_002.png

Out:

Opening raw data file /home/ubuntu/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.
Current compensation grade : 0
72 matching events found
Created an SSP operator (subspace dimension = 3)
3 projection items activated
Loading data for 72 events and 421 original time points ...
    Rejecting  epoch based on MAG : [u'MEG 1711']
1 bad epochs dropped
Reading forward solution from /home/ubuntu/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif...
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Desired named matrix (kind = 3523) not available
    Read MEG forward solution (7498 sources, 306 channels, free orientations)
    Desired named matrix (kind = 3523) not available
    Read EEG forward solution (7498 sources, 60 channels, free orientations)
    MEG and EEG forward solutions combined
    Source spaces transformed to the forward solution coordinate frame
    Converting to surface-based source orientations...
    Average patch normals will be employed in the rotation to the local surface coordinates....
[done]
Computing cross-spectral density from epochs...
    using multitaper spectrum estimation with 3 DPSS windows
[done]
Computing cross-spectral density from epochs...
    using multitaper spectrum estimation with 3 DPSS windows
[done]
    305 out of 366 channels remain after picking
Computing DICS source power...
    computing DICS spatial filter 1 out of 2
    computing DICS spatial filter 2 out of 2
[done]
Updating smoothing matrix, be patient..
Smoothing matrix creation, step 1
Smoothing matrix creation, step 2
Smoothing matrix creation, step 3
Smoothing matrix creation, step 4
Smoothing matrix creation, step 5
Smoothing matrix creation, step 6
Smoothing matrix creation, step 7
Smoothing matrix creation, step 8
Smoothing matrix creation, step 9
Smoothing matrix creation, step 10
colormap: fmin=1.53e+00 fmid=1.61e+00 fmax=2.15e+00 transparent=1
Updating smoothing matrix, be patient..
Smoothing matrix creation, step 1
Smoothing matrix creation, step 2
Smoothing matrix creation, step 3
Smoothing matrix creation, step 4
Smoothing matrix creation, step 5
Smoothing matrix creation, step 6
Smoothing matrix creation, step 7
Smoothing matrix creation, step 8
Smoothing matrix creation, step 9
Smoothing matrix creation, step 10
colormap: fmin=1.53e+00 fmid=1.61e+00 fmax=2.15e+00 transparent=1

Total running time of the script: ( 0 minutes 58.787 seconds)

Generated by Sphinx-Gallery