resting_state_inverse2.py

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()