简介这份资源围绕格兰杰因果框架下的部分定向相干PDC方法展开面向从事脑电、脑电肌电信号分析与神经科学研究的读者用于定量评估多变量时间序列之间的定向影响弥补普通相关性分析难以揭示因果方向的不足。压缩包共13个文件以10个m脚本为主辅以2个txt说明文档和1个mat数据文件整体约80KB涵盖AR模型系数估计、协方差计算、模型定阶、PDC与DTF矩阵计算以及短时连接分析等模块并附带仿真模型与示例脑电数据便于直接运行验证。目前已有723人学习下载。读者可据此搭建从数据预处理、模型拟合到连接指标计算的完整流程理解PDC值在0到1区间内如何反映脑区之间或大脑与肌肉之间的影响强度并借助示例脚本快速复现实验、排查参数设置问题为认知任务与神经疾病相关的网络连接研究提供可复用的分析工具。1. 从脑电到肌电为什么格兰杰-部分定向相干法值得你花时间做脑电肌电联合分析的人迟早会撞上一个问题感觉运动皮层的脑电信号和肌肉的肌电信号之间到底是谁在驱动谁相关性分析能告诉你两者同步但同步不等于因果可能是共同输入也可能是肌肉反馈回皮层。格兰杰因果和部分定向相干PDC就是用来回答方向性问题的工具。PDC 在频域上把多通道信号之间的直接因果影响拆出来排除其他通道的中介效应这对脑电肌电这种多节点耦合场景特别关键。你如果正在做运动控制、康复评估或神经反馈这套方法能让你从“两者有关系”推进到“皮层在 beta 频段驱动肌肉”这种可解释的结论。我下面会从数学定义讲到代码实现再到参数怎么调、坑在哪全部按能复现的标准来写。2. 格兰杰因果与 PDC 的数学骨架从时域到频域的推导2.1 格兰杰因果的定义与直觉格兰杰因果的核心思想很朴素如果加入 X 的历史信息能显著提升对 Y 当前值的预测精度超过只用 Y 自身历史预测的效果就说 X 是 Y 的格兰杰原因。数学上对两个平稳时间序列 X(t) 和 Y(t)先建立自回归模型Y(t) Σ a_i · Y(t-i) ε1(t)再建立联合回归模型Y(t) Σ a_i · Y(t-i) Σ b_i · X(t-i) ε2(t)如果 ε2 的方差显著小于 ε1 的方差就认为 X 对 Y 有格兰杰因果影响。这里的关键前提是序列必须平稳否则 F 检验的分布假设不成立结果就是玄学。对于脑电肌电场景X 通常是某通道脑电如 C3/C4Y 是对侧肌肉的肌电包络。采样率一般在 1000 Hz 左右分析前需要降采样到 256 Hz 或 512 Hz 以减少计算量同时保留 beta 频段13-30 Hz的信息。2.2 从时域格兰杰到频域 PDC 的转换时域格兰杰因果只能告诉你“有没有”不能告诉你“在哪个频段”。PDC 的价值就在这里。对 P 阶多元自回归模型MVARX(t) Σ_{k1}^{P} A(k) · X(t-k) E(t)其中 X(t) 是 N 维信号向量A(k) 是 N×N 的系数矩阵E(t) 是白噪声。做傅里叶变换后得到频域表示A(f) I - Σ_{k1}^{P} A(k) · exp(-j2πfk)PDC 的定义为PDC_{ij}(f) |A_{ij}(f)| / sqrt(Σ_{n1}^{N} |A_{nj}(f)|²)这个公式的含义是从 j 到 i 的直接因果影响除以所有指向 i 的因果影响之和。分母的归一化保证了 PDC 值在 0 到 1 之间而且排除了间接路径的影响。这就是“部分”二字的来源——它只保留直接连接把通过其他通道中转的间接影响剔除掉了。2.3 脑电肌电场景下的特殊考量脑电和肌电的信号特性差异很大。脑电幅度在微伏级肌电在毫伏级直接做 MVAR 会出问题。常见做法是先对肌电做整流和低通滤波得到包络再对两路信号分别做 z-score 标准化。另外脑电肌电之间的传导延迟通常在 20-50 ms 量级这意味着在 MVAR 模型中阶数 P 的选择要能覆盖这个延迟。以 256 Hz 采样率算50 ms 对应约 13 个采样点所以 P 至少取 15 以上才稳妥。还有一个容易翻车的地方容积传导。脑电电极可能直接拾取到肌肉电活动造成虚假的因果方向。解决办法是用 Laplacian 蒙太奇或独立成分分析先做源分离把肌肉伪迹从脑电里去掉。3. 用 Python 跑通 PDC 计算从数据预处理到频域因果图3.1 环境准备与依赖安装我一般用 Python 做这套分析核心依赖是 numpy、scipy 和 statsmodels。statsmodels 里的 VAR 模块可以直接拟合 MVAR 模型省去手写最小二乘的麻烦。pip install numpy scipy statsmodels matplotlib版本方面numpy 1.24 以上、scipy 1.10 以上、statsmodels 0.14 以上都能跑。不建议用太老的版本statsmodels 在 0.12 之前 VAR 模块的稳定性有问题。3.2 数据预处理滤波、降采样与标准化假设你已经有两路信号脑电通道 eeg单位微伏和肌电包络 emg单位毫伏采样率 1000 Hz。import numpy as np from scipy.signal import butter, filtfilt, hilbert, resample def preprocess(eeg, emg, fs_orig1000, fs_target256): # 带通滤波 1-45 Hz去掉漂移和高频噪声 b, a butter(4, [1/(fs_orig/2), 45/(fs_orig/2)], btypeband) eeg_f filtfilt(b, a, eeg) emg_f filtfilt(b, a, emg) # 肌电取包络希尔伯特变换后取模再低通到 10 Hz emg_env np.abs(hilbert(emg_f)) b2, a2 butter(4, 10/(fs_orig/2), btypelow) emg_env filtfilt(b2, a2, emg_env) # 降采样到 256 Hz n_target int(len(eeg_f) * fs_target / fs_orig) eeg_ds resample(eeg_f, n_target) emg_ds resample(emg_env, n_target) # z-score 标准化消除量纲差异 eeg_ds (eeg_ds - np.mean(eeg_ds)) / np.std(eeg_ds) emg_ds (emg_ds - np.mean(emg_ds)) / np.std(emg_ds) return eeg_ds, emg_ds, fs_target滤波用 filtfilt 而不是 lfilter因为 filtfilt 是零相位滤波不会引入时间延迟。肌电包络的低通截止频率选 10 Hz 是常见做法能保留肌电的激活轮廓。降采样前一定要先滤波否则会混叠。标准化这一步不能省脑电和肌电的幅度差三个数量级不标准化的话 MVAR 系数会被大量级信号主导。3.3 拟合 MVAR 模型并计算 PDCfrom statsmodels.tsa.api import VAR def compute_pdc(eeg, emg, fs, order20, freq_range(1, 45)): # 构建多通道数据矩阵形状为 (时间点, 通道数) data np.column_stack([eeg, emg]) # 拟合 VAR 模型 model VAR(data) result model.fit(order) # 提取系数矩阵 A(k)k1...order coefs result.coefs # 形状 (order, 2, 2) # 计算频率向量 freqs np.linspace(freq_range[0], freq_range[1], 100) pdc_12 np.zeros(len(freqs)) # 脑电 - 肌电 pdc_21 np.zeros(len(freqs)) # 肌电 - 脑电 for idx, f in enumerate(freqs): # 计算 A(f) I - sum(A(k) * exp(-j2πfk/fs)) A_f np.eye(2, dtypecomplex) for k in range(order): A_f - coefs[k] * np.exp(-1j * 2 * np.pi * f * (k1) / fs) # PDC 归一化 denom np.sqrt(np.sum(np.abs(A_f)**2, axis0)) pdc_12[idx] np.abs(A_f[0, 1]) / denom[1] pdc_21[idx] np.abs(A_f[1, 0]) / denom[0] return freqs, pdc_12, pdc_21order 参数就是 MVAR 模型的阶数前面说了脑电肌电场景建议 15-25我一般先用 20 试。freq_range 覆盖 1-45 Hzbeta 频段在 13-30 Hz 之间看结果时重点盯这个区间。A_f 的计算里指数项的符号是负的因为傅里叶变换的定义是 exp(-j2πft)这里 t 用采样点索引除以采样率来近似。3.4 显著性检验打乱试验的代理数据方法PDC 值本身没有统计显著性需要做代理数据检验。常用方法是随机打乱试验间的时间对齐关系重复计算 PDC得到零分布。def surrogate_test(eeg_trials, emg_trials, fs, order20, n_surrogate200): # eeg_trials 和 emg_trials 形状为 (试验数, 时间点) n_trials, n_points eeg_trials.shape _, pdc_real, _ compute_pdc(eeg_trials.flatten(), emg_trials.flatten(), fs, order) pdc_surr np.zeros((n_surrogate, len(pdc_real))) for s in range(n_surrogate): # 随机打乱试验顺序破坏脑电肌电的跨试验对应关系 perm np.random.permutation(n_trials) eeg_shuffled eeg_trials[perm].flatten() _, pdc_s, _ compute_pdc(eeg_shuffled, emg_trials.flatten(), fs, order) pdc_surr[s] pdc_s # 计算每个频率点的 p 值 p_values np.mean(pdc_surr pdc_real, axis0) return pdc_real, p_valuesn_surrogate 取 200 是精度和耗时的折中条件允许可以上 1000。打乱的是试验顺序而不是时间点因为时间点打乱会破坏信号的自相关结构导致零分布偏移。p 值小于 0.05 的频率点才认为因果连接显著。4. 参数调优与结果解读阶数、频段和方向的判断标准4.1 MVAR 阶数怎么选AIC、BIC 与经验值的权衡阶数 P 是 PDC 分析里最关键的参数。太小模型欠拟合因果方向可能被漏掉太大模型过拟合噪声被当成因果。statsmodels 的 VAR 模块提供了 select_order 方法基于 AIC 或 BIC 自动选阶。model VAR(data) order_result model.select_order(maxlags40) print(order_result.summary())AIC 倾向于选大一点的阶数BIC 倾向于选小的。我的经验是如果 AIC 和 BIC 给出的阶数差在 5 以内取两者的中间值如果差很多优先信 BIC因为脑电肌电信号的信噪比通常不高过拟合的风险更大。另外选出来的阶数要跟生理延迟对照一下——如果 BIC 选出来 P5对应 256 Hz 下约 20 ms刚好在传导延迟范围内那可以接受如果选出来 P2对应 8 ms明显短于已知的皮层到肌肉传导时间那就得手动加大到 15 以上。4.2 频段聚焦beta 和 gamma 频段的因果模式差异脑电肌电相干性最经典的发现是在 beta 频段13-30 Hz感觉运动皮层和肌肉之间在这个频段有显著的因果耦合。PDC 分析通常也会在这个频段看到脑电到肌电方向的显著峰值。gamma 频段30-45 Hz有时也能看到但更容易受肌电伪迹污染解读时要谨慎。实际操作中我会把 PDC 谱分成几个频段分别做统计delta1-4 Hz、theta4-8 Hz、alpha8-13 Hz、beta13-30 Hz、gamma30-45 Hz。每个频段取 PDC 的均值再做组间比较。这样比逐频率点比较更稳健也更容易跟临床指标对应。4.3 方向性判断脑电到肌电 vs 肌电到脑电PDC 的两个方向要同时看。如果脑电到肌电的 PDC 在 beta 频段显著高于肌电到脑电说明皮层驱动占主导这是运动执行时的典型模式。如果反过来肌电到脑电更强可能反映的是感觉反馈通路在运动想象或被动运动时更常见。有个容易搞混的地方PDC 的归一化是在“指向目标通道的所有连接”上做的所以脑电到肌电的 PDC 和肌电到脑电的 PDC 不是直接可比的它们各自归一化在不同的分母上。要比较方向强弱应该看原始 A(f) 矩阵的对应元素模值或者用未归一化的部分定向相干。这一点很多论文里都没写清楚审稿人有时候也会揪。5. 避坑与排查PDC 分析中最容易翻车的五个地方5.1 信号不平稳导致因果方向随机翻转现象同一组数据分两段做 PDC脑电到肌电的峰值频率和方向都不一致甚至完全相反。原因格兰杰因果和 PDC 都建立在信号平稳的假设上。脑电肌电信号在运动开始、持续、结束各阶段的统计特性变化很大如果直接把整段数据扔进去MVAR 模型拟合出来的系数是各阶段平均的结果没有物理意义。解决做 PDC 之前先做平稳性检验。常用的是 ADF 检验和 KPSS 检验两者结论一致才认为平稳。如果不平稳做一阶差分或者分段处理。我一般按运动事件把数据切成 500 ms 的窗口每个窗口单独做 PDC再看时间演化。5.2 容积传导造成虚假的脑电到肌电因果现象脑电到肌电的 PDC 在所有频段都显著没有频率选择性而且脑电通道越靠近肌肉越明显。原因脑电电极拾取到了肌肉电活动这部分信号跟肌电是同一个源当然会有“因果”。但这是伪迹不是神经通路。解决用 Laplacian 蒙太奇做空间滤波或者用 ICA 把肌肉成分从脑电里分离出去。判断标准是真正的皮层肌肉耦合有频率选择性beta 频段峰值容积传导没有。如果 PDC 谱是平的基本可以断定是伪迹。5.3 阶数选择不当导致因果方向反转现象P5 时脑电到肌电显著P25 时变成肌电到脑电显著。原因阶数太低时模型无法捕捉到真实的延迟把间接路径当成了直接路径阶数太高时过拟合的噪声可能掩盖真实连接。方向反转说明模型不稳定。解决用 AIC/BIC 选阶后在选出的阶数附近做敏感性分析比如选 P15、20、25 各跑一遍看方向是否一致。如果方向随阶数剧烈变化说明数据质量不够需要增加试验次数或改善信噪比。5.4 试验次数不足导致代理检验不显著现象PDC 谱看起来有峰值但代理检验 p 值都在 0.1 以上怎么都到不了 0.05。原因代理检验的零分布估计需要足够的试验次数。如果只有 20 次试验打乱后的零分布很粗糙p 值的最小分辨率就是 1/2000.005但实际功效很低。解决试验次数至少 50 次以上最好 100 次。如果试验次数不够可以用 bootstrap 方法替代 permutation或者降低显著性阈值到 0.01 但增加 n_surrogate 到 1000。另一个办法是把多个被试的数据合并做混合效应模型但这要求被试间条件一致。5.5 肌电包络参数影响 beta 频段因果强度现象肌电包络低通截止频率从 10 Hz 改到 5 Hzbeta 频段的 PDC 峰值掉了 30%。原因肌电包络的截止频率决定了包络里保留了多少 beta 频段信息。截止频率太低beta 振荡被滤掉因果强度自然下降。解决肌电包络的低通截止频率不要低于 10 Hz如果重点分析 beta 频段可以放宽到 20 Hz。但也不能太高否则包络里混入高频噪声。我一般先用 10 Hz如果 beta 频段 PDC 不显著再试 20 Hz 对比。6. 进阶技巧用滑动窗口 PDC 追踪因果连接的时变特性静态 PDC 给出的是整段数据的平均因果强度但脑电肌电耦合在运动过程中是动态变化的。滑动窗口 PDC 能追踪这种时变特性代价是每个窗口的数据点变少MVAR 拟合的方差变大。窗口长度的选择是个权衡太短模型不稳定太长时间分辨率不够。我的经验是窗口至少包含 10 倍于阶数的数据点。如果阶数 P20窗口至少 200 个采样点256 Hz 下约 780 ms。def sliding_pdc(eeg, emg, fs, order20, win_len256, step64): n_points len(eeg) n_wins (n_points - win_len) // step 1 pdc_beta np.zeros(n_wins) time_axis np.zeros(n_wins) for w in range(n_wins): start w * step end start win_len eeg_win eeg[start:end] emg_win emg[start:end] freqs, pdc_12, _ compute_pdc(eeg_win, emg_win, fs, order) # 取 beta 频段 13-30 Hz 的均值 beta_idx (freqs 13) (freqs 30) pdc_beta[w] np.mean(pdc_12[beta_idx]) time_axis[w] (start win_len/2) / fs return time_axis, pdc_betawin_len 取 256 对应 1 秒窗口step 取 64 对应 75% 重叠。重叠是为了让时间曲线更平滑但会引入自相关做统计检验时要注意校正。如果发现 beta 频段 PDC 在运动开始后 200-300 ms 出现峰值跟皮层肌肉传导延迟对得上那这个结果就比较可信。如果峰值出现在运动开始前要么是预激活要么是伪迹得回去查数据质量。我做了这么多年脑电肌电分析最大的教训就是PDC 的结果好看不等于正确。每次跑出显著因果我都会问自己三个问题——信号平稳吗容积传导排除了吗阶数敏感性查了吗这三个问题有一个答不上来结果就不敢往论文里写。希望帮到你。本文还有配套的精品资源点击获取