非官方中文翻译
本页译自 MNE-Python 官方文档：
查看英文原文 ↗。
文档版权归 MNE Developers 所有（BSD-3-Clause）。
说明
转到末尾下载完整示例代码。
使用 DICS 波束形成器计算源功率#
从单试次活动计算动态相干源成像（DICS）[1] 滤波器，以估计某频段上的源功率。本示例演示如何对 somato 数据集中 beta 频段活动的事件相关同步（ERS）做源定位。
# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
# Roman Goj <roman.goj@gmail.com>
# Denis Engemann <denis.engemann@gmail.com>
# Stefan Appelhoff <stefan.appelhoff@mailbox.org>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import numpy as np
import mne
from mne.beamformer import apply_dics_csd, make_dics
from mne.datasets import somato
from mne.time_frequency import csd_morlet
print(__doc__)
读取原始数据并创建分段：
data_path = somato.data_path()
subject = "01"
task = "somato"
raw_fname = data_path / f"sub-{subject}" / "meg" / f"sub-{subject}_task-{task}_meg.fif"
# Use a shorter segment of raw just for speed here
raw = mne.io.read_raw_fif(raw_fname)
raw.crop(0, 120) # one minute for speed (looks similar to using all ~800 s)
# Read epochs
events = mne.find_events(raw)
epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, preload=True)
del raw
# Paths to forward operator and FreeSurfer subject directory
fname_fwd = (
data_path / "derivatives" / f"sub-{subject}" / f"sub-{subject}_task-{task}-fwd.fif"
)
subjects_dir = data_path / "derivatives" / "freesurfer" / "subjects"
Opening raw data file /home/circleci/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif...
Range : 237600 ... 506999 = 791.189 ... 1688.266 secs
Ready.
Finding events on: STI 014
15 events found on stim channel STI 014
Event IDs: [1]
Not setting metadata
15 matching events found
Setting baseline interval to [-1.498464098687909, 0.0] s
Applying baseline correction (mode: mean)
0 projection items activated
Loading data for 15 events and 1052 original time points ...
0 bad epochs dropped
我们关注 beta 频段：用对数刻度定义 12 到 30 Hz 的频率范围。
freqs = np.logspace(np.log10(12), np.log10(30), 9)
为不同时间区间计算 beta 频段的互谱密度矩阵。为加快计算，本示例使用 decim=20，代价是精度下降。
csd = csd_morlet(epochs, freqs, tmin=-1, tmax=1.5, decim=20)
csd_baseline = csd_morlet(epochs, freqs, tmin=-1, tmax=0, decim=20)
# ERS activity starts at 0.5 seconds after stimulus onset
csd_ers = csd_morlet(epochs, freqs, tmin=0.5, tmax=1.5, decim=20)
info = epochs.info
del epochs
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/15 [00:00<?, ?it/s]
7%|▋ | CSD epoch blocks : 1/15 [00:00<00:00, 15.71it/s]
13%|█▎ | CSD epoch blocks : 2/15 [00:00<00:00, 15.69it/s]
20%|██ | CSD epoch blocks : 3/15 [00:00<00:00, 15.68it/s]
27%|██▋ | CSD epoch blocks : 4/15 [00:00<00:00, 15.66it/s]
33%|███▎ | CSD epoch blocks : 5/15 [00:00<00:00, 15.69it/s]
40%|████ | CSD epoch blocks : 6/15 [00:00<00:00, 15.68it/s]
47%|████▋ | CSD epoch blocks : 7/15 [00:00<00:00, 15.70it/s]
53%|█████▎ | CSD epoch blocks : 8/15 [00:00<00:00, 15.70it/s]
60%|██████ | CSD epoch blocks : 9/15 [00:00<00:00, 15.71it/s]
67%|██████▋ | CSD epoch blocks : 10/15 [00:00<00:00, 15.72it/s]
73%|███████▎ | CSD epoch blocks : 11/15 [00:00<00:00, 15.71it/s]
80%|████████ | CSD epoch blocks : 12/15 [00:00<00:00, 15.68it/s]
87%|████████▋ | CSD epoch blocks : 13/15 [00:00<00:00, 15.67it/s]
93%|█████████▎| CSD epoch blocks : 14/15 [00:00<00:00, 15.66it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:00<00:00, 15.66it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:00<00:00, 15.66it/s]
[done]
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/15 [00:00<?, ?it/s]
7%|▋ | CSD epoch blocks : 1/15 [00:00<00:00, 21.89it/s]
13%|█▎ | CSD epoch blocks : 2/15 [00:00<00:00, 21.85it/s]
20%|██ | CSD epoch blocks : 3/15 [00:00<00:00, 21.87it/s]
27%|██▋ | CSD epoch blocks : 4/15 [00:00<00:00, 21.88it/s]
33%|███▎ | CSD epoch blocks : 5/15 [00:00<00:00, 21.90it/s]
40%|████ | CSD epoch blocks : 6/15 [00:00<00:00, 21.91it/s]
47%|████▋ | CSD epoch blocks : 7/15 [00:00<00:00, 21.89it/s]
53%|█████▎ | CSD epoch blocks : 8/15 [00:00<00:00, 21.91it/s]
60%|██████ | CSD epoch blocks : 9/15 [00:00<00:00, 21.90it/s]
67%|██████▋ | CSD epoch blocks : 10/15 [00:00<00:00, 21.91it/s]
73%|███████▎ | CSD epoch blocks : 11/15 [00:00<00:00, 21.89it/s]
80%|████████ | CSD epoch blocks : 12/15 [00:00<00:00, 21.89it/s]
87%|████████▋ | CSD epoch blocks : 13/15 [00:00<00:00, 21.88it/s]
93%|█████████▎| CSD epoch blocks : 14/15 [00:00<00:00, 21.83it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:00<00:00, 21.83it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:00<00:00, 21.84it/s]
[done]
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/15 [00:00<?, ?it/s]
7%|▋ | CSD epoch blocks : 1/15 [00:00<00:00, 20.38it/s]
13%|█▎ | CSD epoch blocks : 2/15 [00:00<00:00, 20.42it/s]
20%|██ | CSD epoch blocks : 3/15 [00:00<00:00, 20.45it/s]
27%|██▋ | CSD epoch blocks : 4/15 [00:00<00:00, 20.44it/s]
33%|███▎ | CSD epoch blocks : 5/15 [00:00<00:00, 20.47it/s]
40%|████ | CSD epoch blocks : 6/15 [00:00<00:00, 20.48it/s]
47%|████▋ | CSD epoch blocks : 7/15 [00:00<00:00, 20.43it/s]
53%|█████▎ | CSD epoch blocks : 8/15 [00:00<00:00, 20.44it/s]
60%|██████ | CSD epoch blocks : 9/15 [00:00<00:00, 20.45it/s]
67%|██████▋ | CSD epoch blocks : 10/15 [00:00<00:00, 20.46it/s]
73%|███████▎ | CSD epoch blocks : 11/15 [00:00<00:00, 20.46it/s]
80%|████████ | CSD epoch blocks : 12/15 [00:00<00:00, 20.45it/s]
87%|████████▋ | CSD epoch blocks : 13/15 [00:00<00:00, 20.43it/s]
93%|█████████▎| CSD epoch blocks : 14/15 [00:00<00:00, 20.45it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:00<00:00, 20.46it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:00<00:00, 20.45it/s]
[done]
若要计算某频段的源功率而非各频率分别计算，可对 CSD 对象跨频率求平均。
csd = csd.mean()
csd_baseline = csd_baseline.mean()
csd_ers = csd_ers.mean()
使用在整个时间过程上计算的 CSD 计算 DICS 空间滤波器。
fwd = mne.read_forward_solution(fname_fwd)
filters = make_dics(
info,
fwd,
csd,
noise_csd=csd_baseline,
pick_ori="max-power",
reduce_rank=True,
real_filter=True,
)
del fwd
Reading forward solution from /home/circleci/mne_data/MNE-somato-data/derivatives/sub-01/sub-01_task-somato-fwd.fif...
Reading a source space...
[done]
Reading a source space...
[done]
2 source spaces read
Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available
Read MEG forward solution (8155 sources, 306 channels, free orientations)
Source spaces transformed to the forward solution coordinate frame
Identifying common channels ...
Dropped the following channels:
['EOG 061', 'STI 015', 'STI 016', 'STI 006', 'STI 002', 'STI 004', 'STI 014', 'STI 001', 'STI 003', 'STI 005']
Identifying common channels ...
Computing inverse operator with 306 channels.
306 out of 306 channels remain after picking
Selected 306 channels
Creating the depth weighting matrix...
Whitening the forward solution.
Computing rank from covariance with rank=None
Using tolerance 3.1e-14 (2.2e-16 eps * 306 dim * 0.46 max singular value)
Estimated rank (mag + grad): 64
MEG: rank 64 computed from 306 data channels with 0 projectors
Setting small MEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing rank from covariance with rank=None
Using tolerance 3.1e-14 (2.2e-16 eps * 306 dim * 0.46 max singular value)
Estimated rank (mag + grad): 64
MEG: rank 64 computed from 306 data channels with 0 projectors
Computing rank from covariance with rank=None
Using tolerance 3.2e-14 (2.2e-16 eps * 306 dim * 0.47 max singular value)
Estimated rank (mag + grad): 64
MEG: rank 64 computed from 306 data channels with 0 projectors
Computing DICS spatial filters...
Computing beamformer filters for 8155 sources
Filter computation complete
把 DICS 空间滤波器分别应用于由基线和 ERS 活动期间计算的 CSD。
Computing DICS source power...
[done]
Computing DICS source power...
[done]
可视化 ERS 活动期间相对于基线功率的源功率。
stc = beta_source_power / baseline_source_power
message = "DICS source power in the 12-30 Hz frequency band"
brain = stc.plot(
hemi="both",
views="axial",
subjects_dir=subjects_dir,
subject=subject,
time_label=message,
)
Using control points [1.40188307 1.5102619 2.34636424]
参考文献#
脚本总运行时间：（0 分 6.431 秒）