Skip to content

Instantly share code, notes, and snippets.

@dengemann
Last active October 29, 2017 18:06
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save dengemann/511923a3044c2c0562df14f742516d1b to your computer and use it in GitHub Desktop.
Save dengemann/511923a3044c2c0562df14f742516d1b to your computer and use it in GitHub Desktop.
# Authors: Denis A. Engemann <denis.engemann@gmail.com>
#
# License: BSD (3-clause)
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
import mne
data_path = mne.datasets.somato.data_path()
raw_fname = data_path + '/MEG/somato/sef_raw_sss.fif'
event_id, tmin, tmax = 1, -1., 3.
# Setup for reading the raw data
raw = mne.io.Raw(raw_fname, preload=True)
raw.filter(1, 30, method='iir')
baseline = (None, 0)
events = mne.find_events(raw, stim_channel='STI 014')
# picks MEG
picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=True, stim=False)
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
baseline=baseline, reject=dict(mag=4e-12, grad=4000e-13),
preload=True)
scalings = dict(mag=4e15, grad=1e13)
noise_cov = mne.cov.compute_covariance(epochs, tmax=0)
evoked = epochs.average()
picks = mne.pick_types(evoked.info, meg=True)
picks_mag = mne.pick_types(evoked.info, meg='mag')
picks_grad = mne.pick_types(evoked.info, meg='grad')
# Delete MAG-Grad Covariance (Xterms)
noise_cov_b = noise_cov.copy()
noise_cov_b['data'][np.ix_(picks_mag, picks_grad)] = 0.0
noise_cov_b['data'][np.ix_(picks_grad, picks_mag)] = 0.0
# plot
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
evoked_white = mne.cov.whiten_evoked(evoked, noise_cov, picks)
evoked_white.plot(picks=picks, unit=False, hline=[-2, 2], axes=axes[:, 0])
evoked_white = mne.cov.whiten_evoked(evoked, noise_cov_b, picks)
evoked_white.plot(picks=picks, unit=False, hline=[-2, 2], axes=axes[:, 1])
axes[0, 0].set_title(axes[0, 0].get_title() + ' + Xterms')
axes[0, 1].set_title(axes[0, 1].get_title() + ' no Xterms')
axes[1, 0].set_title(axes[1, 0].get_title() + ' + Xterms')
axes[1, 1].set_title(axes[1, 1].get_title() + ' no Xterms')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment