dics_sim.py

import os.path as op
import numpy as np
import mne
from mne.time_frequency import csd_morlet
from mne.beamformer import make_dics, apply_dics_csd

# Simulate bilateral auditory data
subject = 'sample'
data_path = mne.datasets.testing.data_path()
fname_raw = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif')
fname_fwd = op.join(
    data_path, 'MEG', 'sample',
    'sample_audvis_trunc-meg-eeg-oct-6-fwd.fif'
)
subjects_dir = op.join(data_path, 'subjects')
info = mne.io.read_raw_fif(fname_raw).del_proj().pick('grad').info
labels = mne.read_labels_from_annot(
    subject, 'aparc.a2009s', 'both', regexp='G_temp_sup-G_T',
    subjects_dir=subjects_dir)
assert len(labels) == 2
duration = 10
assert isinstance(duration, int)
fwd = mne.read_forward_solution(fname_fwd)
sim = mne.simulation.SourceSimulator(fwd['src'], 1. / info['sfreq'], duration)
n = int(round(info['sfreq']))
tot = n * duration
f = 10.
t = np.arange(n // 2) / info['sfreq']  # half the window
samps = np.arange(0, tot, n)
assert len(samps) == duration
events = np.zeros((len(samps), 3), int)
events[:, 2] = 1
events[:, 0] = samps
for li, label in enumerate(labels):
    # introduce a phase shift to make LCMV happy
    ph = li * np.pi / 4.
    data = np.sin(2 * np.pi * f * t + ph)
    data *= np.hanning(data.shape[-1]) * 1e-9
    sim.add_data(label, data, events)
raw = mne.simulation.simulate_raw(info, sim, forward=fwd)
cov = mne.make_ad_hoc_cov(info)
mne.simulation.add_noise(raw, cov)
tmin, tmax = -0.2, 0.5
epochs = mne.Epochs(
    raw, events, event_id=1, tmin=tmin, tmax=tmax, baseline=None, preload=True)
assert len(epochs) == duration - 1
evoked = epochs.average()
evoked.plot_joint()

# DICS
freqs = [f]
tmin, tmax = -0.2, 0.5
csd = csd_morlet(epochs, freqs, tmin=tmin, tmax=tmax, verbose=True)
csd_baseline = csd_morlet(epochs, freqs, tmin=tmin, tmax=0.05, verbose=True)
csd_ers = csd_morlet(epochs, freqs, tmin=0, tmax=tmax, verbose=True)
info = epochs.info
csd = csd.mean()
csd_baseline = csd_baseline.mean()
csd_ers = csd_ers.mean()
filters = make_dics(info, fwd, csd, noise_csd=csd_baseline, reduce_rank=True,
                    pick_ori='max-power', verbose=True)
baseline_source_power, freqs = apply_dics_csd(csd_baseline, filters)
beta_source_power, freqs = apply_dics_csd(csd_ers, filters)
stc_dics = beta_source_power / baseline_source_power

# dSPM
method = 'dSPM'
inv = mne.minimum_norm.make_inverse_operator(info, fwd, cov, verbose=True)
stc_mne = mne.minimum_norm.apply_inverse(evoked, inv, method=method)

# LCMV
cov = mne.compute_covariance(epochs, tmin=0, tmax=tmax)
lcmv = mne.beamformer.make_lcmv(info, forward=fwd, data_cov=cov, verbose=True)
stc_lcmv = mne.beamformer.apply_lcmv(evoked, lcmv)

kwargs = dict(hemi='split', views=['lat', 'med'],
              subjects_dir=subjects_dir)
brain_dics = stc_dics.plot(title='DICS', **kwargs)
brain_mne = stc_mne.plot(title=method, initial_time=0.260, **kwargs)
brain_lcmv = stc_lcmv.plot(title='LCMV', initial_time=0.240, **kwargs)