【神经数据分析Python包介绍】MNE-Python及其扩展包(上)

MNE-Python 是 Python 生态中最成熟的神经电生理信号处理框架。本教程专注于 LFP（Local Field Potential）分析，所有代码示例均以多通道 LFP 数据为场景。

你不需要用 MNE 做所有事。MNE 的核心价值在于：

1. 统一的数据结构：Raw → Epochs → Evoked 流水线
2. 可靠的滤波：经过充分验证的 FIR/IIR 实现
3. 统计检验：排列簇检验等非参数方法
4. 元数据追踪：自动记录所有处理步骤

对于 FOOOF、PAC、Granger 因果等高级分析，用 MNE 做预处理，然后导出 numpy 数组给专门的包（如 fooof、tensorpac、spectral_connectivity）处理，是完全合理的工作流。

0. 安装

pip install mne

验证：

import mneprint(mne.__version__)mne.sys_info()

LFP 分析不需要 mne[full]

mne[full] 安装的 pyvista / vtk / nibabel 是为 EEG/MEG 源定位的 3D 可视化准备的。LFP 分析用不到，装 mne 就够了。

1. 核心数据结构

MNE 有三个你需要关心的数据结构：

Raw（连续数据）──切割──▶ Epochs（分段数据）──平均──▶ Evoked（平均响应）

LFP 分析中，Source Estimate 不会用到（那是 EEG/MEG 从头皮反推脑内源的工具）。

1.1 Info：元数据容器

Info 是贯穿 MNE 所有数据结构的元数据字典。它记录采样率、通道名、通道类型等信息。

MNE 的滤波、PSD、时频等函数都从 Info 里读取采样率和通道信息，不需要你每次手动传参。

import mneimport numpy as np# === 为 LFP 数据手动创建 Info ===# 假设你有 32 通道 LFP，采样率 1000 Hzn_channels = 32sfreq = 1000.0# 通道名：按脑区命名（这是 LFP 的典型做法）ch_names = [f'mPFC-{i}'for i in range(1, 9)] + \           [f'HPC-{i}'for i in range(1, 9)] + \           [f'AMG-{i}'for i in range(1, 9)] + \           [f'VTA-{i}'for i in range(1, 9)]# 通道类型：全部设为 'seeg'（立体脑电/深部电极）# MNE 中 LFP 通道推荐用 'seeg'（stereo-EEG）或 'ecog'（皮层电极）ch_types = ['seeg'] * n_channelsinfo = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types=ch_types)print(info)

MNE 没有专门的 'lfp' 类型。根据你的电极：

• 深部探针（硅探针、四极管、微丝阵列）→ 'seeg'
• 皮层表面电极（ECoG grid）→ 'ecog'
• DBS 导联 → 'dbs'
• 都不是 / 不确定 → 'misc'（不影响分析，只影响默认绘图设置）

**不建议把 LFP 标成 'eeg'**。MNE 不会因为通道类型是 'eeg' 就自动套用 average reference，但 'eeg' 会触发头皮 EEG 的默认缩放、绘图和部分工作流分支；LFP 一般应标成 'seeg'、'ecog'、'dbs' 或 'misc'。

备忘一下字段：

1.2 Raw：连续数据

Raw 是 MNE 中存储连续时间序列的容器。对 LFP 分析来说，你最常做的事是从 numpy 数组创建 Raw 对象。

# === 从 numpy 数组创建 Raw ===# 假设 lfp_data 是你从录制系统导出的数据# shape: (n_channels, n_samples) = (32, 300000)  即 32 通道，300 秒lfp_data = np.random.randn(32, 300000) * 1e-4# 模拟数据，单位：伏特# 注意：MNE 内部默认单位是伏特 (V)# 如果你的数据是 µV，需要先乘以 1e-6# 如果是 mV，乘以 1e-3raw = mne.io.RawArray(lfp_data, info)print(raw)# <RawArray | 32 x 300000 (300.0 s), ~73.2 MB, data loaded>

这里注意一下默认单位问题。

MNE 内部一律用 SI 单位：

• 电压：**伏特 (V)**，不是 µV
• 如果你的 LFP 数据是 µV（大多数系统导出的默认单位），创建 Raw 前必须 data *= 1e-6
• 否则 MNE 的滤波阈值、拒绝标准、绘图 y 轴全部会出错

Raw 对象的核心操作：

# === 时间切片 ===raw_seg = raw.copy().crop(tmin=10.0, tmax=60.0)  # 截取 10-60 秒# ⚠ crop 是 in-place 操作，不 copy() 会修改原始数据！# === 通道选择 ===raw_hpc = raw.copy().pick([f'HPC-{i}'for i in range(1, 9)])  # 只要海马通道# === 提取 numpy 数组 ===data, times = raw.get_data(return_times=True)# data.shape = (32, 300000)# times.shape = (300000,)  单位：秒# 提取特定通道 + 时间段data_hpc = raw.get_data(picks=['HPC-1', 'HPC-2'], tmin=10.0, tmax=20.0)# === 标记坏道 ===raw.info['bads'] = ['mPFC-3', 'AMG-7']  # 手动标记# 坏道不会自动从对象里删除，但很多函数默认会把 info['bads'] 排除在 picks 之外# 如果你显式传入这些通道名，它们仍然会出现在结果里# === 添加时间标注 ===annot = mne.Annotations(    onset=[15.0, 45.0, 120.0],           # 开始时间 (秒)    duration=[2.0, 3.0, 1.0],            # 持续时间 (秒)    description=['bad_noise', 'bad_movement', 'stim_on'])raw.set_annotations(annot)# 'bad_*' 开头的标注在 epoching 时自动被丢弃

1.3 Events 与 Epochs：分段数据

LFP 实验通常有外部事件（刺激、行为标记），你需要围绕这些事件切分数据。

事件 (Events) 是一个 (n_events, 3) 的整数数组：

[sample_index, 0, event_id]

# === 方法 1：从 numpy 数组手动创建事件 ===# 假设你有刺激时间戳（秒），先转换为采样点stim_times = np.array([5.0, 15.0, 25.0, 35.0, 45.0])  # 5 次刺激stim_samples = (stim_times * sfreq).astype(int)# 构造 events 数组events = np.column_stack([    stim_samples,    np.zeros(len(stim_samples), dtype=int),    np.ones(len(stim_samples), dtype=int)   # event_id = 1])print(events)# [[ 5000     0     1]#  [15000     0     1]#  [25000     0     1]#  ...]# === 方法 2：从 Annotations 创建事件 ===events_from_annot, event_id_map = mne.events_from_annotations(raw)# === 方法 3：多种事件 ===# 假设有两类刺激：光刺激 (id=1) 和声音刺激 (id=2)light_times = np.array([5.0, 25.0, 45.0])sound_times = np.array([15.0, 35.0, 55.0])events_light = np.column_stack([    (light_times * sfreq).astype(int),    np.zeros(len(light_times), dtype=int),    np.ones(len(light_times), dtype=int)])events_sound = np.column_stack([    (sound_times * sfreq).astype(int),    np.zeros(len(sound_times), dtype=int),    np.full(len(sound_times), 2, dtype=int)])events_all = np.vstack([events_light, events_sound])events_all = events_all[events_all[:, 0].argsort()]  # 按时间排序event_dict = {'light_stim': 1, 'sound_stim': 2}

创建 Epochs：

# LFP 的 reject 标准和 EEG 完全不同# seeg 通道的单位是 V，典型 LFP 振幅 ~100-500 µVreject_criteria = dict(seeg=500e-6)  # 拒绝超过 500 µV 的 epochepochs = mne.Epochs(    raw,    events_all,    event_id=event_dict,    tmin=-0.5,                # epoch 起始：事件前 500 ms    tmax=1.0,                 # epoch 结束：事件后 1000 ms    baseline=(-0.5, 0),       # 基线校正：用事件前 500ms 的均值    reject=reject_criteria,    preload=True)print(epochs)# === Epochs 索引 ===epochs_light = epochs['light_stim']       # 只取光刺激epochs_sound = epochs['sound_stim']       # 只取声音刺激# === 提取数据 ===data_3d = epochs.get_data()  # shape: (n_epochs, n_channels, n_times)# === 丢弃特定 epoch ===epochs.drop([0, 3])           # 手动丢弃第 0 和第 3 个epochs.drop_bad()             # 按 reject 标准自动丢弃print(epochs.drop_log)        # 查看每个 epoch 的丢弃原因

1.4 Evoked：平均响应

evoked_light = epochs['light_stim'].average()evoked_sound = epochs['sound_stim'].average()print(evoked_light)# <Evoked | 'light_stim' (average, N=3), ...># 提取数据evoked_data = evoked_light.get_data()  # shape: (n_channels, n_times)evoked_times = evoked_light.times       # 时间轴 (秒)# 计算差异波evoked_diff = mne.combine_evoked([evoked_light, evoked_sound], weights=[1, -1])

2. 数据 I/O

2.1 从 numpy 数组创建（最常用）

大多数 LFP 分析场景，你的数据已经被上游工具（Spike2、NeuroExplorer、自定义 MATLAB 脚本）读取并导出为 .mat、.npy 或 .csv。这时直接从 numpy 创建 Raw 是最可靠的：

import numpy as npimport scipy.io as sioimport mne# === 从 .mat 文件 ===mat = sio.loadmat('lfp_recording.mat')lfp_data = mat['lfp']          # shape: (n_channels, n_samples)sfreq = float(mat['fs'][0, 0]) # 采样率# 确保单位是 Vif lfp_data.max() > 0.01:      # 数据可能是 µV    lfp_data = lfp_data * 1e-6# 创建 Info 和 Rawch_names = [f'ch{i+1}'for i in range(lfp_data.shape[0])]info = mne.create_info(ch_names, sfreq, ch_types='seeg')raw = mne.io.RawArray(lfp_data, info)# === 从 .npy 文件 ===lfp_data = np.load('lfp_32ch.npy')  # (32, 600000)raw = mne.io.RawArray(lfp_data * 1e-6, info)# === 从 .csv（少量数据时可用）===import pandas as pddf = pd.read_csv('lfp_export.csv')lfp_data = df.iloc[:, 1:].values.T  # 转置为 (n_channels, n_samples)

2.2 读取录制系统原始格式

MNE 核心包内置了对常见神经电生理格式的读取支持：

# === Neuralynx ===raw = mne.io.read_raw_neuralynx('path/to/neuralynx_folder/', preload=True)# 会自动读取文件夹中所有 .ncs 文件作为各通道# === EDF ===raw = mne.io.read_raw_edf('recording.edf', preload=True)# EDF 是很多系统的通用导出格式# === FIF（MNE 原生格式）===raw = mne.io.read_raw_fif('processed_lfp.fif', preload=True)

不管你的录制系统是什么，最可靠的路线是：

1. 用系统自带工具或 Neo（pip install neo）把数据读进 Python
2. 提取为 numpy 数组 (n_channels, n_samples) + 采样率
3. 用 mne.io.RawArray() 创建 Raw 对象

这比依赖 MNE 的格式读取器更不容易出错（尤其是 Blackrock、Plexon 等格式）。

2.3 保存与导出

# === 保存为 MNE 原生 .fif 格式 ===raw.save('lfp_preprocessed-raw.fif', overwrite=True)# MNE 的 .fif 格式保留所有元数据（坏道、标注、滤波历史等）# 保存 Epochsepochs.save('lfp_epochs-epo.fif', overwrite=True)# 保存 Evokedmne.write_evokeds('lfp_evoked-ave.fif', [evoked_light, evoked_sound])# === 导出为 numpy ===np.save('lfp_filtered.npy', raw.get_data())# === 导出为 EDF ===raw.export('lfp_export.edf', overwrite=True)

3. 通道管理

LFP 分析中，通道管理的核心任务是按脑区分组（而非 EEG 那样按头皮位置排列 montage）。

3.1 通道选择与重命名

# === 按名字选择 ===raw_hpc = raw.copy().pick([f'HPC-{i}'for i in range(1, 9)])# === 按类型选择 ===raw_seeg = raw.copy().pick('seeg')  # 只保留 seeg 类型通道# === 重命名通道 ===# 如果原始通道名是 ch1, ch2, ... 需要改成有意义的名字rename_map = {'ch1': 'mPFC-superficial','ch2': 'mPFC-deep','ch3': 'dHPC-CA1','ch4': 'dHPC-CA3',}raw.rename_channels(rename_map)# === 修改通道类型 ===# 把某些通道从 seeg 改为其他类型raw.set_channel_types({'EMG-1': 'emg',     # 肌电'STIM': 'stim',     # 刺激通道'RESP': 'bio',      # 呼吸信号})

3.2 按脑区分组操作

# 定义脑区分组regions = {'mPFC': [f'mPFC-{i}'for i in range(1, 9)],'HPC':  [f'HPC-{i}'for i in range(1, 9)],'AMG':  [f'AMG-{i}'for i in range(1, 9)],'VTA':  [f'VTA-{i}'for i in range(1, 9)],}# 按脑区提取数据for region, chs in regions.items():    data_region = raw.get_data(picks=chs)    print(f'{region}: {data_region.shape}')# mPFC: (8, 300000)# HPC:  (8, 300000)# ...# 对 Epochs 同样适用for region, chs in regions.items():    data_epoch_region = epochs.get_data(picks=chs)    print(f'{region}: {data_epoch_region.shape}')# mPFC: (n_epochs, 8, n_times)

3.3 通道插值（坏道修复）

EEG 中常用球面样条插值修复坏道，但 LFP 不太适用（通道间距太近，相邻通道信号可能完全不同）。**直接排除坏道通常是更好的选择。

# 标记坏道raw.info['bads'] = ['HPC-3']# 插值坏道（用周围通道加权平均）# 对 LFP 来说，这需要通道有位置信息# 如果没有电极位置，可以手动设置# 或者直接排除坏道而不插值raw_interp = raw.copy().interpolate_bads()# 更简单的做法：直接排除坏道raw_clean = raw.copy().pick(    [ch for ch in raw.ch_names if ch notin raw.info['bads']])

4. 预处理：滤波

滤波是 LFP 预处理的核心。MNE 内置了 FIR 和 IIR 滤波器。

4.1 带通滤波

# === 基本带通：1-300 Hz（去除直流漂移 + 高频噪声）===raw_filt = raw.copy().filter(l_freq=1.0, h_freq=300.0)# === 提取特定频段 ===raw_theta = raw.copy().filter(l_freq=4.0, h_freq=12.0)    # theta (4-12 Hz)raw_gamma = raw.copy().filter(l_freq=30.0, h_freq=100.0)   # gamma (30-100 Hz)raw_ripple = raw.copy().filter(l_freq=100.0, h_freq=250.0) # ripple (100-250 Hz)

4.2 滤波器参数详解

# MNE 默认使用 FIR 滤波器 (method='fir')raw_filt = raw.copy().filter(    l_freq=1.0,          # 高通截止频率 (Hz)，None = 不做高通    h_freq=300.0,        # 低通截止频率 (Hz)，None = 不做低通    method='fir',        # 'fir' 或 'iir'    fir_design='firwin', # FIR 设计方法    phase='zero',        # 'zero'=零相位(filtfilt), 'minimum'=因果滤波    picks='seeg',        # 只滤波 seeg 通道    n_jobs=-1,           # 并行处理)# === 使用 IIR 滤波器 ===# 当需要更陡的过渡带时（比如去除 50/60 Hz 工频）iir_params = dict(order=4, ftype='butter')raw_filt_iir = raw.copy().filter(    l_freq=1.0, h_freq=300.0,    method='iir',    iir_params=iir_params)

FIR vs IIR 怎么选？

	FIR（MNE 默认）	IIR
相位	零相位（不引入延迟）	需要 filtfilt 才能零相位
过渡带	较宽（需要更多系数）	较窄（低阶即可）
稳定性	永远稳定	可能不稳定
推荐场景	离线分析（默认选 FIR 就行）	需要极窄过渡带

对 LFP 离线分析，用 MNE 默认的 FIR 就够了。

4.3 陷波滤波（去除工频干扰）

# 去除 50 Hz 及其谐波（中国电网频率）raw_notch = raw.copy().notch_filter(    freqs=[50, 100, 150, 200, 250],  # 50 Hz + 谐波    method='fir',                      # 或 'spectrum_fit'    picks='seeg')# 'spectrum_fit' 方法更精准（先拟合频谱再减去）raw_notch2 = raw.copy().notch_filter(    freqs=[50, 100, 150, 200, 250],    method='spectrum_fit',    picks='seeg')

50 Hz 还是 60 Hz？

• 中国大陆、欧洲、日本（东）：50 Hz
• 美国、日本（西）：60 Hz
• 一定要根据实验所在地选择正确的工频频率

4.4 查看滤波器特性

# 可视化滤波器频率响应# 这对理解滤波器做了什么非常有用filter_params = mne.filter.create_filter(    raw.get_data(), sfreq=sfreq,    l_freq=4.0, h_freq=12.0,    method='fir')mne.viz.plot_filter(filter_params, sfreq=sfreq)

4.5 对 Epochs 滤波

# 也可以对 Epochs 滤波（不推荐，应在 Raw 阶段做）epochs_filt = epochs.copy().filter(l_freq=1.0, h_freq=40.0)

先滤波还是先分段？必须先滤波再分段。 如果先 epoch 再滤波，epoch 边缘会产生严重的滤波伪迹（edge artifact），尤其是慢波段（theta, delta）。这是 MNE 文档反复强调的原则。

5. 预处理：重参考

LFP 信号是两点之间的电位差。重参考 = 换一个参考点。

5.1 常见参考方案

# === 1. 共同平均参考 (CAR) ===# 对 seeg 通道做平均参考raw_car, _ = mne.set_eeg_reference(    raw.copy(), ref_channels='average', ch_type='seeg')# === 2. 指定通道参考 ===# 用某个特定通道作为参考（比如一个安静的白质通道）raw_ref, _ = mne.set_eeg_reference(    raw.copy(), ref_channels=['VTA-8'], ch_type='seeg')# === 3. 区域平均参考 ===# 每个脑区内部做 CAR# MNE 没有内置函数，但很容易实现raw_regional = raw.copy()for region, chs in regions.items():    data_region = raw_regional.get_data(picks=chs)    mean_ref = data_region.mean(axis=0)for ch in chs:        ch_idx = raw_regional.ch_names.index(ch)        raw_regional._data[ch_idx] -= mean_ref

5.2 双极参考（LFP 最常用）

双极参考 = 相邻通道相减。它能有效去除共模噪声，增强局部信号。

# === MNE 内置双极参考 ===# 计算相邻通道差：ch1-ch2, ch2-ch3, ch3-ch4, ...anode = [f'HPC-{i}'for i in range(1, 8)]    # HPC-1 到 HPC-7cathode = [f'HPC-{i}'for i in range(2, 9)]  # HPC-2 到 HPC-8raw_bipolar = mne.set_bipolar_reference(    raw.copy(),    anode=anode,    cathode=cathode,    ch_name=[f'HPC-{i}{i+1}'for i in range(1, 8)]  # 新通道名)print(raw_bipolar.ch_names)# ['mPFC-1', ..., 'HPC-12', 'HPC-23', ..., 'HPC-78', ...]

LFP 重参考方案怎么选？

方案	适用场景	优点	缺点
双极	线性探针（Neuropixels, 硅探针）	去共模噪声、增强局部信号	通道数减 1
CAR	同脑区内多电极	简单	脑区间信号差异大时效果差
白质参考	有白质电极时	参考信号稳定	需要确认参考电极确实在白质
不参考	录制时已有良好参考	不引入新伪迹	—

6. 预处理：伪迹检测与去除

6.1 基于幅度的自动拒绝

# === 在 Epochs 创建时设置 reject 阈值 ===reject = dict(seeg=500e-6)   # 拒绝振幅超过 500 µV 的 epochflat = dict(seeg=1e-7)       # 拒绝振幅小于 0.1 µV 的 epoch（可能是死通道）epochs = mne.Epochs(raw, events, event_id=event_dict,                    tmin=-0.5, tmax=1.0,                    reject=reject, flat=flat,                    preload=True)print(f'保留 {len(epochs)} / {len(events)} 个 epoch')print(epochs.drop_log_stats())  # 丢弃比例

6.2 基于 Annotations 的手动标记

# === 手动标记伪迹段 ===# 比如动物运动导致的噪声段bad_segments = mne.Annotations(    onset=[30.0, 85.0, 150.0],    duration=[5.0, 10.0, 3.0],    description=['bad_movement', 'bad_saturation', 'bad_noise'])raw.set_annotations(bad_segments)# 在后续 epoching 时，与 'bad_*' 重叠的 epoch 自动被丢弃

6.3 ICA 去伪迹

ICA（独立成分分析）在 LFP 中不如在 EEG 中常用，但在某些场景下有用（比如去除心电伪迹、呼吸伪迹）。

from mne.preprocessing import ICA# === 设置 ICA ===ica = ICA(    n_components=15,       # 提取的成分数（小于通道数）    method='fastica',      # 'fastica', 'infomax', 'picard'    random_state=42)# 先滤波再做 ICA（ICA 对低频漂移敏感）raw_for_ica = raw.copy().filter(l_freq=1.0, h_freq=None)ica.fit(raw_for_ica)# === 查看 ICA 成分 ===# 时间序列ica.plot_sources(raw_for_ica, show_scrollbars=False)# 各成分的频谱ica.plot_properties(raw_for_ica, picks=[0, 1, 2])# === 标记并排除伪迹成分 ===# 方法 1：手动检查后标记ica.exclude = [0, 3]  # 排除第 0 和第 3 个成分# 方法 2：用参考信号自动检测# 如果你有同步记录的心电或呼吸信号ecg_indices, ecg_scores = ica.find_bads_ecg(raw, ch_name='ECG-1')ica.exclude.extend(ecg_indices)# === 应用 ICA（去除伪迹成分）===raw_clean = ica.apply(raw.copy())

7. 重采样

高采样率（如 30 kHz）对 LFP 分析通常是多余的（LFP 频段一般关注 <300 Hz）。重采样可以大幅减小数据量。

# === 降采样到 1000 Hz ===raw_downsampled = raw.copy().resample(sfreq=1000.0, n_jobs=-1)print(f'降采样后: {raw_downsampled.n_times} 点, 采样率 {raw_downsampled.info["sfreq"]} Hz')# === 先低通滤波再降采样（推荐）===# MNE 的 resample() 内部会自动做抗混叠低通滤波# 但如果你的目标采样率较低，建议手动先滤波raw_for_ds = raw.copy().filter(l_freq=None, h_freq=400.0)  # 低通 400 Hzraw_ds = raw_for_ds.resample(sfreq=1000.0)# === 对 Epochs 降采样 ===epochs_ds = epochs.copy().resample(sfreq=250.0)# 如果事件时序精度很关键，官方更推荐先 epoch 再降采样；# 若直接对 Raw 降采样，可用 events=... 同步重采样事件采样点

目标采样率应至少是你关心的最高频率的 2.5 倍（Nyquist + 余量）：

• 只关心 theta/gamma (<100 Hz) → 250 Hz 够用
• 关心 ripple (<250 Hz) → 1000 Hz
• 关心 fast ripple (<600 Hz) → 1500 Hz
• Spike sorting → 不要降采样，保持原始采样率

8. 频谱分析

8.1 计算功率谱密度 (PSD)

MNE 提供两种 PSD 估计方法：Welch 法和 Multitaper 法。

# === Welch 法 ===spectrum = raw.compute_psd(    method='welch',    fmin=1.0,          # 最低频率    fmax=300.0,        # 最高频率    n_fft=2048,        # FFT 长度    n_overlap=1024,    # 重叠点数    picks='seeg')# spectrum 是一个 Spectrum 对象psds, freqs = spectrum.get_data(return_freqs=True)# psds.shape = (n_channels, n_freqs)  单位：V²/Hz# freqs.shape = (n_freqs,)# === Multitaper 法 ===spectrum_mt = raw.compute_psd(    method='multitaper',    fmin=1.0,    fmax=300.0,    bandwidth=2.0,     # 频率平滑带宽 (Hz)    picks='seeg')

	Welch	Multitaper
原理	分段取平均	多个正交窗取平均
频率分辨率	取决于 n_fft	取决于 bandwidth
方差	需要较多数据	同等数据长度下方差更小
适用	长时间连续记录	短数据段或需要精确估计

8.2 PSD 可视化

# === 内置绘图 ===spectrum.plot(picks='seeg', amplitude=False)  # amplitude=False → 画功率不画幅度# === 自定义绘图（更灵活）===import matplotlib.pyplot as pltpsds_db = 10 * np.log10(psds)  # 转换为 dBfig, ax = plt.subplots(figsize=(10, 5))for i, ch_name in enumerate(spectrum.ch_names):if ch_name.startswith('HPC'):        ax.plot(freqs, psds_db[i], label=ch_name, alpha=0.7)ax.set(xlabel='Frequency (Hz)', ylabel='Power (dB, V²/Hz)',       title='HPC channels PSD', xlim=(1, 100))ax.legend(fontsize=8)ax.set_xscale('log')plt.tight_layout()plt.show()

8.3 Epochs 的 PSD

# 对每个 epoch 分别计算 PSDspectrum_epochs = epochs.compute_psd(    method='welch',    fmin=1.0, fmax=100.0,    n_fft=512,    picks='seeg')# spectrum_epochs.get_data().shape = (n_epochs, n_channels, n_freqs)# 按条件比较 PSDspectrum_light = epochs['light_stim'].compute_psd(method='welch', fmin=1, fmax=100)spectrum_sound = epochs['sound_stim'].compute_psd(method='welch', fmin=1, fmax=100)psd_light, freqs = spectrum_light.get_data(return_freqs=True)psd_sound, _ = spectrum_sound.get_data(return_freqs=True)# 按通道平均后比较psd_light_mean = psd_light.mean(axis=(0, 1))  # 跨 epoch 和通道平均psd_sound_mean = psd_sound.mean(axis=(0, 1))fig, ax = plt.subplots()ax.semilogy(freqs, psd_light_mean, label='Light stim')ax.semilogy(freqs, psd_sound_mean, label='Sound stim')ax.set(xlabel='Frequency (Hz)', ylabel='PSD (V²/Hz)')ax.legend()plt.show()

8.4 提取频段功率

# 定义频段bands = {'delta': (1, 4),'theta': (4, 12),'beta':  (12, 30),'gamma': (30, 100),}# 计算各频段的平均功率psds, freqs = spectrum.get_data(return_freqs=True)for band_name, (fmin, fmax) in bands.items():    freq_mask = (freqs >= fmin) & (freqs <= fmax)    band_power = psds[:, freq_mask].mean(axis=1)  # 每个通道的频段平均功率    print(f'{band_name}: {band_power.mean():.2e} V²/Hz')

9. 时频分析

9.1 Morlet 小波时频分解

import mne.time_frequency as tfr# 定义频率轴freqs_tfr = np.logspace(np.log10(4), np.log10(100), 30)  # 4-100 Hz, 30 个频率n_cycles = freqs_tfr / 2.0# 每个频率的小波周期数（频率/2 是常用选择）# === 对 Epochs 做时频分解 ===power = tfr.tfr_morlet(    epochs,    freqs=freqs_tfr,    n_cycles=n_cycles,    return_itc=False,     # True 则同时返回 inter-trial coherence    picks='seeg',    n_jobs=-1,    average=True# 跨 trial 平均)# power 是一个 AverageTFR 对象# power.data.shape = (n_channels, n_freqs, n_times)print(power)

n_cycles 的选择n_cycles 控制时间-频率分辨率的权衡：

• 固定值（如 n_cycles=7）：所有频率使用相同周期数，低频时间分辨率差
• 频率成比例（如 n_cycles=freqs/2）：低频用少周期（好的时间分辨率），高频用多周期（好的频率分辨率）
• LFP 分析中推荐 freqs/2 到 `freqs/4

9.2 时频图可视化

# === 内置绘图 ===power.plot(picks='HPC-1', baseline=(-0.5, -0.1),           mode='percent',     # 基线校正模式           title='HPC-1 时频图')# mode 选项：# 'mean'     → (power - baseline_mean)# 'ratio'    → power / baseline_mean# 'logratio' → log10(power / baseline_mean)# 'percent'  → (power - baseline_mean) / baseline_mean × 100# 'zscore'   → (power - baseline_mean) / baseline_std# === 自定义绘图 ===ch_idx = power.ch_names.index('HPC-1')tfr_data = power.data[ch_idx]  # (n_freqs, n_times)# 基线校正baseline_mask = power.times < 0baseline_mean = tfr_data[:, baseline_mask].mean(axis=1, keepdims=True)tfr_db = 10 * np.log10(tfr_data / baseline_mean)fig, ax = plt.subplots(figsize=(10, 5))im = ax.pcolormesh(power.times, freqs_tfr, tfr_db,                   cmap='RdBu_r', vmin=-3, vmax=3, shading='auto')ax.set(xlabel='Time (s)', ylabel='Frequency (Hz)', title='HPC-1 Time-Frequency')ax.set_yscale('log')ax.axvline(0, color='k', linestyle='--', label='Stimulus')fig.colorbar(im, ax=ax, label='Power (dB re baseline)')plt.tight_layout()plt.show()

9.3 Multitaper 时频分解

# Multitaper 比 Morlet 的频率估计方差更小power_mt = tfr.tfr_multitaper(    epochs,    freqs=freqs_tfr,    n_cycles=freqs_tfr / 2.0,    time_bandwidth=4.0,    # NW 时间带宽积 → 使用 2*NW-1 = 7 个锥度    return_itc=False,    picks='seeg',    n_jobs=-1,    average=True)

9.4 单 trial 时频（不平均）

# average=False → 返回每个 trial 的时频power_single = tfr.tfr_morlet(    epochs,    freqs=freqs_tfr,    n_cycles=freqs_tfr / 2.0,    return_itc=False,    picks='seeg',    average=False,    n_jobs=-1)# power_single.data.shape = (n_epochs, n_channels, n_freqs, n_times)# 用处：# - 检查 trial-to-trial 变异性# - 做单 trial 分析（如与行为数据关联）# - 先 subtract_evoked() 再做 TFR，得到 induced power

9.5 诱发 (Evoked) vs 诱导 (Induced) 功率

# Evoked power = 先平均再做 TFR（锁相成分）evoked = epochs.average()power_evoked = tfr.tfr_morlet(    evoked, freqs=freqs_tfr, n_cycles=freqs_tfr/2,    return_itc=False, picks='seeg')# Total power = 先做 TFR 再平均（所有成分）power_total = tfr.tfr_morlet(    epochs, freqs=freqs_tfr, n_cycles=freqs_tfr/2,    return_itc=False, picks='seeg', average=True)# Induced power = 先从每个 trial 中减去 evoked，再做 TFRepochs_induced = epochs.copy().subtract_evoked()power_induced = tfr.tfr_morlet(    epochs_induced, freqs=freqs_tfr, n_cycles=freqs_tfr/2,    return_itc=False, picks='seeg', average=True)

10. 统计检验

10.1 配对设计的簇排列检验

from mne.stats import permutation_cluster_1samp_test# === 比较两个条件的 PSD（配对设计：同一批 epoch）===psd_light, freqs = epochs['light_stim'].compute_psd(    method='welch', fmin=1, fmax=100).get_data(return_freqs=True)psd_sound, _ = epochs['sound_stim'].compute_psd(    method='welch', fmin=1, fmax=100).get_data(return_freqs=True)# 取某个通道ch_idx = 0X = psd_light[:, ch_idx, :] - psd_sound[:, ch_idx, :]# X.shape = (n_epochs, n_freqs)# 对“条件差值”做单样本簇排列检验T_obs, clusters, cluster_p, H0 = permutation_cluster_1samp_test(    X,    n_permutations=1000,    threshold=2.0,        # t 值阈值    tail=0,               # 0=双侧, 1=单侧(A>B), -1=单侧(A<B)    out_type='mask',    n_jobs=-1)# 找显著簇for i_cl, (cl, p) in enumerate(zip(clusters, cluster_p)):if p < 0.05:        freq_range = freqs[cl]        print(f'簇 {i_cl}: {freq_range[0]:.1f}-{freq_range[-1]:.1f} Hz, p={p:.4f}')

10.2 时频数据的簇检验

# === 对时频数据做配对簇检验 ===# 需要 single-trial TFR 数据power_light = tfr.tfr_morlet(    epochs['light_stim'], freqs=freqs_tfr, n_cycles=freqs_tfr/2,    return_itc=False, picks=['HPC-1'], average=False).data[:, 0, :, :]  # (n_epochs, n_freqs, n_times)power_sound = tfr.tfr_morlet(    epochs['sound_stim'], freqs=freqs_tfr, n_cycles=freqs_tfr/2,    return_itc=False, picks=['HPC-1'], average=False).data[:, 0, :, :]# 配对设计：先做条件差值，再做 2D 单样本簇检验（频率 × 时间）power_diff = power_light - power_soundT_obs, clusters, cluster_p, H0 = permutation_cluster_1samp_test(    power_diff,    n_permutations=1000,    threshold=2.0,    tail=0,    out_type='mask',    n_jobs=-1)# 可视化显著簇sig_mask = np.zeros_like(T_obs, dtype=bool)for cl, p in zip(clusters, cluster_p):if p < 0.05:        sig_mask |= clpower_times = epochs['light_stim'].timesfig, axes = plt.subplots(1, 2, figsize=(14, 5))axes[0].pcolormesh(power_times, freqs_tfr, T_obs, cmap='RdBu_r', shading='auto')axes[0].set_title('T-statistics')axes[1].pcolormesh(power_times, freqs_tfr, sig_mask.astype(float),                   cmap='Reds', shading='auto')axes[1].set_title('Significant clusters (p<0.05)')for ax in axes:    ax.set(xlabel='Time (s)', ylabel='Frequency (Hz)')plt.tight_layout()plt.show()

10.3 单样本检验（对基线）

# 检验某个条件的功率是否显著不同于基线power_data = tfr.tfr_morlet(    epochs['light_stim'], freqs=freqs_tfr, n_cycles=freqs_tfr/2,    return_itc=False, picks=['HPC-1'], average=False).data[:, 0, :, :]  # (n_epochs, n_freqs, n_times)# 基线归一化：每个 trial 减去基线均值power_times = epochs['light_stim'].timesbaseline_mask = power_times < 0baseline_mean = power_data[:, :, baseline_mask].mean(axis=2, keepdims=True)power_norm = power_data / baseline_mean - 1# percent change# 单样本排列簇检验（检验是否显著不同于 0）T_obs, clusters, cluster_p, H0 = permutation_cluster_1samp_test(    power_norm,    n_permutations=1000,    threshold=2.0,    tail=1,          # 只检验增强（>0）    out_type='mask',    n_jobs=-1)

11. 可视化

11.1 连续数据浏览

# === 交互式浏览 ===raw.plot(    duration=10.0,        # 显示窗口宽度 (秒)    n_channels=16,        # 同时显示的通道数    scalings=dict(seeg=100e-6),  # y 轴缩放 (100 µV)    title='LFP recording',    show_scrollbars=True,    block=True# 阻塞直到关闭窗口)# 交互操作：# - 点击通道名 → 标记为坏道# - 拖选时间段 → 标记为 bad_* annotation# - 方向键 → 滚动

11.2 Epochs 浏览

epochs.plot(    n_channels=16,    n_epochs=5,           # 同时显示几个 epoch    scalings=dict(seeg=100e-6),    title='LFP Epochs',    block=True)# 点击某个 epoch → 标记为 bad（丢弃）

11.3 Evoked 绘图

# === 所有通道叠加 ===evoked_light.plot(    spatial_colors=False,    titles='LFP Evoked (light stim)',    time_unit='ms')# === 指定通道 ===evoked_light.plot(picks=['HPC-1', 'HPC-2', 'mPFC-1'])# === 多条件对比 ===mne.viz.plot_compare_evokeds(    {'light': evoked_light, 'sound': evoked_sound},    picks=['HPC-1'],    title='HPC-1: Light vs Sound')

11.4 PSD 绘图

# === 内置方法 ===spectrum.plot(picks='seeg')# === 按脑区分面板 ===fig, axes = plt.subplots(2, 2, figsize=(12, 8))for ax, (region, chs) in zip(axes.flat, regions.items()):    psds_region, freqs = spectrum.get_data(picks=chs, return_freqs=True)    psds_db = 10 * np.log10(psds_region)    ax.plot(freqs, psds_db.T, alpha=0.5)    ax.plot(freqs, psds_db.mean(axis=0), 'k-', linewidth=2, label='Mean')    ax.set(title=region, xlabel='Frequency (Hz)', ylabel='PSD (dB)')    ax.set_xlim(1, 100)    ax.set_xscale('log')    ax.legend()plt.tight_layout()plt.show()

11.5 时频图

# 参见 §9.2 的示例# === 多通道时频对比 ===fig, axes = plt.subplots(2, 2, figsize=(14, 10))for ax, ch_name in zip(axes.flat, ['mPFC-1', 'HPC-1', 'AMG-1', 'VTA-1']):    ch_idx = power.ch_names.index(ch_name)    tfr_data = power.data[ch_idx]    baseline_mask = power.times < 0    baseline_mean = tfr_data[:, baseline_mask].mean(axis=1, keepdims=True)    tfr_db = 10 * np.log10(tfr_data / baseline_mean)    im = ax.pcolormesh(power.times, freqs_tfr, tfr_db,                       cmap='RdBu_r', vmin=-3, vmax=3, shading='auto')    ax.set(title=ch_name, xlabel='Time (s)', ylabel='Frequency (Hz)')    ax.set_yscale('log')    ax.axvline(0, color='k', linestyle='--')    fig.colorbar(im, ax=ax, shrink=0.6)plt.tight_layout()plt.show()

12. 实战 Pipeline 模板

一个完整的 LFP 分析流水线示例：

"""LFP 分析 Pipeline 模板━━━━━━━━━━━━━━━━━━━━━实验：光遗传刺激大鼠 mPFC，记录 mPFC + dHPC LFP电极：32 通道硅探针 (mPFC 16ch + dHPC 16ch)采样率：30 kHz分析目标：刺激前后 theta/gamma 功率变化"""import numpy as npimport mneimport mne.time_frequency as tfrfrom mne.stats import permutation_cluster_1samp_test# ============================================================# 1. 加载数据# ============================================================# 从 .mat 文件读取（你的数据格式可能不同）import scipy.io as siomat = sio.loadmat('optostim_session01.mat')lfp_raw = mat['lfp'] * 1e-6# (32, 18000000) 10 分钟 @ 30 kHz, 转为 Vsfreq_orig = 30000.0stim_times = mat['stim_times'].flatten()  # 光刺激时间戳 (秒)# 通道名ch_names = [f'mPFC-{i}'for i in range(1, 17)] + \           [f'dHPC-{i}'for i in range(1, 17)]info = mne.create_info(ch_names, sfreq_orig, ch_types='seeg')raw = mne.io.RawArray(lfp_raw, info)print(f'原始数据: {raw.n_times} 点, {raw.times[-1]:.1f} s, {sfreq_orig} Hz')# ============================================================# 2. 预处理# ============================================================# 2a. 坏道标记raw.info['bads'] = ['mPFC-8', 'dHPC-12']  # 通过目视检查确定# 2b. 陷波滤波 (去工频)raw.notch_filter(freqs=[50, 100, 150, 200, 250], picks='seeg')# 2c. 带通滤波raw.filter(l_freq=0.5, h_freq=300.0, picks='seeg')# 2d. 降采样 (30 kHz → 1000 Hz)raw.resample(sfreq=1000.0, n_jobs=-1)sfreq = raw.info['sfreq']# 2e. 双极参考 (mPFC 和 dHPC 分别做)for prefix, n_ch in [('mPFC', 16), ('dHPC', 16)]:    good_chs = [f'{prefix}-{i}'for i in range(1, n_ch+1)iff'{prefix}-{i}'notin raw.info['bads']]if len(good_chs) > 1:        anode = good_chs[:-1]        cathode = good_chs[1:]        new_names = [f'{a}-{c.split("-")[1]}'for a, c in zip(anode, cathode)]        raw = mne.set_bipolar_reference(raw, anode, cathode, ch_name=new_names)# ============================================================# 3. 分段 (Epoching)# ============================================================events = np.column_stack([    (stim_times * sfreq).astype(int),    np.zeros(len(stim_times), dtype=int),    np.ones(len(stim_times), dtype=int)])event_dict = {'opto_stim': 1}epochs = mne.Epochs(    raw, events, event_id=event_dict,    tmin=-2.0, tmax=3.0,    baseline=None,          # 不在这里做基线校正（留给后续分析）    reject=dict(seeg=500e-6),    preload=True)print(f'保留 {len(epochs)} / {len(stim_times)} 个 epoch')# ============================================================# 4. 频谱分析# ============================================================# 4a. PSD 对比（刺激前 vs 刺激后）epochs_pre = epochs.copy().crop(tmin=-2.0, tmax=0.0)epochs_post = epochs.copy().crop(tmin=0.5, tmax=2.5)psd_pre, freqs = epochs_pre.compute_psd(    method='welch', fmin=1, fmax=100, n_fft=1024).get_data(return_freqs=True)psd_post, _ = epochs_post.compute_psd(    method='welch', fmin=1, fmax=100, n_fft=1024).get_data(return_freqs=True)# 4b. 时频分析freqs_tfr = np.logspace(np.log10(4), np.log10(100), 30)power = tfr.tfr_morlet(    epochs, freqs=freqs_tfr, n_cycles=freqs_tfr/2,    return_itc=False, picks='seeg', average=True, n_jobs=-1)# ============================================================# 5. 统计检验# ============================================================# 频段功率变化的显著性检验bands = {'theta': (4, 12), 'gamma': (30, 100)}for band_name, (fmin, fmax) in bands.items():    freq_mask = (freqs >= fmin) & (freqs <= fmax)    power_pre_band = psd_pre[:, :, freq_mask].mean(axis=2)   # (n_epochs, n_channels)    power_post_band = psd_post[:, :, freq_mask].mean(axis=2)# 配对检验：每个通道from scipy import statsfor ch_idx, ch_name in enumerate(epochs.ch_names):        t, p = stats.ttest_rel(power_post_band[:, ch_idx],                                power_pre_band[:, ch_idx])if p < 0.05:            ratio = power_post_band[:, ch_idx].mean() / power_pre_band[:, ch_idx].mean()            print(f'{ch_name}{band_name}: {ratio:.2f}x change, p={p:.4f}')# ============================================================# 6. 保存# ============================================================epochs.save('optostim_session01-epo.fif', overwrite=True)power.save('optostim_session01-tfr.h5', overwrite=True)

13. API 速查表

创建数据

读取数据

预处理

分段与平均

频谱分析

统计

可视化