非官方中文翻译 本页译自 MNE-Python 官方文档： 查看英文原文 ↗。 文档版权归 MNE Developers 所有（BSD-3-Clause）。

说明

转到末尾下载完整示例代码。

计算互谱密度（CSD）矩阵#

互谱密度（CSD）矩阵类似协方差矩阵，只是位于时频域。它是计算传感器间相干性或 DICS 波束形成器的第一步。

本脚本演示 MNE-Python 提供的三种 CSD 计算方法：

  1. 使用短时傅里叶变换：mne.time_frequency.csd_fourier()

  2. 使用多锥度方法：mne.time_frequency.csd_multitaper()

  3. 使用 Morlet 小波：mne.time_frequency.csd_morlet()

# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import mne
from mne.datasets import sample
from mne.time_frequency import csd_fourier, csd_morlet, csd_multitaper

print(__doc__)

在下面的示例中，CSD 矩阵的计算可以使用多核。把 n_jobs 设为大于 1 的值即可指定使用的核数。

n_jobs = 1

加载示例数据集。

data_path = sample.data_path()
meg_path = data_path / "MEG" / "sample"
fname_raw = meg_path / "sample_audvis_raw.fif"
fname_event = meg_path / "sample_audvis_raw-eve.fif"
raw = mne.io.read_raw_fif(fname_raw)
events = mne.read_events(fname_event)

Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.

默认情况下，CSD 矩阵使用全部 MEG/EEG 通道计算。解释混合传感器类型的 CSD 矩阵时要注意，各传感器的测量单位不同，因而缩放也不同。本示例为速度和清晰起见只选择一种通道类型：梯度计。

picks = mne.pick_types(raw.info, meg="grad")

# Make some epochs, based on events with trigger code 1
epochs = mne.Epochs(
    raw,
    events,
    event_id=1,
    tmin=-0.2,
    tmax=1,
    picks=picks,
    baseline=(None, 0),
    reject=dict(grad=4000e-13),
    preload=True,
)

Not setting metadata
72 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] s
Applying baseline correction (mode: mean)
0 projection items activated
Loading data for 72 events and 722 original time points ...
0 bad epochs dropped

使用短时傅里叶变换与（自适应）多锥度计算 CSD 矩阵非常直接：

csd_fft = csd_fourier(epochs, fmin=15, fmax=20, n_jobs=n_jobs)
csd_mt = csd_multitaper(epochs, fmin=15, fmax=20, adaptive=True, n_jobs=n_jobs)

Computing cross-spectral density from epochs...

  0%|          | CSD epoch blocks : 0/72 [00:00<?,       ?it/s]
  7%|▋         | CSD epoch blocks : 5/72 [00:00<00:00,  299.40it/s]
 14%|█▍        | CSD epoch blocks : 10/72 [00:00<00:00,  303.77it/s]
 21%|██        | CSD epoch blocks : 15/72 [00:00<00:00,  306.43it/s]
 28%|██▊       | CSD epoch blocks : 20/72 [00:00<00:00,  307.53it/s]
 36%|███▌      | CSD epoch blocks : 26/72 [00:00<00:00,  309.05it/s]
 43%|████▎     | CSD epoch blocks : 31/72 [00:00<00:00,  308.53it/s]
 50%|█████     | CSD epoch blocks : 36/72 [00:00<00:00,  309.14it/s]
 58%|█████▊    | CSD epoch blocks : 42/72 [00:00<00:00,  309.85it/s]
 67%|██████▋   | CSD epoch blocks : 48/72 [00:00<00:00,  310.35it/s]
 74%|███████▎  | CSD epoch blocks : 53/72 [00:00<00:00,  310.45it/s]
 81%|████████  | CSD epoch blocks : 58/72 [00:00<00:00,  309.64it/s]
 88%|████████▊ | CSD epoch blocks : 63/72 [00:00<00:00,  308.21it/s]
 94%|█████████▍| CSD epoch blocks : 68/72 [00:00<00:00,  307.11it/s]
100%|██████████| CSD epoch blocks : 72/72 [00:00<00:00,  306.13it/s]
[done]
    Using multitaper spectrum estimation with 7 DPSS windows
Computing cross-spectral density from epochs...

  0%|          | CSD epoch blocks : 0/72 [00:00<?,       ?it/s]
  1%|▏         | CSD epoch blocks : 1/72 [00:00<00:02,   28.57it/s]
  3%|▎         | CSD epoch blocks : 2/72 [00:00<00:02,   29.51it/s]
  4%|▍         | CSD epoch blocks : 3/72 [00:00<00:02,   29.70it/s]
  6%|▌         | CSD epoch blocks : 4/72 [00:00<00:02,   29.72it/s]
  7%|▋         | CSD epoch blocks : 5/72 [00:00<00:02,   29.46it/s]
  8%|▊         | CSD epoch blocks : 6/72 [00:00<00:02,   29.58it/s]
 10%|▉         | CSD epoch blocks : 7/72 [00:00<00:02,   29.50it/s]
 11%|█         | CSD epoch blocks : 8/72 [00:00<00:02,   29.55it/s]
 12%|█▎        | CSD epoch blocks : 9/72 [00:00<00:02,   29.57it/s]
 14%|█▍        | CSD epoch blocks : 10/72 [00:00<00:02,   29.50it/s]
 15%|█▌        | CSD epoch blocks : 11/72 [00:00<00:02,   29.36it/s]
 17%|█▋        | CSD epoch blocks : 12/72 [00:00<00:02,   29.34it/s]
 18%|█▊        | CSD epoch blocks : 13/72 [00:00<00:02,   29.35it/s]
 19%|█▉        | CSD epoch blocks : 14/72 [00:00<00:01,   29.46it/s]
 21%|██        | CSD epoch blocks : 15/72 [00:00<00:01,   29.52it/s]
 22%|██▏       | CSD epoch blocks : 16/72 [00:00<00:01,   29.50it/s]
 24%|██▎       | CSD epoch blocks : 17/72 [00:00<00:01,   29.56it/s]
 25%|██▌       | CSD epoch blocks : 18/72 [00:00<00:01,   29.63it/s]
 26%|██▋       | CSD epoch blocks : 19/72 [00:00<00:01,   29.61it/s]
 28%|██▊       | CSD epoch blocks : 20/72 [00:00<00:01,   29.60it/s]
 29%|██▉       | CSD epoch blocks : 21/72 [00:00<00:01,   29.58it/s]
 31%|███       | CSD epoch blocks : 22/72 [00:00<00:01,   29.57it/s]
 32%|███▏      | CSD epoch blocks : 23/72 [00:00<00:01,   29.56it/s]
 33%|███▎      | CSD epoch blocks : 24/72 [00:00<00:01,   29.48it/s]
 35%|███▍      | CSD epoch blocks : 25/72 [00:00<00:01,   29.44it/s]
 36%|███▌      | CSD epoch blocks : 26/72 [00:00<00:01,   29.40it/s]
 38%|███▊      | CSD epoch blocks : 27/72 [00:00<00:01,   29.36it/s]
 39%|███▉      | CSD epoch blocks : 28/72 [00:00<00:01,   29.28it/s]
 40%|████      | CSD epoch blocks : 29/72 [00:00<00:01,   29.30it/s]
 42%|████▏     | CSD epoch blocks : 30/72 [00:01<00:01,   29.31it/s]
 43%|████▎     | CSD epoch blocks : 31/72 [00:01<00:01,   29.26it/s]
 44%|████▍     | CSD epoch blocks : 32/72 [00:01<00:01,   29.28it/s]
 46%|████▌     | CSD epoch blocks : 33/72 [00:01<00:01,   29.30it/s]
 47%|████▋     | CSD epoch blocks : 34/72 [00:01<00:01,   29.32it/s]
 49%|████▊     | CSD epoch blocks : 35/72 [00:01<00:01,   29.33it/s]
 50%|█████     | CSD epoch blocks : 36/72 [00:01<00:01,   29.34it/s]
 51%|█████▏    | CSD epoch blocks : 37/72 [00:01<00:01,   29.29it/s]
 53%|█████▎    | CSD epoch blocks : 38/72 [00:01<00:01,   29.25it/s]
 54%|█████▍    | CSD epoch blocks : 39/72 [00:01<00:01,   28.95it/s]
 56%|█████▌    | CSD epoch blocks : 40/72 [00:01<00:01,   28.94it/s]
 57%|█████▋    | CSD epoch blocks : 41/72 [00:01<00:01,   28.89it/s]
 58%|█████▊    | CSD epoch blocks : 42/72 [00:01<00:01,   28.91it/s]
 60%|█████▉    | CSD epoch blocks : 43/72 [00:01<00:01,   28.92it/s]
 61%|██████    | CSD epoch blocks : 44/72 [00:01<00:00,   28.97it/s]
 62%|██████▎   | CSD epoch blocks : 45/72 [00:01<00:00,   29.07it/s]
 64%|██████▍   | CSD epoch blocks : 46/72 [00:01<00:00,   28.92it/s]
 65%|██████▌   | CSD epoch blocks : 47/72 [00:01<00:00,   29.00it/s]
 67%|██████▋   | CSD epoch blocks : 48/72 [00:01<00:00,   29.04it/s]
 68%|██████▊   | CSD epoch blocks : 49/72 [00:01<00:00,   29.08it/s]
 69%|██████▉   | CSD epoch blocks : 50/72 [00:01<00:00,   28.40it/s]
 71%|███████   | CSD epoch blocks : 51/72 [00:01<00:00,   28.51it/s]
 72%|███████▏  | CSD epoch blocks : 52/72 [00:01<00:00,   28.59it/s]
 74%|███████▎  | CSD epoch blocks : 53/72 [00:01<00:00,   28.69it/s]
 75%|███████▌  | CSD epoch blocks : 54/72 [00:01<00:00,   28.74it/s]
 76%|███████▋  | CSD epoch blocks : 55/72 [00:01<00:00,   28.82it/s]
 78%|███████▊  | CSD epoch blocks : 56/72 [00:01<00:00,   28.89it/s]
 79%|███████▉  | CSD epoch blocks : 57/72 [00:01<00:00,   28.96it/s]
 81%|████████  | CSD epoch blocks : 58/72 [00:01<00:00,   29.03it/s]
 82%|████████▏ | CSD epoch blocks : 59/72 [00:02<00:00,   29.09it/s]
 83%|████████▎ | CSD epoch blocks : 60/72 [00:02<00:00,   29.15it/s]
 85%|████████▍ | CSD epoch blocks : 61/72 [00:02<00:00,   29.20it/s]
 86%|████████▌ | CSD epoch blocks : 62/72 [00:02<00:00,   29.25it/s]
 88%|████████▊ | CSD epoch blocks : 63/72 [00:02<00:00,   29.29it/s]
 89%|████████▉ | CSD epoch blocks : 64/72 [00:02<00:00,   29.35it/s]
 90%|█████████ | CSD epoch blocks : 65/72 [00:02<00:00,   29.39it/s]
 92%|█████████▏| CSD epoch blocks : 66/72 [00:02<00:00,   29.41it/s]
 93%|█████████▎| CSD epoch blocks : 67/72 [00:02<00:00,   29.45it/s]
 94%|█████████▍| CSD epoch blocks : 68/72 [00:02<00:00,   29.50it/s]
 96%|█████████▌| CSD epoch blocks : 69/72 [00:02<00:00,   29.53it/s]
 97%|█████████▋| CSD epoch blocks : 70/72 [00:02<00:00,   29.56it/s]
 99%|█████████▊| CSD epoch blocks : 71/72 [00:02<00:00,   29.58it/s]
100%|██████████| CSD epoch blocks : 72/72 [00:02<00:00,   29.61it/s]
100%|██████████| CSD epoch blocks : 72/72 [00:02<00:00,   29.31it/s]
[done]

使用 Morlet 小波计算 CSD 时，需要指定确切的计算频率。对每个频率会构造对应的小波并与信号卷积，得到时频分解。

CSD 通过计算所有传感器对之间的时频表示相关性来构建。时频分解原本与信号采样率相同，本例约 600 Hz，这意味着时间上过采样；计算 CSD 时不必使用全部采样点，只要足以得到可靠的相关统计即可。设置 decim=10 表示每 10 个采样点取一个，可显著加快计算，对 CSD 影响很小。

frequencies = [16, 17, 18, 19, 20]
csd_wav = csd_morlet(epochs, frequencies, decim=10, n_jobs=n_jobs)

Computing cross-spectral density from epochs...

  0%|          | CSD epoch blocks : 0/72 [00:00<?,       ?it/s]
  1%|▏         | CSD epoch blocks : 1/72 [00:00<00:01,   42.09it/s]
  3%|▎         | CSD epoch blocks : 2/72 [00:00<00:01,   41.82it/s]
  4%|▍         | CSD epoch blocks : 3/72 [00:00<00:01,   41.94it/s]
  6%|▌         | CSD epoch blocks : 4/72 [00:00<00:01,   42.08it/s]
  7%|▋         | CSD epoch blocks : 5/72 [00:00<00:01,   42.16it/s]
  8%|▊         | CSD epoch blocks : 6/72 [00:00<00:01,   42.22it/s]
 10%|▉         | CSD epoch blocks : 7/72 [00:00<00:01,   42.28it/s]
 11%|█         | CSD epoch blocks : 8/72 [00:00<00:01,   42.26it/s]
 12%|█▎        | CSD epoch blocks : 9/72 [00:00<00:01,   42.23it/s]
 14%|█▍        | CSD epoch blocks : 10/72 [00:00<00:01,   42.29it/s]
 15%|█▌        | CSD epoch blocks : 11/72 [00:00<00:01,   42.33it/s]
 17%|█▋        | CSD epoch blocks : 12/72 [00:00<00:01,   42.29it/s]
 18%|█▊        | CSD epoch blocks : 13/72 [00:00<00:01,   42.30it/s]
 19%|█▉        | CSD epoch blocks : 14/72 [00:00<00:01,   42.31it/s]
 21%|██        | CSD epoch blocks : 15/72 [00:00<00:01,   42.34it/s]
 22%|██▏       | CSD epoch blocks : 16/72 [00:00<00:01,   42.32it/s]
 24%|██▎       | CSD epoch blocks : 17/72 [00:00<00:01,   42.34it/s]
 25%|██▌       | CSD epoch blocks : 18/72 [00:00<00:01,   42.34it/s]
 26%|██▋       | CSD epoch blocks : 19/72 [00:00<00:01,   42.37it/s]
 28%|██▊       | CSD epoch blocks : 20/72 [00:00<00:01,   42.35it/s]
 29%|██▉       | CSD epoch blocks : 21/72 [00:00<00:01,   42.35it/s]
 31%|███       | CSD epoch blocks : 22/72 [00:00<00:01,   42.37it/s]
 32%|███▏      | CSD epoch blocks : 23/72 [00:00<00:01,   42.38it/s]
 33%|███▎      | CSD epoch blocks : 24/72 [00:00<00:01,   42.40it/s]
 35%|███▍      | CSD epoch blocks : 25/72 [00:00<00:01,   42.38it/s]
 36%|███▌      | CSD epoch blocks : 26/72 [00:00<00:01,   42.39it/s]
 38%|███▊      | CSD epoch blocks : 27/72 [00:00<00:01,   42.40it/s]
 39%|███▉      | CSD epoch blocks : 28/72 [00:00<00:01,   42.40it/s]
 40%|████      | CSD epoch blocks : 29/72 [00:00<00:01,   42.36it/s]
 42%|████▏     | CSD epoch blocks : 30/72 [00:00<00:00,   42.37it/s]
 43%|████▎     | CSD epoch blocks : 31/72 [00:00<00:00,   42.39it/s]
 44%|████▍     | CSD epoch blocks : 32/72 [00:00<00:00,   42.40it/s]
 46%|████▌     | CSD epoch blocks : 33/72 [00:00<00:00,   42.40it/s]
 47%|████▋     | CSD epoch blocks : 34/72 [00:00<00:00,   42.42it/s]
 49%|████▊     | CSD epoch blocks : 35/72 [00:00<00:00,   42.43it/s]
 50%|█████     | CSD epoch blocks : 36/72 [00:00<00:00,   42.45it/s]
 51%|█████▏    | CSD epoch blocks : 37/72 [00:00<00:00,   42.44it/s]
 53%|█████▎    | CSD epoch blocks : 38/72 [00:00<00:00,   42.43it/s]
 54%|█████▍    | CSD epoch blocks : 39/72 [00:00<00:00,   42.43it/s]
 56%|█████▌    | CSD epoch blocks : 40/72 [00:00<00:00,   42.43it/s]
 57%|█████▋    | CSD epoch blocks : 41/72 [00:00<00:00,   42.44it/s]
 58%|█████▊    | CSD epoch blocks : 42/72 [00:00<00:00,   42.42it/s]
 60%|█████▉    | CSD epoch blocks : 43/72 [00:01<00:00,   42.42it/s]
 61%|██████    | CSD epoch blocks : 44/72 [00:01<00:00,   42.39it/s]
 62%|██████▎   | CSD epoch blocks : 45/72 [00:01<00:00,   42.39it/s]
 64%|██████▍   | CSD epoch blocks : 46/72 [00:01<00:00,   42.37it/s]
 65%|██████▌   | CSD epoch blocks : 47/72 [00:01<00:00,   42.37it/s]
 67%|██████▋   | CSD epoch blocks : 48/72 [00:01<00:00,   42.38it/s]
 68%|██████▊   | CSD epoch blocks : 49/72 [00:01<00:00,   42.39it/s]
 69%|██████▉   | CSD epoch blocks : 50/72 [00:01<00:00,   42.37it/s]
 71%|███████   | CSD epoch blocks : 51/72 [00:01<00:00,   42.40it/s]
 72%|███████▏  | CSD epoch blocks : 52/72 [00:01<00:00,   42.41it/s]
 74%|███████▎  | CSD epoch blocks : 53/72 [00:01<00:00,   42.42it/s]
 75%|███████▌  | CSD epoch blocks : 54/72 [00:01<00:00,   42.42it/s]
 76%|███████▋  | CSD epoch blocks : 55/72 [00:01<00:00,   42.43it/s]
 78%|███████▊  | CSD epoch blocks : 56/72 [00:01<00:00,   42.44it/s]
 79%|███████▉  | CSD epoch blocks : 57/72 [00:01<00:00,   42.45it/s]
 81%|████████  | CSD epoch blocks : 58/72 [00:01<00:00,   42.46it/s]
 82%|████████▏ | CSD epoch blocks : 59/72 [00:01<00:00,   42.45it/s]
 83%|████████▎ | CSD epoch blocks : 60/72 [00:01<00:00,   42.46it/s]
 85%|████████▍ | CSD epoch blocks : 61/72 [00:01<00:00,   42.48it/s]
 86%|████████▌ | CSD epoch blocks : 62/72 [00:01<00:00,   42.49it/s]
 88%|████████▊ | CSD epoch blocks : 63/72 [00:01<00:00,   42.47it/s]
 89%|████████▉ | CSD epoch blocks : 64/72 [00:01<00:00,   42.47it/s]
 90%|█████████ | CSD epoch blocks : 65/72 [00:01<00:00,   42.47it/s]
 92%|█████████▏| CSD epoch blocks : 66/72 [00:01<00:00,   42.48it/s]
 93%|█████████▎| CSD epoch blocks : 67/72 [00:01<00:00,   42.47it/s]
 94%|█████████▍| CSD epoch blocks : 68/72 [00:01<00:00,   42.46it/s]
 96%|█████████▌| CSD epoch blocks : 69/72 [00:01<00:00,   42.47it/s]
 97%|█████████▋| CSD epoch blocks : 70/72 [00:01<00:00,   42.46it/s]
 99%|█████████▊| CSD epoch blocks : 71/72 [00:01<00:00,   42.45it/s]
100%|██████████| CSD epoch blocks : 72/72 [00:01<00:00,   42.45it/s]
100%|██████████| CSD epoch blocks : 72/72 [00:01<00:00,   42.41it/s]
[done]

得到的 mne.time_frequency.CrossSpectralDensity 对象带有绘图函数，可用于比较不同方法的结果。这里绘制跨频率的平均 CSD。mne.time_frequency.CrossSpectralDensity.plot() 返回所创建图形的列表；本例每个列表只含一个图形，因此用 Python 技巧在变量名后加逗号，把图形（而非列表）赋给 fig：

plot_dict = {
    "Short-time Fourier transform": csd_fft,
    "Adaptive multitapers": csd_mt,
    "Morlet wavelet transform": csd_wav,
}
for title, csd in plot_dict.items():
    (fig,) = csd.mean().plot()
    fig.suptitle(title)
  • Short-time Fourier transform, 15.8-20.0 Hz.
  • Adaptive multitapers, 15.8-20.0 Hz.
  • Morlet wavelet transform, 16.0-20.0 Hz.

脚本总运行时间：（0 分 4.964 秒）

图库由 Sphinx-Gallery 生成