Skip to content

Instantly share code, notes, and snippets.

@dokato
Created June 8, 2018 15:59
Show Gist options
  • Save dokato/84298cd25983288ad16df821034d8329 to your computer and use it in GitHub Desktop.
Save dokato/84298cd25983288ad16df821034d8329 to your computer and use it in GitHub Desktop.
import os
import mne
from mne.datasets.brainstorm import bst_resting
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as ss
########################## INPUT:
folder = os.path.join(bst_resting.data_path(), 'MEG/bst_resting')
file = 'subj002_spontaneous_20111102_01_AUX_raw.fif'
subjects_dir = os.path.join(bst_resting.data_path(), 'subjects/')
subject = 'bst_resting'
trans = None
#################################
raw = mne.io.read_raw_fif(os.path.join(folder, file),
preload=True)
# According to http://neuroimage.usc.edu/brainstorm/DatasetResting
raw.set_channel_types({'EEG057': 'ecg', 'EEG058': 'eog'})
new_fs = 256
raw.resample(new_fs, npad='auto')
channel_picks = mne.pick_types(raw.info, meg=True)
# Filtering the data
raw_filtered = raw.filter(1.5, None, n_jobs=8, picks = channel_picks,
filter_length='auto', fir_design='firwin', method='fir')
raw_filtered = raw_filtered.notch_filter(60, picks=channel_picks, filter_length='auto',
phase='zero')
raw_filtered = raw_filtered.notch_filter(120, picks=channel_picks, filter_length='auto',
phase='zero')
raw_filtered.plot_psd()
del raw
epoch_length = 2
reject = dict(mag=4e-12)
events = mne.event.make_fixed_length_events(raw_filtered, 1, duration = epoch_length)
epochs = mne.Epochs(raw_filtered, events, 1, 0, epoch_length,
proj=True, reject=reject)
rs_cov = mne.compute_covariance( epochs, method = ['shrunk', 'empirical'])
src = mne.setup_source_space(subject, spacing = 'oct6',
subjects_dir=subjects_dir, add_dist=False)
model = mne.make_bem_model(subject=subject, ico=4,
conductivity=(0.3,),
subjects_dir=subjects_dir)
bem = mne.make_bem_solution(model)
################### Warning: no trans matrix!
fwd = mne.make_forward_solution(raw_filtered.info, trans=trans, src=src, bem=bem,
meg=True, eeg=False, mindist=5., n_jobs=2)
###################
# Converts forward solution between different source orientations.
fwd = mne.convert_forward_solution(fwd, surf_ori = True)
fwd = mne.pick_types_forward(fwd, meg=True, eeg=False)
raw_filtered._data = np.random.randn(*raw_filtered._data.shape)
loose=0.2
depth=0.8
inverse_operator = mne.minimum_norm.make_inverse_operator(raw_filtered.info, fwd, rs_cov,
loose=loose, depth=depth)
snr = 1.
lambda2 = 1. / snr ** 2
stc_lor = mne.minimum_norm.apply_inverse_raw(raw_filtered, inverse_operator, lambda2,
method="sLORETA", pick_ori="normal")
stc_beam = mne.beamformer.lcmv_raw(raw_filtered, fwd, None, rs_cov, reg=0.05,
pick_ori="normal", weight_norm='unit-noise-gain',
max_ori_out='signed')
labels = mne.read_labels_from_annot(subject, parc='aparc',
subjects_dir=subjects_dir)
label_colors = [label.color for label in labels]
label_ts_lor = mne.extract_label_time_course(stc_lor, labels, inverse_operator['src'], mode='pca_flip',
return_generator=False)
label_ts_beam = mne.extract_label_time_course(stc_beam, labels, inverse_operator['src'], mode='pca_flip',
return_generator=False)
psds_l, freqs = mne.time_frequency.psd_array_welch(label_ts_lor, raw_filtered.info['sfreq'])
psds_b, freqs = mne.time_frequency.psd_array_welch(label_ts_beam, raw_filtered.info['sfreq'])
def plot_mne_psd(freqs, psds, title = 'PSD'):
f, ax = plt.subplots()
psds = 10 * np.log10(psds)
psds_mean = psds.mean(0)
psds_std = psds.std(0)
ax.plot(freqs, psds_mean, color='k')
ax.fill_between(freqs, psds_mean - psds_std, psds_mean + psds_std,
color='k', alpha=.5)
ax.set(title=title, xlabel='Frequency',
ylabel='Power Spectral Density (dB)')
plt.show()
plot_mne_psd(freqs, psds_l,'sLORETA PSD')
plot_mne_psd(freqs, psds_b,'Beamforming PSD')
def make_env_corr(data):
"""
data - (channels x samples) array
"""
hdata = np.abs(ss.hilbert(data))
return np.corrcoef(hdata)
crrm = make_env_corr(label_ts_beam)
plt.matshow(crrm)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment