用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(duration5, fs1000): 生成具有跨频耦合特性的模拟信号 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()表模拟信号参数说明参数值作用duration5秒信号时长fs1000Hz采样率low_freq4Hz相位载波频率high_freq32Hz振幅载波频率modulation_depth0.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, lowcut3, highcut7, fs1000) amp_filt butter_bandpass_filter(eeg_sim, lowcut30, highcut34, fs1000) plv_value plv_analysis(phase_filt, amp_filt) print(fPLV计算结果: {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(fMVL计算结果: {mvl_value:.3f})表PLV与MVL对比指标计算重点适用场景敏感度PLV纯相位关系低信噪比数据相位同步MVL相位振幅高质量数据振幅调制2.3 调制指数(MI)全流程MI通过熵分析评估耦合强度是最常用的方法def mi_analysis(phase_signal, amp_signal, n_bins18): 计算调制指数 phase np.angle(hilbert(phase_signal)) amp np.abs(hilbert(amp_signal)) # 相位分箱 bins np.linspace(-np.pi, np.pi, n_bins1) mean_amps [] for i in range(n_bins): # 获取每个相位区间的振幅均值 mask (phase bins[i]) (phase bins[i1]) 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(fMI计算结果: {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_iter500): 置换检验评估显著性 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(fPLV显著性: z{plv_z:.2f}, p{plv_p:.4f})注意实际应用中建议n_iter≥1000这里为演示设置为500次3.2 结果解读指南PLV/MVL值0.3通常认为存在显著耦合MI值0.05可能具有生理意义z分数1.96对应p0.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, preloadTrue) # 选取电极并滤波 picks mne.pick_types(raw.info, eegTrue) raw.filter(4, 30, pickspicks) # θ-β频段4.2 全频段耦合分析自动化分析不同频段组合from tqdm import tqdm def cfc_matrix(data, fs, phase_range(4,8), amp_range(30,80), n_bins18): 计算频段间耦合矩阵 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, pf1, fs) phase np.angle(hilbert(phase_sig)) for j, af in enumerate(amp_freqs): # 提取振幅信号 amp_sig butter_bandpass_filter(data, af-5, af5, 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, aspectauto, extent[30,80,8,4], cmapjet) plt.colorbar(labelModulation 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(nopythonTrue) def fast_mi(phase, amp, n_bins): # 实现略... return mi # 并行化处理多个通道 from joblib import Parallel, delayed def parallel_cfc(data): return Parallel(n_jobs4)(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_size1.0, step0.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:startwindow_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(picksi) data_j raw.get_data(picksj) 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)