别再只用PLV了！用Python从零实现EEG相位同步指数(PSI)，附完整代码与避坑指南
别再只用PLV了！用Python从零实现EEG相位同步指数(PSI)，附完整代码与避坑指南
在神经科学研究中，相位同步分析是理解大脑功能连接的重要工具。虽然相位锁定值(PLV)被广泛使用，但相位同步指数(PSI)提供了更精细的时变相位关系刻画能力。本文将带你用Python从零实现PSI计算，解决实际EEG分析中的关键问题。
1. PSI与PLV的核心差异解析
相位同步指数(Phase Synchronization Index)与相位锁定值(Phase Locking Value)虽然都用于测量神经信号间的相位关系，但在计算逻辑和应用场景上存在本质区别：
PSI的独特优势：
- 时变敏感性：PSI能捕捉瞬时相位关系变化，而PLV反映的是时间窗内的统计一致性
- 抗噪能力：PSI对信号中的随机相位扰动更具鲁棒性
- 分辨率差异：PSI可以区分0°和180°相位差，而PLV会将它们视为相同
典型应用场景对比：
| 指标 | 适用场景 | 局限性 |
|---|---|---|
| PLV | 稳态相位关系分析 | 无法反映快速相位变化 |
| PSI | 动态连接研究 | 计算复杂度较高 |
在实际EEG分析中，我们常遇到这样的问题：当两个脑区表现出间歇性耦合时，PLV可能会低估真实的相位同步程度，而PSI则能更好地捕捉这种动态特性。
2. Python实现PSI的完整流程
2.1 环境准备与数据加载
首先确保安装必要的Python库：
pip install numpy scipy matplotlib mne seaborn加载EEG数据的标准做法：
import mne import numpy as np # 加载示例数据集 raw = mne.io.read_raw_fif('sample_audvis_raw.fif', preload=True) picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False) data, times = raw[picks, :] # 获取EEG数据和时间点2.2 希尔伯特变换与瞬时相位计算
获得精确的瞬时相位是PSI计算的关键步骤。常见的错误包括直接对原始信号进行希尔伯特变换而忽略带通滤波：
from scipy.signal import hilbert, butter, filtfilt def bandpass_filter(data, low, high, fs, order=4): nyq = 0.5 * fs low = low / nyq high = high / nyq b, a = butter(order, [low, high], btype='band') return filtfilt(b, a, data) def compute_instant_phase(data, fs, band): filtered = bandpass_filter(data, band[0], band[1], fs) analytic_signal = hilbert(filtered) return np.angle(analytic_signal) # 获取相位角注意：必须分频段计算相位！全频段相位值会混淆不同节律的神经活动特征。
2.3 PSI核心算法实现
基于瞬时相位差计算PSI的标准流程：
def compute_psi(phase1, phase2): """计算两个信号间的相位同步指数""" phase_diff = phase1 - phase2 # 关键步骤：相位差转换与归一化 acot_diff = np.arctan(1/np.abs(phase_diff)) # 反余切变换 psi = (acot_diff - 0.1578) / (1.5708 - 0.1578) # 归一化到[0,1] return np.clip(psi, 0, 1) # 确保值域范围常见错误修正：
- 直接使用相位差而不进行反余切变换会导致数值不稳定
- 归一化参数选择不当会使结果失真（最大值1.5708对应π/2，最小值0.1578对应2π）
3. 多通道PSI分析与可视化
3.1 全脑连接矩阵计算
高效计算所有通道对的PSI值：
def compute_all_psi(phase_data): """phase_data形状为(n_channels, n_samples)""" n_channels = phase_data.shape[0] psi_matrix = np.zeros((n_channels, n_channels)) for i in range(n_channels): for j in range(i+1, n_channels): psi = compute_psi(phase_data[i], phase_data[j]) psi_matrix[i,j] = np.mean(psi) # 取时间平均 psi_matrix[j,i] = psi_matrix[i,j] # 对称矩阵 return psi_matrix3.2 动态PSI热力图绘制
使用Seaborn绘制专业级热力图：
import seaborn as sns import matplotlib.pyplot as plt def plot_psi_heatmap(psi_matrix, ch_names): plt.figure(figsize=(12,10)) mask = np.triu(np.ones_like(psi_matrix, dtype=bool)) # 掩码上三角 sns.heatmap(psi_matrix, mask=mask, cmap='viridis', xticklabels=ch_names, yticklabels=ch_names, annot=True, fmt='.2f') plt.title('PSI Connectivity Matrix') plt.tight_layout()4. 实战技巧与性能优化
4.1 分频段处理的必要性
不同脑节律（θ、α、β等）的相位同步具有独立生理意义。推荐处理流程：
- 定义频带范围：
freq_bands = { 'theta': [4, 7], 'alpha': [8, 12], 'beta': [13, 30] }- 分频段计算PSI：
psi_results = {} for band_name, band_range in freq_bands.items(): phases = np.array([compute_instant_phase(ch_data, sfreq, band_range) for ch_data in data]) psi_results[band_name] = compute_all_psi(phases)4.2 大数据量处理技巧
当处理长时间EEG记录或多被试数据时，可采用以下优化策略：
- 内存映射技术：使用NumPy的memmap处理大型数组
psi_arr = np.memmap('temp_psi.dat', dtype='float32', mode='w+', shape=(n_channels, n_channels))- 并行计算：利用Joblib加速循环
from joblib import Parallel, delayed def parallel_psi(i, j, phase_data): return compute_psi(phase_data[i], phase_data[j]) results = Parallel(n_jobs=4)( delayed(parallel_psi)(i,j,phase_data) for i in range(n_channels) for j in range(i+1, n_channels) )4.3 结果解读的注意事项
- PSI值接近1表示强相位同步，接近0表示无同步
- 要结合生理意义解释连接模式，避免过度解读随机波动
- 建议同时计算PLV作为参照，比较两种指标的结果差异
# 示例：PSI与PLV对比分析 plt.figure(figsize=(10,4)) plt.subplot(121) plt.hist(psi_matrix.flatten(), bins=30, alpha=0.7, label='PSI') plt.subplot(122) plt.hist(plv_matrix.flatten(), bins=30, alpha=0.7, label='PLV') plt.legend()在实际项目中，我发现PSI特别适合分析任务诱发的动态连接变化。例如在认知任务中，前额叶与顶叶的θ频段PSI常在任务开始后300-500ms出现显著升高，这种瞬态连接模式用传统PLV往往难以捕捉。