脑电功能连接指标：相干与相位同步

一、从"哪里激活"到"谁和谁对话"

前五个模块回答的都是"哪里、什么时候"激活：ERP 看电位幅度（模块 3），时频看功率（模块 4）。 但大脑是分布式系统，认知加工靠脑区之间的通信完成。模块 6 换一个问题： 哪些脑区在对话？对话得多紧密？

本模块聚焦功能连接：计算简单、物理直觉清晰，是地基。方向性指标（FieldTrip 里的 cfg.method = 'granger' 一族）建立在自回归模型上，地基不牢时极易把体积传导伪迹 读成"因果方向"，本课程不展开。

二、相位：连接分析的通用语言

频域连接指标几乎都建立在相位（phase）上。把一段窄带信号写成 x(t) = A(t) · cos( φ(t) )：A(t) 是包络（这一拍敲多重）， φ(t) 是瞬时相位（这一拍敲到循环的哪个位置）。两个信号的相位差稳定， 说明它们始终"对得上拍"。

实信号没有现成的相位，标准做法是构造解析信号：给实信号 x 配上一个虚部， 变成复信号 z(t) = x(t) + i·H(x)(t)，其中 H(x) 是希尔伯特变换。 复数用模长和辐角表示——模长就是包络 A(t)，辐角就是瞬时相位 φ(t)。 Matlab 里一行 angle(hilbert(x)) 就能拿到（模块 4 讲时频时见过同一招）。

直觉类比：两个鼓手掌心相对轮流击掌。每次击掌是一个周期，相位就是手掌在 "举起—拍下—分开"循环里的位置。两人每次都同时拍下 = 相位锁定； 一个总慢半拍 = 恒定相位差；各拍各的 = 相位差随机。本节课的所有指标， 都是在给"这两位鼓手的合拍程度"打分。

% 提取瞬时相位：滤波 -> 希尔伯特 -> angle
fs = 500; t = (0:1/fs:2)';                      % 采样率 500Hz，2 秒

x1 = sin(2*pi*10*t) + 0.3*randn(size(t));       % 10Hz 振荡 + 噪声
x2 = sin(2*pi*10*t + pi/2) + 0.3*randn(size(t)); % 同频，恒定超前 pi/2

% 1) 窄带滤波（9-11Hz）：相位只在窄带内可解释
[b, a] = butter(2, [9 11]/(fs/2), 'bandpass');
x1f = filtfilt(b, a, x1);                       % filtfilt：零相位滤波，不引入额外延迟
x2f = filtfilt(b, a, x2);

% 2) 解析信号 -> 瞬时相位（弧度）
phi1 = angle(hilbert(x1f));
phi2 = angle(hilbert(x2f));

% 3) 相位差：理论上应稳定在 -pi/2 附近
dphi = phi1 - phi2;
figure;
subplot(2,1,1); plot(t, phi1, t, phi2); legend('phi_1', 'phi_2');
title('两条瞬时相位曲线');
subplot(2,1,2); plot(t, dphi); hold on; yline(-pi/2, 'r--');
ylim([-pi pi]); title('相位差（两端的滤波边缘抖动先忽略）');

三、相干 coherence：最老牌、也最容易被骗

相干在频域定义（模块 4 的谱估计知识直接用上）：对频率 f，先算两信号的互谱 Sxy(f)——一个同时含幅值与相位信息的复数——和各自的自谱 Sxx(f)、Syy(f)，然后：

Coh(f) = | Sxy(f) | / sqrt( Sxx(f) · Syy(f) )

文字解释：分子是"两信号在这个频率上共同的成分有多强"，分母用各自能量做归一化， 所以取值 0–1。1 表示该频率上两信号完全线性相关，0 表示完全独立。它是脑电连接文献的 "开场白指标"——也是最常被误读的指标，因为有两个坑。

体积传导为什么专坑"看重相似性"的指标：源信号传到两个电极几乎同时到达（介质中的传导延迟 是微秒级，脑电采样率分辨不出），于是两电极信号零时间滞后地相似。相干、互相关 这类指标无法区分"真交互"与"同一源的两份拷贝"——这就是第五节 PLI 要下刀的地方。

四、PLV：只看相位的同步强度

PLV（phase-locking value）把注意力收窄到相位：对每个时间点（或每个试次）算相位差 Δφ(t) = φx(t) − φy(t)，把它看成单位圆上的一个向量，再取所有向量的平均：

PLV = | (1/N) · Σ e^( i·Δφ(t) ) |

文字解释：相位差完全锁定时，所有向量指向同一方向，平均向量长度 = 1；相位差均匀随机时， 向量四面八方互相抵消，平均长度 = 0。PLV 天然落在 0–1，衡量的是 相位差分布的集中度——像指南针：大家指向一致，指针长；各指各的，指针缩成一个点。

• 相比相干：PLV 剔除了幅值信息，只回答"相位是否同步"，解释更干净；
• 统计前提：N 次重复观测（多试次或多段时间窗）。单试次 PLV 没有意义—— 没有"分布"就谈不上"集中度"；
• 它不抗体积传导：同一源拷贝到两个电极，相位差恒为 0，PLV 照样 = 1。 换指标换不掉这个坑，要靠下一节的 PLI。

五、PLI：只认相位差的"符号"

PLI（phase lag index）的思路极其大胆：相位差的其他信息全部扔掉，只留正负号：

PLI = | (1/N) · Σ sign( Δφ(t) ) |

文字解释：统计"φx 比 φy 超前的次数"与"落后的次数"是否失衡。恒定超前 → 全是 +1， PLI 趋近 1；真随机 → 正负各半，PLI 趋近 0；而零滞后（Δφ ≈ 0）时， 相位差在 0 附近被噪声推来推去，正号负号对半出现，PLI 同样趋近 0。换句话说： 体积传导制造的假连接恰好落在 PLI 的"盲区"里，被自动扔掉；真实交互带着毫秒级延迟 → 相位差系统性偏离 0 → 符号失衡 → PLI 升高。

代价也很直接：正负号是信息量极低的粗糙统计量。当真实相位差本来就很小（快频段、短延迟） 或试次数少时，噪声在 0 附近频繁翻转符号，PLI 估计值剧烈抖动、系统性偏低。 这正是下一节 wPLI 要修的问题。

六、wPLI：给符号配上"权重"

wPLI（weighted phase lag index）的改进一句话讲完：不让每个相位差一票一权， 而是按可信度加权。权重取交叉谱虚部的幅度 |Im(Sxy)|——相位差离 0 越远（滞后越明确）， 虚部越大，这一票分量越重；贴着 0 的相位差权重接近 0，就算被噪声翻转了符号也几乎不影响总和：

wPLI = | Σ |w(t)| · sign( w(t) ) | / Σ |w(t)|，其中 w(t) = Im( Sxy(t) )

文字解释：分子是加权后的符号净失衡，分母归一化。效果：既保留 PLI"扔掉零滞后"的优点， 又把噪声敏感度大幅压低。同族还有去偏版本：FieldTrip 的 wpli_debiased （修正少试次时的正偏差）等。新手可以直接记结论：头皮电极水平的连接分析， wPLI（或去偏版）通常比相干更稳妥。

七、四大指标对比表

八、Matlab 手算 PLV：把公式落成代码

下面这段代码不依赖任何工具箱，从模拟数据到 PLV 全流程手写。跑通它， 你对"相位同步"的理解会超过 90% 只会调函数的人。

%% 手算 PLV：filter -> hilbert -> angle -> 相位差 -> PLV
clear; clc; close all;
rng(6);                                   % 固定随机种子，结果可复现

% --- 模拟数据：60 个试次，y 相对 x 恒定滞后 pi/2（真实交互的样子） ---
fs = 500; T = 2; t = (0:1/fs:T-1/fs)';
nTrials = 60;
x = zeros(nTrials, numel(t));             % 试次 x 时间
y = zeros(nTrials, numel(t));
for k = 1:nTrials
    p0 = 2*pi*rand;                       % 每试次初始相位随机（试次之间不同步）
    x(k,:) = sin(2*pi*10*t + p0) + 0.2*randn(size(t));
    y(k,:) = sin(2*pi*10*t + p0 + pi/2) + 0.2*randn(size(t));  % 恒定相位差
end

% --- 1) 窄带滤波：相位只在窄带内可解释 ---
[b, a] = butter(2, [9 11]/(fs/2), 'bandpass');
xf = filtfilt(b, a, x, [], 2);            % 第 5 个参数 2 = 沿第 2 维（时间）滤波
yf = filtfilt(b, a, y, [], 2);

% --- 2) 解析信号 -> 瞬时相位（弧度） ---
phix = angle(hilbert(xf, [], 2));         % hilbert(x, [], dim) 沿指定维做变换
phiy = angle(hilbert(yf, [], 2));

% --- 3) 相位差 ---
dphi = phix - phiy;

% --- 4) PLV：单位圆上相位差向量的平均长度 ---
plv = squeeze(abs(mean(exp(1i*dphi), 1)));   % 沿试次维平均 -> 1x时间的 PLV 曲线

% --- 顺手算个 PLI 对比（只用符号，见第五节） ---
pli = squeeze(abs(mean(sign(dphi), 1)));

figure;
subplot(2,1,1);
plot(t, plv, 'LineWidth', 1.5); hold on; plot(t, pli);
ylim([0 1.05]); legend('PLV', 'PLI');
title('稳态区（忽略两端 0.3s 滤波边缘）两者都应接近 1');
subplot(2,1,2);
plot(dphi(1,:), '.');                     % 看一眼相位差本身：应稳定贴着 -pi/2
ylabel('\Delta\phi (rad)'); xlabel('时间 (s)');

九、从模拟数据到真实数据

把模拟数据换成你模块 2 预处理好的真实数据，只需替换中间几行：EEG.data 是 通道 × 时间 × 试次的三维矩阵，挑两个电极（如 F3 与 O1），沿第三维做相位差即可——练习里 就做这件事。单被试单指标的数值别急着写进论文：6.2 解决"这条边是真的吗"（统计与校正）， 6.3 解决"整套流程怎么工程化"（FieldTrip）。