Skip to content

Instantly share code, notes, and snippets.

@dengemann
Last active March 19, 2016 01:35
Show Gist options
  • Save dengemann/1d01871aae2effd1397f to your computer and use it in GitHub Desktop.
Save dengemann/1d01871aae2effd1397f to your computer and use it in GitHub Desktop.
import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt
import mne
from mne import io
from mne.connectivity import spectral_connectivity
from mne.datasets import sample
from swish.surrogates import theilerize_raw
from mne.minimum_norm import (apply_inverse, apply_inverse_epochs,
read_inverse_operator)
from mne.connectivity import seed_target_indices, spectral_connectivity
###############################################################################
# Set parameters
data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
print(__doc__)
data_path = sample.data_path()
subjects_dir = data_path + '/subjects'
fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
fname_raw = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_event = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
label_name_lh = 'Aud-lh'
fname_label_lh = data_path + '/MEG/sample/labels/%s.label' % label_name_lh
fwd_fname = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)
event_id, tmin, tmax = 1, -0.2, 0.5
method = "dSPM" # use dSPM method (could also be MNE or sLORETA)
# Load data
inverse_operator = read_inverse_operator(fname_inv)
# Setup for reading the raw data
raw = io.Raw(raw_fname, preload=True)
# for pick in picks:
# Add a bad channel
raw.info['bads'] += ['MEG 2443']
picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=False,
exclude='bads')
# rng = np.random.RandomState(42)
# raw._data[picks] = raw._data[rng.permutation(picks)]
# randomize_phase_raw(raw)
# theilerize_raw(raw, random_seed=1, pi_factor=200)
events = mne.read_events(event_fname)
# Pick MEG gradiometers
picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True,
exclude='bads')
# Create epochs for the visual condition
event_id, tmin, tmax = 3, -0.2, 0.5
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6))
stc_gen = apply_inverse_epochs(
epochs, inverse_operator=inverse_operator, return_generator=True)
new_epochs = epochs.copy() # XXX hack
for ii, stc in enumerate(stc_gen):
# new_epochs.append(
# XXX insert theilerization here
new_epochs.data[ii] = mne.apply_forward(
info=epochs.info, forward=fwd, stc=stc).data
# epochs.average().plot_joint()
# stop
# Compute connectivity for band containing the evoked response.
# We exclude the baseline period
fmin, fmax = 10., 20.
sfreq = raw.info['sfreq'] # the sampling frequency
tmin = 0.0 # exclude the baseline period
con, freqs, times, n_epochs, n_tapers = spectral_connectivity(
epochs, method='pli', mode='multitaper', sfreq=sfreq, fmin=fmin, fmax=fmax,
faverage=True, tmin=tmin, mt_adaptive=False, n_jobs=1)
# the epochs contain an EOG channel, which we remove now
ch_names = epochs.ch_names
idx = [ch_names.index(name) for name in ch_names if name.startswith('MEG')]
con = con[idx][:, idx]
# con is a 3D array where the last dimension is size one since we averaged
# over frequencies in a single band. Here we make it 2D
con = con[:, :, 0]
show_con = con + con.T
plt.matshow(show_con, cmap='viridis')
plt.colorbar()
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment