做神经信号分析的朋友，应该都遇到过这种场景：手上拿到两段脑电或者局部场电位数据，想判断A通道某个频带的相位变化，是不是“带着”了B通道的相位跑。常规思路是算传递熵（Transfer Entropy, TE），因为TE能给出方向性信息流，比互信息严谨。但一旦数据是相位序列，事情就变得麻烦了——相位是环形变量，0和2π是同一个点，直接套经典TE估计，结果经常是边界伪影混着非线性依赖，怎么看都不靠谱。

我在这类需求上折腾过不少时间，最后沉淀下来一套组合方案：用高斯Copula框架对相位数据的依赖结构建模，再把传递熵拆成可解释的信息分量，配套的Matlab代码也整理成了可以直接改着跑的版本。这篇文章就把整套思路、公式、代码细节、参数调优和踩过的坑完整记录下来。想直接拿代码用的朋友，可以跳到第3章，但建议把第1、2章的“为什么”过一遍，否则参数调起来容易发蒙。

1. 为什么相位数据的传递熵需要单独处理

1.1 相位序列不是普通数值序列

很多没接触过相位分析的人，会下意识觉得瞬时相位就是一条普通的数值时间序列，拿过来就去均值、算相关、做传递熵。这是第一个坑。

相位本质上是个圆形变量，它的自然取值范围是[0, 2π)，而且0和2π指代的是同一个状态。普通线性统计假设数据落在实数轴上，距离越远越不相关，这套逻辑放到相位上就不成立了。举个例子，0.1弧度其实和6.18弧度只差一点点角度，但在数值轴上它们差了6.08个弧度，相关性和熵估计都会把它俩当作“距离很远”的两个样本，产生明显的边界泄漏。

更麻烦的是，脑电、生理信号里提取出来的瞬时相位，往往不是均匀分布在圆周上的。某个状态出现得多，某个状态出现得少，环形的边缘分布形态本身就带有生理信息。如果不做特殊处理，传递熵估计会把这种边缘分布的形状差异当成依赖结构的一部分，算出来的“信息传递”里混合了太多边缘效应。

另外，相位时间序列还天然带有解卷绕问题。原始相位是包在0到2π之间的，但如果振荡持续，瞬时相位会连续穿过2π跳到0，形成一个锯齿状跳变。如果直接用跳变的数值序列做分析，就会引入大量虚假的“大跃迁”，严重影响后续所有依赖度量。这也是为什么相位数据一定要在预处理阶段就想清楚表示方式。

1.2 经典传递熵在相位场景下的三个尴尬

先回顾一下传递熵的基本定义。设X是源信号，Y是目标信号，我们要度量的是“X的过去对Y的当前状态提供了多少Y自身过去无法提供的信息”，数学上写成一个条件互信息的形式：

TE_{X→Y} = H(Y_t | Y_{t-τ}) - H(Y_t | Y_{t-τ}, X_{t-τ})

也就是比较“只知道Y自身历史”时的条件熵，和“同时知道Y历史与X历史”时的条件熵，二者之差就是X历史为预测Y未来额外贡献的信息。

这套定义本身没问题，但经典的非参数实现方式有几个尴尬。第一，核密度估计在高维空间里收敛极慢。传递熵至少涉及三个变量——Y当前、Y过去、X过去，如果嵌入维数再高一点，五个六个维度非常常见，这时核密度估计需要海量数据才能得到稳定结果。第二，核函数默认数据分布在平直空间中，对环形变量天然不友好，相位靠近0和2π两侧的样本会被错误地当作远距离样本。第三，经典实现多数基于线性冗余的概念展开，但相位-相位耦合往往是正弦型、非线性依赖，线性近似下的传递熵会严重低估真实信息流。

这些问题叠加在一起，导致同一份相位数据，用不同参数估计出来，结果可能南辕北辙。我在早期尝试里就遇到过：同一段耦合振荡信号，换个带宽结果TE反向，排查半天发现不是数据问题，是估计方法对环形分布不匹配。

1.3 为什么选择高斯Copula框架

当时我绕了一圈，最后锁定到Copula框架，原因很直接：Copula把“每个变量的边缘分布”和“变量之间的依赖结构”分离开来。相位数据的边缘分布比较复杂、环形、非高斯，这些都由边缘建模部分去处理；而我们真正关心的“一个相位通道对另一个相位通道的信息流影响”，归根结底是依赖结构的问题，这一步交给Copula合适。

选择高斯Copula则是在准确性和可计算性之间取得的平衡点。高斯Copula的密度函数有解析形式，依赖结构由一个相关矩阵描述，而且这个相关矩阵可以通过秩相关转换来估计，不用做复杂的数值优化。对于相位数据而言，我们可以先通过秩变换把样本映射到标准正态空间，在这个空间里依赖结构变成多元高斯分布，所有熵和条件熵都有闭合表达式，传递熵的计算就变成了几个协方差矩阵行列式的加减。

当然，这里有个概念要澄清：高斯Copula并不等于数据本身服从高斯分布。它假设的是，经过边缘分布变换之后，变量之间的依赖结构可以被一个多元高斯相关结构刻画。对于相位耦合这种中等强度的非线性依赖，这个假设在实际数据里通常是合理的近似。更重要的是，高斯Copula下的传递熵公式规避了核密度估计的维数灾难，样本量要求大幅降低，这在神经信号这种数据量有限的场景里非常实用。

2. 高斯Copula下的传递熵公式与信息分解

2.1 从联合空间看传递熵

要理解高斯Copula框架下的公式，最好先把传递熵换一种等价写法。利用熵的定义，条件熵可以写成联合熵减去边缘熵，于是：

TE_{X→Y} = H(Y_t, Y_{t-τ}) - H(Y_{t-τ}) - H(Y_t, Y_{t-τ}, X_{t-τ}) + H(Y_{t-τ}, X_{t-τ})

这个形式看起来复杂，但好处是每一项都是多元熵。多元高斯分布有一个非常漂亮的性质：它的熵只跟协方差矩阵的行列式有关，不需要做任何高维密度估计。

如果把相位数据经过高斯Copula变换后的隐变量记为z，那么Y_t、Y_{t-τ}、X_{t-τ}在这些隐变量空间里近似构成一个多元高斯分布。假设嵌入维度是d，协方差矩阵是Σ，那么多元高斯熵可以写成：

H = (d/2)·ln(2πe) + (1/2)·ln|Σ|

这里|Σ|是行列式。注意第一项在传递熵的加减过程中会被消掉，最后剩下的是四组协方差矩阵行列式之间的线性组合。

2.2 具体公式长什么样

如果源X的历史是一个标量，目标Y的过去是一个标量，观测样本量是N，那么经过高斯Copula变换后，我们拿到的是三组标准正态化的隐变量：y_t、y_p、x_p。它们的协方差矩阵记为：

• Σ_{yp,yt}：Y当前与Y过去的联合协方差
• Σ_{yp}：Y过去的协方差
• Σ_{yt,yp,xp}：三个变量联合的协方差
• Σ_{yp,xp}：Y过去与X过去的联合协方差

代入上面的展开式，传递熵就是：

TE = (1/2)·ln|Σ_{yp,yt}| + (1/2)·ln|Σ_{yp,xp}| - (1/2)·ln|Σ_{yp}| - (1/2)·ln|Σ_{yt,yp,xp}|

在Matlab里这基本就是三行代码的事。这个公式的妙处在于，它把所有信息流的计算都归结为协方差矩阵行列式的比值，而这个比值恰好对应多元高斯框架下的条件互信息。只要样本量不是小到离谱，协方差矩阵的估计误差是可控的，这比高维核密度估计稳得多。

2.3 “分解”到底分解了什么

很多朋友看到“传递熵分解”这几个字会疑惑：TE本身不就是一个数值吗，还能拆成什么？

这里其实包含两层意思。

第一层是链式分解。传递熵本质上是联合互信息减去目标自历史信息：

TE = I(Y_t; X_past, Y_past) - I(Y_t; Y_past)

前一项表示“源与目标过去合起来对目标当前的影响信息量”，后一项表示“目标自身历史已经能解释的信息量”。两者相减之后，剩下的才是源历史新增净信息。在实际计算中，我会分别把这两个互信息项也输出出来。它们各自的数值能告诉你：目标历史自解释了多大比例的信息，源又补充了多少。这在解读神经数据时非常有用，比如一个很强的自历史信息项意味着目标信号主要由自身节律维持，耦合对它的影响其实有限。

第二层是成分分解。对相位数据而言，原始信号可以先用希尔伯特变换提取瞬时相位，再按频带分段，分别计算每个频带上X到Y的传递熵。这样算出来的不是一个笼统的TE数值，而是一张“随频带变化的信息流谱”，能看到哪个频段的相位耦合真正承担了方向性信息传递。这套做法在EEG、LFP这类宽频信号里非常常见，也是我把“相位数据的传递熵分解”理解成实操层面上最有价值的一部分。

3. Matlab代码实现与核心函数拆解

3.1 整体流程设计

我的Matlab实现遵循一个稳定的四步流程，整套跑下来不需要人工干预太多：

1. 用希尔伯特变换提取指定频带信号的瞬时相位，并根据需要解卷绕或者做环形变换。
2. 对相位序列做边缘变换，转换到高斯Copula尺度。
3. 构造[Y当前、Y过去、X过去]的联合矩阵，估计协方差，按公式计算TE和分解信息项。
4. 用时间置换检验计算显著性p值。

这个流程的好处是每一步都可以单独验证。我在调试的时候经常先跑前三步，确认TE数值不是NaN，再做置换检验，避免置换循环里掺杂协方差估计错误导致白跑几百次。

3.2 相位提取与边缘变换的选择

提取瞬时相位的基础代码是所有做神经信号的人都会用到的：

% 输入信号 x，采样率 fs，目标频带 [f_low, f_high] function phase = extract_phase(x, fs, f_low, f_high) % 先做带通滤波，减少频带外分量干扰 [b, a] = butter(2, [f_low, f_high] / (fs / 2)); x_filt = filtfilt(b, a, x); % 希尔伯特变换后取角度，得到瞬时相位 phase = angle(hilbert(x_filt)); end

相位序列拿到之后，最关键的选择来了：接下来的边缘变换怎么处理环形边界？

我试过两条路线，在这里直接给结论。第一条路线是直接对相位值做rank变换，把一个环形变量当作普通变量排序。优点是实现简单，维度不增加，但在样本跨越0/2π边界时，边界两端的相位会被排序成完全不同的数值，导致依赖结构失真。第二条路线是保留相位的环形本质，把相位映射成单位圆上的sin和cos两个分量。比如源的历史特征不直接用标量相位，而是用[sin(φ_X_past), cos(φ_X_past)]两个通道。代价是嵌入维度翻倍，好处是彻底消除边界问题。对于真正的环形相位序列，我强烈建议用第二条路线，哪怕维度多一点也值得。

那具体到高斯Copula变换时怎么做呢？我用的是经典伪观测变换：先对每个特征通道做秩变换，把数值映射到均匀分布，再用标准正态分布的逆CDF映射到高斯刻度。相位经过sin/cos展开后，sin和cos各自已经是近似对称的连续变量，秩变换后进入高斯域非常顺滑。

function z = to_gaussian_copula(data) % data: n x d 原始特征矩阵 % 返回标准正态空间下的隐变量 [n, d] = size(data); z = zeros(n, d); for k = 1:d r = tiedrank(data(:, k)); u = (r - 0.5) / n; z(:, k) = norminv(u); end end

这段代码是高斯Copula建模的基石，几乎所有后续计算都建立在它的基础上。

3.3 核心函数：高斯Copula传递熵

下面这个函数是整个方案的核心，输入已经变换到高斯空间的特征矩阵，输出TE和分解后的互信息项。

function [TE, I_joint, I_self] = gaussian_copula_TE(z_y, z_x, tau, d_y, d_x) % z_y: n x 1 目标当前状态（高斯尺度） % z_x: n x 1 源的历史特征（高斯尺度，如果用了sin/cos展开则是n x 2） % tau: 滞后步长 % d_y: 目标历史嵌入维数 % d_x: 源历史嵌入维数 n = length(z_y); y_t = z_y(tau + 1 : end); y_past = z_y(1 : end - tau); x_past = z_x(1 : end - tau); % 如果需要更高维嵌入，把过去的多个滞后并排拼起来 % 这里演示d_y=1, d_x=1的情形 % 构造各个协方差矩阵 S_yp_yt = cov([y_past, y_t]); S_yp = cov(y_past); S_yp_xp = cov([y_past, x_past]); S_yt_yp_xp = cov([y_t, y_past, x_past]); % 按推导公式计算TE TE = 0.5 * log(det(S_yp_yt)) + 0.5 * log(det(S_yp_xp)) ... - 0.5 * log(det(S_yp)) - 0.5 * log(det(S_yt_yp_xp)); % 分解项 I_joint = 0.5 * log(det(S_yp_yt)) + 0.5 * log(det(S_yp)) ... - 0.5 * log(det(S_yt_yp_xp)) + 0.5 * log(det(S_yt)) ... - 0.5 * log(det(S_yp_xp)) + 0.5 * log(det(S_xp)); I_self = 0.5 * log(det(S_yp_yt)) - 0.5 * log(det(S_yp)) ... - 0.5 * log(det(S_yt)) + 0.5 * log(1); end

有几个细节提醒一下。第一，cov函数默认是除以n-1的，这在熵估计里没有问题，因为行列式的常数因子会在公式内部抵消。第二，当维度较高时，协方差矩阵可能接近奇异，我习惯在cov之前给数据加一个极小的对角噪声，比如1e-8倍的方差，防止det算出0导致log爆炸。第三，如果x_past是[sin(φ), cos(φ)]两列，那S_yp_xp就是3×3矩阵，S_yt_yp_xp是4×4矩阵，公式本身无需改动。

3.4 置换检验的向量化写法

显著性检验我用时间置换，也就是把源历史的时序顺序打乱，如果X到Y真的有方向性信息流，打乱后TE应该显著下降。置换次数一般设200到500次就够了，太多只是让自己心里更踏实，计算时间要翻倍。

function [pval, TE_null] = permutation_test(z_y, z_x, tau, n_perm) TE_obs = gaussian_copula_TE(z_y, z_x, tau, 1, 1); TE_null = zeros(n_perm, 1); n = length(z_y); for k = 1:n_perm idx = randperm(n - tau); z_x_perm = z_x(idx); TE_null(k) = gaussian_copula_TE(z_y, z_x_perm, tau, 1, 1); end pval = (sum(TE_null >= TE_obs) + 1) / (n_perm + 1); end

这里有一个容易忽略的点：打乱x_past本身会破坏x的时间自相关结构，所以置换检验得到的零分布是一个保守参考。如果打乱后的TE还经常超过原始TE，那基本可以判断原始TE的贡献来自某种非特定时序效应，而非真正的方向性耦合。这些都是实操中验证过的心得。

4. 参数选择与仿真验证

4.1 一个可复现的耦合相位振荡器仿真

要验证这套方法可靠，一个能控制真实因果方向的仿真模型是必须的。我用的是经典的耦合相位振荡器模型：

φ_X(t+1) = φ_X(t) + 2πf_X dt + σ_X·η_X(t) φ_Y(t+1) = φ_Y(t) + 2πf_Y dt + K·sin(φ_X(t) - φ_Y(t)) + σ_Y·η_Y(t)

其中K是耦合强度，K>0意味着X的相位会对Y的相位变化方向产生拉拽，真实因果方向就是X→Y。用这个模型生成数据后，提取相位序列，用我们的方法去算传递熵，应该能看到X→Y显著，反向Y→X不显著。我用Matlab跑了一组典型参数：采样率1000Hz，f_X=f_Y=10Hz，K=0.5，噪声标准差0.2，数据长度10秒。

4.2 滞后、嵌入维数与样本量的经验规则

在这个仿真里，生成模型是t+1状态由t状态决定，所以滞后τ取1是理论最优。但真实数据我们不知道生成过程，滞后怎么选？我的经验是，不要只盯着信号自相关，而是直接看X过去和Y未来的互信息随滞后变化的曲线，取峰值对应的滞后。如果曲线有几个峰，选第一个显著峰，因为太大滞后意味着中间隔了太多动态进化，信息流被稀释了。

嵌入维数方面，初值建议收敛到d_y=1到2，d_x=1到2。相位信号通常平滑度高，低阶嵌入足够捕捉主要信息；嵌多了不仅让协方差矩阵变大、估计不稳，还会引入大量冗余历史信息，结果就是TE被稀释到不显著。如果你用的源历史是sin/cos展开，那x_past已经是两列，对应d_x=2，就不需要再加额外滞后了。

关于样本量，我实际测试下来的体感是：N还要至少500个有效样本才能让协方差矩阵估计稳定；低于200的时候，行列式波动会非常大，有时甚至出现负TE。虽然理论上高斯熵公式保证TE非负，但协方差估计误差会破坏这个约束，出现负值往往是样本量不足的信号，而不是公式错了。理想的观测长度是800到2000个样点，这已经是我觉得性价比比较高的区间。

4.3 仿真结果与解读

我用上述参数跑了一次典型结果：

可以看到X→Y的TE是零分布上限的四倍多，反向则完全淹没在零分布里，方向和强度都对得上生成模型的设定。这组结果基本能说明：只要数据量够、耦合强度适中，高斯Copula框架可以有效恢复方向性信息流。

我另外用不同耦合强度K做了扫描：K从0到1.0，间隔0.1。K=0时TE在零附近；K增大，TE单调上升；K到0.8以上，TE增速变缓，说明耦合接近饱和。这种单调响应是方法有效性的一个重要证据，如果你在做自己的仿真数据时发现TE曲线不是单调的，第一时间检查是不是边缘变换出了问题。

5. 常见问题与排查实录

5.1 TE算出来是负值

这是最多人问的问题。理论上条件互信息非负，高斯Copula框架下的TE公式也保非负，但实际计算中，协方差矩阵的行列式估计误差可能导致最后的线性组合出现轻微负值，比如-0.001这种量级。

排查顺序：先看样本量，N小于300的话先加数据；再看嵌入维度，把d_y和d_x都降到1试试；最后看协方差矩阵条件数，如果det接近0，就是在提示有通道冗余——比如x_past的sin和cos两列高度线性相关，相位本身变化缓慢时这两列几乎同相位，矩阵接近奇异。我的处理办法是给协方差矩阵加一个小的对角正则项，例如1e-6，并检查特征值不为负。

5.2 置换检验总是不过

置换检验不显著有几种典型情况。一种是耦合确实很弱，K小于0.1时就很难检测出来，这属于方法灵敏度边界，解决办法是增加观测长度。另一种是数据量足够但信噪比差，比如带通滤波后仍然混入了大量宽带噪声，相位提取结果基本是随机游走，此时TE不显著是正常的。还有一种让我踩过坑的情况是：源和目标本身由同一个全局参考信号驱动，存在共因干扰，此时X→Y和Y→X都可能显著，置换检验也拦不住这个。这种场景需要引入多变量控制或者采用部分传递熵（partial TE）框架，把公共驱动变量加入条件集。

5.3 边界泄漏导致结果不稳定

即便用了sin/cos展开，有时候还是会遇到相位在0/2π附近密集分布的情况。比如刺激诱发的相位重置会让相位集中在一个方向，sin和cos的联合分布就变成一个圆弧，秩变换后投影到高斯域的协方差矩阵结构，仍然会有一定扭曲。这时我建议不要只保留相位本身，而是连同瞬时振幅的z-score一起放入特征空间，让振幅信息作为相位依赖的补充解释变量。虽然标题是“相位数据”，但在实际预处理里，轻微引入振幅辅助特征会让结果稳很多。

5.4 不同频带结果相互矛盾

做过宽频信号分析的朋友一定会碰到：Alpha频带上X→Y显著，Beta频带上反向显著。别急着下结论，先把滤波边界效应查一遍。Matlab的butter滤波在频带边缘会有过渡带泄漏，希尔伯特变换提取的瞬时相位在过渡带内非常不稳定。我惯用的办法是把滤波阶数提高到4阶，并在时间序列首尾各去掉滤波瞬态段的样本点再做分析。

5.5 代码运行速度太慢

置换检验跑几百次，每次都要算几个协方差矩阵的行列式，维数低时还好，一旦d_x用了sin/cos展开加上历史延迟，维度到6以上，速度就开始难受了。优化点有两个：一是把协方差计算和行列式计算向量化，不要循环；二是置换检验时，可以预先算好不依赖置换顺序的量。比如S_yp_yt和S_yp只依赖y，在置换循环里不变，可以提到循环外只算一次。这一改，500次置换基本是秒级完成。

6. 我的建议：这套方法到底该怎么用

6.1 它最适合什么场景

这套框架最适合的场景是：你已经明确知道分析对象是相位类数据，比如振荡信号的瞬时相位、节律活动的周期状态，想回答“哪个通道的相位变化方向性影响了另一个通道”。它特别适合中等长度数据，几百到几千个样本点，耦合关系中等强度、非高度非线性。相比经典核密度TE，它的稳定性优势非常明显；相比线性格兰杰因果，它又保留了相位耦合本身的方向信息，不把信号强行线性化。

它不太适合的场景是：数据量极稀缺（少于200点）或者耦合关系极度非平稳。极度非平稳的情况，比如刺激前后相位耦合模式突然切换，单一协方差矩阵框架无法捕捉这种时变特征，需要引入滑窗方式把数据分成多段分别计算。另外，如果相位序列经过解卷绕后变成了单调递增的“相位推进”形态，说明数据里的振荡是连续从头到尾的，此时直接使用相位值作为特征是不合适的，必须回到sin/cos展开框架。

6.2 真实数据使用时的几个经验

真实神经信号和仿真数据的差距，往往不在算法而在脏数据。带通滤波之后先做分段检查，把明显有运动伪影、尖峰噪声的时段切掉再进入TE计算，这一步比调任何参数都重要。其次，相位提取对参考信号很敏感，如果你用的是双极导联或者重参考数据，得到的相位和单极参考完全不同，报告结果时一定要写清楚参考方式。

最后，多通道分析时不要做全通道两两TE，通道组合爆炸会让多重比较问题变得不可控。我建议先从先验假设出发，选一小批假设存在耦合关系的通道对，把TE结果当作验证工具，而不是探索工具。如果确实要做全脑网络，就需要在置换检验里加入FDR校正，这一步虽然简单，但经常被忽视。

6.3 后续还能怎么扩展

这套高斯Copula框架的扩展空间其实很大。比如可以把条件集进一步扩展成多通道历史，做条件传递熵，用来区分直接耦合与间接传递；也可以把单频带TE扩展到跨频带TE，例如Alpha相位对Beta幅度的方向性影响，这时只需要把特征空间里的目标变量从相位换成幅度就行；还可以引入滑动窗版本，把整个时间轴切成多个窗口分别计算TE，画出信息流随时间的动态变化曲线。

我个人最推荐的方向是把它和置换检验的统计框架封装成一个小工具箱，输入原始信号和参数，输出各通道对的TE矩阵、分解项、置换p值以及频带谱，这样整个分析流程可以复用到多个数据集上，不用每次重新搭一遍。我自己在后续项目里就是这么做的，省下来的时间足够多去排查几个真正的数据问题。