说明
转到末尾下载完整示例代码。
Spectrum 与 EpochsSpectrum 类：频域数据#
本教程介绍如何基于连续 Raw、分段 Epochs 或平均后的 Evoked 数据，创建并可视化数据的频域表示。
照例，我们先导入所需的模块，并加载示例数据集：
import numpy as np
import mne
sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = sample_data_folder / "MEG" / "sample" / "sample_audvis_raw.fif"
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False).crop(tmax=60)
三种传感器空间容器（Raw、Epochs 和 Evoked）都提供 compute_psd() 方法，且可用的选项相同。
Effective window size : 3.410 (s)
默认情况下，连续数据采用 Welch 法 [ 1 ] 进行谱估计，而分段或平均数据采用多锥度（multitaper）法 [ 2 ]。通过向 compute_psd() 方法传入 method='welch' 或 method='multitaper' 可以覆盖这一默认设置。
此外还有许多其他可用选项；例如，我们可以针对给定的时间范围、选定的频率范围以及可用通道的一个子集来计算频谱：
raw.compute_psd(method="multitaper", tmin=10, tmax=20, fmin=5, fmax=30, picks="eeg")
Using multitaper spectrum estimation with 7 DPSS windows
你还可以向底层的谱估计函数传递一些参数，例如 Welch 法的 FFT 窗口长度与重叠长度；详情请参阅 mne.time_frequency.Spectrum 的文档字符串（尤其是其 method_kw 参数）以及谱估计函数 psd_array_welch() 和 psd_array_multitaper()。
对于分段数据，谱估计结果所属的类为 mne.time_frequency.EpochsSpectrum，而不是 mne.time_frequency.Spectrum，但这两个类的大部分 API 相同。例如，两者都有 get_data() 方法，并提供一个可返回各频点频率的选项：
with mne.use_log_level("WARNING"): # hide some irrelevant info messages
events = mne.find_events(raw, stim_channel="STI 014")
event_dict = {
"auditory/left": 1,
"auditory/right": 2,
"visual/left": 3,
"visual/right": 4,
}
epochs = mne.Epochs(
raw, events, tmin=-0.3, tmax=0.7, event_id=event_dict, preload=True
)
epo_spectrum = epochs.compute_psd()
psds, freqs = epo_spectrum.get_data(return_freqs=True)
print(f"\nPSDs shape: {psds.shape}, freqs shape: {freqs.shape}")
epo_spectrum
Using multitaper spectrum estimation with 7 DPSS windows
PSDs shape: (77, 364, 301), freqs shape: (301,)
此外，Spectrum 和 EpochsSpectrum 都实现了 __getitem__ 方法，这意味着可以通过方括号索引访问其数据。对于 Spectrum 对象（由 Raw 或 Evoked 数据计算得到），其索引方式类似于 Raw 对象或 NumPy 数组：
evoked = epochs["auditory"].average()
evk_spectrum = evoked.compute_psd()
# the first 3 frequency bins for the first 4 channels:
print(evk_spectrum[:4, :3])
Using multitaper spectrum estimation with 7 DPSS windows
[[5.61863537e-23 1.13487136e-22 9.88010499e-23]
[2.89065956e-23 4.24056397e-23 3.91160668e-23]
[1.56762015e-25 2.65842303e-25 2.39426775e-25]
[5.38758770e-23 1.03771739e-22 9.98519756e-23]]
相比之下，EpochsSpectrum 的索引方式类似于 Epochs 对象：可以使用字符串值按条件名称选择特定 epoch 的谱估计结果，返回的是一个新建的 EpochsSpectrum 实例，而不是由数据值构成的 NumPy 数组。也可以通过分层事件描述符（hierarchical event descriptors，HEDs）进行选择：
# get both "visual/left" and "visual/right" epochs:
epo_spectrum["visual"]
可视化 Spectrum 对象#
Spectrum 和 EpochsSpectrum 对象都提供绘图方法：plot()（频率 × 功率）、plot_topo()（对每个传感器分别绘制频率 × 功率）以及 plot_topomap()（特定频段内功率的插值头皮地形图）。下面演示其中几个绘图选项；完整细节请参阅文档字符串。
evk_spectrum.plot(picks="data", exclude="bads")
evk_spectrum.plot_topo(color="k", fig_facecolor="w", axis_facecolor="w")
Plotting power spectral density (dB=True).
evk_spectrum.plot_topomap(ch_type="eeg", agg_fun=np.median)
迁移旧版代码#
下表快速对照了引入 Spectrum 和 EpochsSpectrum 类前后实现相同功能的代码。
| 旧写法 | 新写法 |
|---|---|
|
|
|
|
|
|
|
|
|
|
警告
函数 mne.time_frequency.psd_welch 和 mne.time_frequency.psd_multitaper 已被移除；新代码应使用 Raw.compute_psd()、Epochs.compute_psd() 和 Evoked.compute_psd() 方法，并将 method='welch' 或 method='multitaper' 作为参数传入。
类方法 Raw.plot_psd()、Epochs.plot_psd()、Raw.plot_psd_topo() 和 Epochs.plot_psd_topomap() 仍保留在 API 中，以支持旧版代码，但在编写新代码时应避免使用。
参考文献#
脚本总运行时间：（0 分钟 5.434 秒）