Plot a vector time-frequency source time course

vector_stc.py

# -*- coding: utf-8 -*-
"""
.. _ex-inverse-dics-epochs:
=======================================================================
Compute source level time-frequency timecourses using a DICS beamformer
=======================================================================
In this example, a Dynamic Imaging of Coherent Sources (DICS)
:footcite:`GrossEtAl2001` beamformer is used to transform sensor-level
time-frequency objects to the source level. We will look at the event-related
synchronization (ERS) of beta band activity in the :ref:`somato dataset
<somato-dataset>`.
"""
# Authors: Marijn van Vliet <w.m.vanvliet@gmail.com>
#          Alex Rockhill <aprockhill@mailbox.org>
#
# License: BSD-3-Clause

import numpy as np
import matplotlib.pyplot as plt
import mne
from mne.datasets import somato
from mne.time_frequency import tfr_morlet, csd_tfr
from mne.beamformer import make_dics, apply_dics_tfr_epochs

print(__doc__)

# %%
# Organize the data that we will use for this example.

data_path = somato.data_path()
subject = '01'
task = 'somato'
raw_fname = (data_path / f'sub-{subject}' / 'meg' /
             f'sub-{subject}_task-{task}_meg.fif')
fname_fwd = (data_path / 'derivatives' / f'sub-{subject}' /
             f'sub-{subject}_task-{task}-fwd.fif')
subjects_dir = data_path / 'derivatives' / 'freesurfer' / 'subjects'

# %%
# First, we load the data and compute for each epoch the time-frequency
# decomposition in sensor space.

# Load raw data and make epochs. For speed, we only use the first 5 events.
raw = mne.io.read_raw_fif(raw_fname)
events = mne.find_events(raw)[:3]
epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2,
                    preload=True)

# We are mostly interested in the beta band since it has been shown to be
# active for somatosensory stimulation
freqs = np.array([15])

# Use Morlet wavelets to compute sensor-level time-frequency (TFR)
# decomposition for each epoch. We must pass ``output='complex'`` if we wish to
# use this TFR later with a DICS beamformer. We also pass ``average=False`` to
# compute the TFR for each individual epoch.
epochs_tfr = tfr_morlet(epochs, freqs, n_cycles=5, return_itc=False,
                        output='complex', average=False)

# %%
# Now, we build a DICS beamformer and project the sensor-level TFR to the
# source level.

# Compute the Cross-Spectral Density (CSD) matrix for the sensor-level TFRs.
# We are interested in increases in power relative to the baseline period, so
# we will make a separate CSD for just that period as well.
csd = csd_tfr(epochs_tfr, tmin=-1, tmax=1.5)
baseline_csd = csd_tfr(epochs_tfr, tmin=-1, tmax=0)

# use the CSDs and the forward model to build the DICS beamformer
fwd = mne.read_forward_solution(fname_fwd)

# compute vector solution
filters = make_dics(epochs.info, fwd, csd, noise_csd=baseline_csd,
                    pick_ori='vector', reduce_rank=True, real_filter=False)


# remove edge artifact
epochs_tfr.crop(tmin=-1, tmax=1.5)

# project the TFR for each epoch to source space
epochs_stcs = apply_dics_tfr_epochs(
    epochs_tfr, filters, return_generator=True)

# %%
# We can view the vector solution for each time-frequency source time course.
# We have a complex solution for the time-frequency; the real portion for
# the sine wave fit and the imaginary portion from the cosine fit. If we take
# the angle, that gives us the phase. Since, we have three independent
# solutions, one for ``x``, one for ``y`` and one for ``z`` we can take the
# cosine of the angle to find when the solution is pointing in the direction
# of each component. For instance, when the phase of the ``x`` component is
# near 0 or 180 degrees, a large positive and negative vector are mapped to
# the ``x`` direction respectively. When the phase is near 90 or 270 degrees,
# the vector is mapped to near zero. If there is power at that time and at
# this chosen vertex, it must be going in the ``y`` or ``z`` direction. We
# don't know which it is, so it's a good thing we have the two other ``y``
# and ``z`` independent source estimates. The complimentary information from
# all three directions, will give us a reasonable estimate of the vector
# and, consequently, not only the strength but the direction of the dipole.

# select one stc for an example
stc = next(epochs_stcs)

stc.data = stc.data.real  # np.abs(stc.data) * np.cos(np.angle(stc.data))

# plot the timecourse phase
brain = stc.plot(
    subjects_dir=subjects_dir,
    hemi='both',
    views='dorsal',
    initial_time=0.875,
)

# simulate data:
# 1) purely x component

stc.data = np.zeros(stc.data.shape, dtype=complex)
t = np.arange(stc.times.size, dtype=np.float64) / epochs.info['sfreq']
# 5 Hz simulate pure sine wave
stc.data[:, 0] = np.sin(np.pi * 2 * 5 * t) + \
  1j * np.sin(np.pi * 2. * 5 * (t + np.pi / 2))
stc.data = stc.data.real

# plot the timecourse phase
brain = stc.plot(
    subjects_dir=subjects_dir,
    hemi='both',
    views='dorsal',
    clim=dict(kind='percent', lims=(0, 50, 100))
)

# 2) oblique component

stc.data = np.zeros(stc.data.shape, dtype=complex)
t = np.arange(stc.times.size, dtype=np.float64) / epochs.info['sfreq']
# 5 Hz simulate pure sine wave
stc.data[:, :] = np.sin(np.pi * 2 * 5 * t) + \
  1j * np.sin(np.pi * 2. * 5 * (t + np.pi / 2))
stc.data[:, 0] *= 2  # even more oblique
stc.data = stc.data.real

# plot the timecourse phase
brain = stc.plot(
    subjects_dir=subjects_dir,
    hemi='both',
    views='dorsal',
    clim=dict(kind='percent', lims=(0, 50, 100))
)


# You can save a movie like the one on our documentation website with:
# brain.save_movie(time_dilation=30, tmin=0.875, tmax=0.95,
#                  interpolation='linear', time_viewer=True)