Note
Click here to download the full example code
Maxwell filter data with movement compensation#
Demonstrate movement compensation on simulated data. The simulated data contains bilateral activation of auditory cortices, repeated over 14 different head rotations (head center held fixed). See the following for details:
# Authors: Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD-3-Clause
import mne
from mne.preprocessing import maxwell_filter
print(__doc__)
data_path = mne.datasets.misc.data_path(verbose=True) / 'movement'
head_pos = mne.chpi.read_head_pos(data_path / 'simulated_quats.pos')
raw = mne.io.read_raw_fif(data_path / 'simulated_movement_raw.fif')
raw_stat = mne.io.read_raw_fif(data_path / 'simulated_stationary_raw.fif')
Opening raw data file /home/circleci/mne_data/MNE-misc-data/movement/simulated_movement_raw.fif...
Range : 25800 ... 34208 = 42.956 ... 56.955 secs
Ready.
Opening raw data file /home/circleci/mne_data/MNE-misc-data/movement/simulated_stationary_raw.fif...
Range : 25800 ... 34208 = 42.956 ... 56.955 secs
Ready.
Visualize the “subject” head movements. By providing the measurement information, the distance to the nearest sensor in each direction (e.g., left/right for the X direction, forward/backward for Y) can be shown in blue, and the destination (if given) shown in red.
mne.viz.plot_head_positions(
head_pos, mode='traces', destination=raw.info['dev_head_t'], info=raw.info)
This can also be visualized using a quiver.
mne.viz.plot_head_positions(
head_pos, mode='field', destination=raw.info['dev_head_t'], info=raw.info)
Getting helmet for system 306m
Process our simulated raw data (taking into account head movements).
# extract our resulting events
events = mne.find_events(raw, stim_channel='STI 014')
events[:, 2] = 1
raw.plot(events=events)
topo_kwargs = dict(times=[0, 0.1, 0.2], ch_type='mag', vmin=-500, vmax=500,
time_unit='s')
14 events found
Event IDs: [ 1 2 3 4 5 6 7 8 9 10 11 12 13 14]
First, take the average of stationary data (bilateral auditory patterns).
evoked_stat = mne.Epochs(raw_stat, events, 1, -0.2, 0.8).average()
evoked_stat.plot_topomap(title='Stationary', **topo_kwargs)
Not setting metadata
14 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] sec
Applying baseline correction (mode: mean)
0 projection items activated
Second, take a naive average, which averages across epochs that have been simulated to have different head positions and orientations, thereby spatially smearing the activity.
epochs = mne.Epochs(raw, events, 1, -0.2, 0.8)
evoked = epochs.average()
evoked.plot_topomap(title='Moving: naive average', **topo_kwargs)
Not setting metadata
14 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] sec
Applying baseline correction (mode: mean)
0 projection items activated
Third, use raw movement compensation (restores pattern).
raw_sss = maxwell_filter(raw, head_pos=head_pos)
evoked_raw_mc = mne.Epochs(raw_sss, events, 1, -0.2, 0.8).average()
evoked_raw_mc.plot_topomap(title='Moving: movement compensated (raw)',
**topo_kwargs)
Maxwell filtering raw data
No bad MEG channels
Processing 203 gradiometers and 102 magnetometers
Automatic origin fit: head of radius 91.0 mm
Using origin -4.1, 16.0, 51.7 mm in the head frame
Using 90/95 harmonic components for 0.000 (75/80 in, 15/15 out)
Appending head position result channels and loading raw data from disk
Processing 1 data chunk
Using 87/95 harmonic components for 0.000 (72/80 in, 15/15 out)
Using 88/95 harmonic components for 1.001 (73/80 in, 15/15 out)
Using 90/95 harmonic components for 2.000 (75/80 in, 15/15 out)
Using 88/95 harmonic components for 3.000 (73/80 in, 15/15 out)
Using 88/95 harmonic components for 3.999 (73/80 in, 15/15 out)
Using 88/95 harmonic components for 5.000 (73/80 in, 15/15 out)
Using 89/95 harmonic components for 6.001 (74/80 in, 15/15 out)
Using 93/95 harmonic components for 6.999 (78/80 in, 15/15 out)
Using 88/95 harmonic components for 8.000 (73/80 in, 15/15 out)
Using 91/95 harmonic components for 9.001 (76/80 in, 15/15 out)
Using 93/95 harmonic components for 10.000 (78/80 in, 15/15 out)
Using 93/95 harmonic components for 11.000 (78/80 in, 15/15 out)
Using 89/95 harmonic components for 11.999 (74/80 in, 15/15 out)
Using 88/95 harmonic components for 13.000 (73/80 in, 15/15 out)
Used 14 head positions for 0.000 - 13.999 sec (#1/1)
[done]
Not setting metadata
14 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] sec
Applying baseline correction (mode: mean)
0 projection items activated
Fourth, use evoked movement compensation. For these data, which contain very large rotations, it does not as cleanly restore the pattern.
evoked_evo_mc = mne.epochs.average_movements(epochs, head_pos=head_pos)
evoked_evo_mc.plot_topomap(title='Moving: movement compensated (evoked)',
**topo_kwargs)
Automatic origin fit: head of radius 91.0 mm
Aligning and averaging up to 14 epochs
No bad MEG channels
Processing 203 gradiometers and 102 magnetometers
Processing epoch 1 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 2 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 3 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 4 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 5 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 6 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 7 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 8 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 9 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 10 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 11 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 12 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 13 (device location: 0.0, 0.0, 40.0 mm)
Processing epoch 14 (device location: 0.0, 0.0, 40.0 mm)
Created Evoked dataset from 14 epochs
Total running time of the script: ( 0 minutes 15.532 seconds)
Estimated memory usage: 10 MB