用Python实战脑电分析：手把手教你计算PLV、MVL、MI跨频耦合指标（附完整代码）
用Python实战脑电分析：手把手教你计算PLV、MVL、MI跨频耦合指标（附完整代码）
当你第一次看到脑电图(EEG)数据时，那些看似杂乱的波形可能令人望而生畏。但正是这些波动中隐藏着大脑活动的奥秘——不同频段的神经振荡如何相互影响。作为神经科学研究的重要工具，跨频耦合分析能揭示θ波与γ波等不同频段间的互动规律。本文将带你用Python从零实现三种主流算法，无需复杂公式推导，直接上手实操。
1. 环境准备与数据模拟
工欲善其事，必先利其器。我们先搭建一个可复现的分析环境：
# 基础环境配置 import numpy as np import matplotlib.pyplot as plt from scipy.signal import hilbert, butter, filtfilt import random # 确保结果可复现 np.random.seed(42)模拟数据是理解算法的绝佳起点。我们构建一个4Hz低频相位调制32Hz高频振幅的合成信号：
def generate_cfc_signal(duration=5, fs=1000): """生成具有跨频耦合特性的模拟信号""" t = np.arange(0, duration, 1/fs) # 低频相位成分 (4Hz) low_freq = 4 phase_signal = np.sin(2*np.pi*low_freq*t) # 高频振幅成分 (32Hz)，受低频调制 high_freq = 32 modulated_amp = 1 + 0.8*np.sin(2*np.pi*low_freq*t) # 调制深度0.8 amp_signal = modulated_amp * np.sin(2*np.pi*high_freq*t) return phase_signal + amp_signal, t eeg_sim, time = generate_cfc_signal()表：模拟信号参数说明
| 参数 | 值 | 作用 |
|---|---|---|
| duration | 5秒 | 信号时长 |
| fs | 1000Hz | 采样率 |
| low_freq | 4Hz | 相位载波频率 |
| high_freq | 32Hz | 振幅载波频率 |
| modulation_depth | 0.8 | 耦合强度 |
提示：实际分析时，建议先对原始EEG进行预处理，包括去噪、去除眼电等步骤，这里为简化流程直接使用模拟数据。
2. 核心算法实现与对比
2.1 锁相值(PLV)计算
PLV通过评估相位同步性来量化耦合强度，其值域为[0,1]，越接近1表示耦合越强：
def plv_analysis(phase_signal, amp_signal): """计算相位-振幅锁相值""" # 提取相位和振幅 phase = np.angle(hilbert(phase_signal)) amp_phase = np.angle(hilbert(amp_signal)) # 核心计算公式 plv = np.abs(np.mean(np.exp(1j*(phase - amp_phase)))) return plv # 应用示例 phase_filt = butter_bandpass_filter(eeg_sim, lowcut=3, highcut=7, fs=1000) amp_filt = butter_bandpass_filter(eeg_sim, lowcut=30, highcut=34, fs=1000) plv_value = plv_analysis(phase_filt, amp_filt) print(f"PLV计算结果: {plv_value:.3f}")PLV特点：
- 对噪声鲁棒性强
- 需要窄带滤波
- 可能低估非线性耦合
2.2 平均向量长度(MVL)实现
MVL将振幅信息纳入计算，适合信噪比较高的数据：
def mvl_analysis(phase_signal, amp_signal): """计算平均向量长度""" phase = np.angle(hilbert(phase_signal)) amp = np.abs(hilbert(amp_signal)) # 振幅归一化 amp_norm = (amp - np.min(amp)) / (np.max(amp) - np.min(amp)) # 向量平均 mvl = np.abs(np.mean(amp_norm * np.exp(1j*phase))) return mvl mvl_value = mvl_analysis(phase_filt, amp_filt) print(f"MVL计算结果: {mvl_value:.3f}")表：PLV与MVL对比
| 指标 | 计算重点 | 适用场景 | 敏感度 |
|---|---|---|---|
| PLV | 纯相位关系 | 低信噪比数据 | 相位同步 |
| MVL | 相位+振幅 | 高质量数据 | 振幅调制 |
2.3 调制指数(MI)全流程
MI通过熵分析评估耦合强度，是最常用的方法：
def mi_analysis(phase_signal, amp_signal, n_bins=18): """计算调制指数""" phase = np.angle(hilbert(phase_signal)) amp = np.abs(hilbert(amp_signal)) # 相位分箱 bins = np.linspace(-np.pi, np.pi, n_bins+1) mean_amps = [] for i in range(n_bins): # 获取每个相位区间的振幅均值 mask = (phase >= bins[i]) & (phase < bins[i+1]) mean_amps.append(np.mean(amp[mask])) # 计算归一化概率分布 prob = mean_amps / np.sum(mean_amps) entropy = -np.sum(prob * np.log(prob + 1e-10)) # 避免log(0) # 最终MI计算 mi = (np.log(n_bins) - entropy) / np.log(n_bins) return mi, mean_amps mi_value, amp_dist = mi_analysis(phase_filt, amp_filt) print(f"MI计算结果: {mi_value:.3f}")可视化能直观验证结果：
plt.figure(figsize=(10,4)) plt.subplot(121) plt.plot(amp_dist) plt.title('振幅相位分布') plt.subplot(122) plt.polar(np.linspace(0, 2*np.pi, 18), amp_dist) plt.title('极坐标显示') plt.tight_layout()3. 统计检验与结果解读
3.1 置换检验实现
为避免假阳性，需评估结果的统计显著性：
def permutation_test(phase, amp, func, n_iter=500): """置换检验评估显著性""" observed = func(phase, amp) null_dist = np.zeros(n_iter) for i in range(n_iter): # 随机打乱振幅序列 shift = np.random.randint(len(amp)//10, len(amp)*9//10) shuffled_amp = np.roll(amp, shift) null_dist[i] = func(phase, shuffled_amp) # 计算z分数 z_score = (observed - np.mean(null_dist)) / np.std(null_dist) p_value = np.sum(null_dist >= observed) / n_iter return z_score, p_value plv_z, plv_p = permutation_test(phase_filt, amp_filt, plv_analysis) print(f"PLV显著性: z={plv_z:.2f}, p={plv_p:.4f}")注意：实际应用中建议n_iter≥1000，这里为演示设置为500次
3.2 结果解读指南
- PLV/MVL值：>0.3通常认为存在显著耦合
- MI值：>0.05可能具有生理意义
- z分数：>1.96对应p<0.05显著性水平
表：典型耦合模式参考
| 耦合类型 | 频段组合 | 常见脑区 | 认知功能 |
|---|---|---|---|
| θ-γ耦合 | 4-8Hz & 30-80Hz | 前额叶 | 工作记忆 |
| α-γ耦合 | 8-12Hz & 30-100Hz | 枕叶 | 视觉注意 |
| β-γ耦合 | 13-30Hz & 50-150Hz | 运动区 | 运动准备 |
4. 真实EEG分析实战
4.1 数据加载与预处理
使用MNE库处理真实EEG数据：
import mne # 示例数据加载 sample_data = mne.datasets.sample.data_path() raw = mne.io.read_raw_fif(sample_data + '/MEG/sample/sample_audvis_raw.fif', preload=True) # 选取电极并滤波 picks = mne.pick_types(raw.info, eeg=True) raw.filter(4, 30, picks=picks) # θ-β频段4.2 全频段耦合分析
自动化分析不同频段组合：
from tqdm import tqdm def cfc_matrix(data, fs, phase_range=(4,8), amp_range=(30,80), n_bins=18): """计算频段间耦合矩阵""" phase_freqs = np.arange(*phase_range) amp_freqs = np.arange(*amp_range) result = np.zeros((len(phase_freqs), len(amp_freqs))) for i, pf in tqdm(enumerate(phase_freqs)): # 提取相位信号 phase_sig = butter_bandpass_filter(data, pf-1, pf+1, fs) phase = np.angle(hilbert(phase_sig)) for j, af in enumerate(amp_freqs): # 提取振幅信号 amp_sig = butter_bandpass_filter(data, af-5, af+5, fs) amp = np.abs(hilbert(amp_sig)) # 计算MI result[i,j], _ = mi_analysis(phase, amp, n_bins) return result # 示例使用 cfc_result = cfc_matrix(raw.get_data()[0], raw.info['sfreq'])4.3 结果可视化技巧
热图展示耦合强度分布：
plt.figure(figsize=(10,6)) plt.imshow(cfc_result, aspect='auto', extent=[30,80,8,4], cmap='jet') plt.colorbar(label='Modulation Index') plt.xlabel('Amplitude Frequency (Hz)') plt.ylabel('Phase Frequency (Hz)') plt.title('Cross-Frequency Coupling Matrix')5. 工程实践中的常见问题
5.1 数据质量检查清单
- 采样率是否足够（≥2×最高分析频率）
- 是否进行了适当的带通滤波
- 信号是否包含明显伪迹（眼动、肌电等）
- 数据长度是否足够（建议≥10个低频周期）
5.2 参数选择经验法则
| 参数 | 推荐值 | 调整建议 |
|---|---|---|
| 相位频带带宽 | ±1-2Hz | 过宽会降低特异性 |
| 振幅频带带宽 | ±5-10Hz | γ频段可适当放宽 |
| 相位分箱数 | 12-36 | 权衡分辨率与稳定性 |
| 置换检验次数 | ≥1000 | 计算资源允许时越多越好 |
5.3 性能优化技巧
# 使用numba加速计算 from numba import jit @jit(nopython=True) def fast_mi(phase, amp, n_bins): # 实现略... return mi # 并行化处理多个通道 from joblib import Parallel, delayed def parallel_cfc(data): return Parallel(n_jobs=4)(delayed(cfc_analysis)(ch) for ch in data)6. 完整项目架构
推荐的项目目录结构：
/eeg_cfc_analysis │── /data # 原始数据 │ ├── raw_eeg.edf │ └── metadata.csv │── /notebooks # Jupyter笔记本 │ └── cfc_demo.ipynb │── /src # Python模块 │ ├── preprocessing.py # 预处理函数 │ ├── cfc_methods.py # 核心算法 │ └── visualization.py # 绘图工具 │── requirements.txt # 依赖库 └── run_analysis.py # 主入口典型分析流程：
- 数据预处理（去噪、滤波）
- 选择目标频段组合
- 计算耦合指标
- 统计显著性检验
- 结果可视化与解释
# 示例主程序框架 def main(): # 1. 加载数据 raw = load_eeg('data/example.edf') # 2. 预处理 raw.filter(1, 50) raw = remove_artifacts(raw) # 3. 分析指定频段 results = [] for phase_band in [(4,8), (8,12)]: for amp_band in [(30,50), (50,80)]: res = analyze_band_combination(raw, phase_band, amp_band) results.append(res) # 4. 可视化 generate_report(results)7. 进阶应用方向
7.1 时变耦合分析
通过滑动窗口研究耦合动态变化：
def time_varying_cfc(data, fs, window_size=1.0, step=0.1): """计算时变跨频耦合""" n_samples = len(data) window_points = int(window_size * fs) step_points = int(step * fs) results = [] for start in range(0, n_samples-window_points, step_points): segment = data[start:start+window_points] mi = mi_analysis(segment, fs) results.append(mi) return np.array(results)7.2 多通道耦合网络
研究不同脑区间的耦合模式：
def network_cfc(raw, ch_pairs): """计算通道间的耦合网络""" n_channels = len(raw.ch_names) adj_matrix = np.zeros((n_channels, n_channels)) for i, j in ch_pairs: data_i = raw.get_data(picks=i) data_j = raw.get_data(picks=j) adj_matrix[i,j] = plv_analysis(data_i, data_j) return adj_matrix7.3 机器学习整合
将CFC特征用于分类任务：
from sklearn.ensemble import RandomForestClassifier def extract_cfc_features(X, fs): """从EEG片段提取CFC特征""" features = [] for phase_band in [(4,8), (8,12)]: for amp_band in [(30,50), (50,80)]: mi = band_analysis(X, fs, phase_band, amp_band) features.append(mi) return np.array(features) # 示例分类流程 X_train = [extract_cfc_features(x) for x in train_data] clf = RandomForestClassifier().fit(X_train, y_train)