Day 23: Beamformers

Today I tested the performance of two beamformer methods on simulated data. In order to do this, I had to generate
simulated epoch data, as opposed to simulated evoked (epochs are the trial by trial data, while evoked is the
average of the epochs). Simulating the epoch data was similar, except I used a lower signal to noise ratio (for
simulation code, see my github repository).

Compared to yesterday, the performance looks worse.

The next steps might be to look at the whole brain volume results, and to set up a script to loop through all
parameters and find the set that performs the best.

In [1]:

%matplotlib inline
%load spm_test_beamformers.py

In [3]:

import numpy as np
import matplotlib.pyplot as plt
import mne
from mne.datasets import spm_face

mne.set_log_level(False)

data_path = spm_face.data_path()

fwd_fname = data_path + \
            '/MEG/spm/SPM_CTF_MEG_example_faces1_3D_oct-6-rh-amyg-vol-fwd.fif'

epo_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D-sim-epo.fif'

fwd = mne.read_forward_solution(fwd_fname)
sim = mne.read_epochs(epo_fname)

noise_cov = mne.compute_covariance(sim.crop(None, 0, copy=True))
data_cov = mne.compute_covariance(sim.crop(0, None, copy=True))

noise_csd = mne.time_frequency.compute_epochs_csd(sim.crop(None, 0, copy=True))
data_csd = mne.time_frequency.compute_epochs_csd(sim.crop(0, None, copy=True))

stcs = []
stc = mne.beamformer.lcmv(sim.average(), fwd, noise_cov, data_cov)
stcs.append(stc)

stc = mne.beamformer.dics(sim.average(), fwd, noise_csd, data_csd)
stcs.append(stc)

for stc in stcs:

    fig = plt.figure()
    t = stc.times

    ax = fig.add_subplot(221)
    ax.plot(t, stc.data[len(stc.vertno[0])+len(stc.vertno[1]):].mean(0))
    ax.set_title('Right Amygdala')

    ax = fig.add_subplot(222)
    ax.plot(t, stc.lh_data.mean(0))
    ax.set_title('Left Hemisphere')

    ax = fig.add_subplot(223)
    ax.plot(t, stc.rh_data.mean(0))
    ax.set_title('Right Hemisphere')

    ax = fig.add_subplot(224)
    x = np.zeros((2, len(t)))
    ws = mne.time_frequency.morlet(sim.info['sfreq'], [3, 10], n_cycles=[1, 1.5])
    x[0][:len(ws[0])] = np.real(ws[0])
    x[1][:len(ws[1])] = np.real(ws[1])
    x[0] = np.roll(x[0], 100)
    x[1] = np.roll(x[1], 200)
    ax.plot(t, x[0])
    ax.plot(t, x[1])

plt.show()