1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
"""
=====================================
Time-frequency beamforming using LCMV
=====================================
Compute LCMV source power [1]_ in a grid of time-frequency windows and
display results.
References
----------
.. [1] Dalal et al. Five-dimensional neuroimaging: Localization of the
time-frequency dynamics of cortical activity.
NeuroImage (2008) vol. 40 (4) pp. 1686-1700
"""
# Author: Roman Goj <roman.goj@gmail.com>
#
# License: BSD (3-clause)
import mne
from mne import compute_covariance
from mne.datasets import sample
from mne.event import make_fixed_length_events
from mne.beamformer import tf_lcmv
from mne.viz import plot_source_spectrogram
print(__doc__)
data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
noise_fname = data_path + '/MEG/sample/ernoise_raw.fif'
event_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'
fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
subjects_dir = data_path + '/subjects'
label_name = 'Aud-lh'
fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
###############################################################################
# Read raw data, preload to allow filtering
raw = mne.io.read_raw_fif(raw_fname, preload=True)
raw.info['bads'] = ['MEG 2443'] # 1 bad MEG channel
# Pick a selection of magnetometer channels. A subset of all channels was used
# to speed up the example. For a solution based on all MEG channels use
# meg=True, selection=None and add grad=4000e-13 to the reject dictionary.
# We could do this with a "picks" argument to Epochs and the LCMV functions,
# but here we use raw.pick_types() to save memory.
left_temporal_channels = mne.read_selection('Left-temporal')
raw.pick_types(meg='mag', eeg=False, eog=False, stim=False, exclude='bads',
selection=left_temporal_channels)
reject = dict(mag=4e-12)
# Re-normalize our empty-room projectors, which should be fine after
# subselection
raw.info.normalize_proj()
# Setting time limits for reading epochs. Note that tmin and tmax are set so
# that time-frequency beamforming will be performed for a wider range of time
# points than will later be displayed on the final spectrogram. This ensures
# that all time bins displayed represent an average of an equal number of time
# windows.
tmin, tmax = -0.5, 0.75 # s
tmin_plot, tmax_plot = -0.3, 0.5 # s
# Read epochs. Note that preload is set to False to enable tf_lcmv to read the
# underlying raw object.
# Filtering is then performed on raw data in tf_lcmv and the epochs
# parameters passed here are used to create epochs from filtered data. However,
# reading epochs without preloading means that bad epoch rejection is delayed
# until later. To perform bad epoch rejection based on the reject parameter
# passed here, run epochs.drop_bad(). This is done automatically in
# tf_lcmv to reject bad epochs based on unfiltered data.
event_id = 1
events = mne.read_events(event_fname)
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
baseline=None, preload=False, reject=reject)
# Read empty room noise, preload to allow filtering, and pick subselection
raw_noise = mne.io.read_raw_fif(noise_fname, preload=True)
raw_noise.info['bads'] = ['MEG 2443'] # 1 bad MEG channel
raw_noise.pick_types(meg='mag', eeg=False, eog=False, stim=False,
exclude='bads', selection=left_temporal_channels)
raw_noise.info.normalize_proj()
# Create artificial events for empty room noise data
events_noise = make_fixed_length_events(raw_noise, event_id, duration=1.)
# Create an epochs object using preload=True to reject bad epochs based on
# unfiltered data
epochs_noise = mne.Epochs(raw_noise, events_noise, event_id, tmin, tmax,
proj=True, baseline=None,
preload=True, reject=reject)
# Make sure the number of noise epochs is the same as data epochs
epochs_noise = epochs_noise[:len(epochs.events)]
# Read forward operator
forward = mne.read_forward_solution(fname_fwd)
# Read label
label = mne.read_label(fname_label)
###############################################################################
# Time-frequency beamforming based on LCMV
# Setting frequency bins as in Dalal et al. 2008 (high gamma was subdivided)
freq_bins = [(4, 12), (12, 30), (30, 55), (65, 299)] # Hz
win_lengths = [0.3, 0.2, 0.15, 0.1] # s
# Setting the time step
tstep = 0.05
# Setting the whitened data covariance regularization parameter
data_reg = 0.001
# Subtract evoked response prior to computation?
subtract_evoked = False
# Calculating covariance from empty room noise. To use baseline data as noise
# substitute raw for raw_noise, epochs.events for epochs_noise.events, tmin for
# desired baseline length, and 0 for tmax_plot.
# Note, if using baseline data, the averaged evoked response in the baseline
# period should be flat.
noise_covs = []
for (l_freq, h_freq) in freq_bins:
raw_band = raw_noise.copy()
raw_band.filter(l_freq, h_freq, n_jobs=1, fir_design='firwin')
epochs_band = mne.Epochs(raw_band, epochs_noise.events, event_id,
tmin=tmin_plot, tmax=tmax_plot, baseline=None,
proj=True)
noise_cov = compute_covariance(epochs_band, method='shrunk', rank=None)
noise_covs.append(noise_cov)
del raw_band # to save memory
# Computing LCMV solutions for time-frequency windows in a label in source
# space for faster computation, use label=None for full solution
stcs = tf_lcmv(epochs, forward, noise_covs, tmin, tmax, tstep, win_lengths,
freq_bins=freq_bins, subtract_evoked=subtract_evoked,
reg=data_reg, label=label, rank=None)
# Plotting source spectrogram for source with maximum activity.
# Note that tmin and tmax are set to display a time range that is smaller than
# the one for which beamforming estimates were calculated. This ensures that
# all time bins shown are a result of smoothing across an identical number of
# time windows.
plot_source_spectrogram(stcs, freq_bins, tmin=tmin_plot, tmax=tmax_plot,
source_index=None, colorbar=True)
|
