阶段五 · 前沿专题 · 模块 8 · 第 25 课 / 共 28 课

CFC 跨频率耦合（下）：Matlab 实现

把上一课的原理落成代码：滤波、希尔伯特、相位分箱、调制指数、置换检验， 一环扣一环；再用仿真数据验证——先让代码抓住你亲手埋进去的耦合。

🕐 60 分钟🎯 难度：高阶✅ 前置：8.3

🎯 本节学习目标

  • 独立实现完整 PAC 管线：滤波 → 希尔伯特 → 分箱 → MI → 置换检验
  • 说清每一步"为什么这么做"：零相位、掐边缘、循环平移而非随机打乱
  • 养成仿真验证的习惯：构造已知耦合的信号检验管线能否还原它
  • 会画并解读 comodulogram 调制图
  • 掌握组级分析骨架与 PAC 论文的报告规范

一、本课地图：从一段信号到一张调制图

① 选定通道 / ROI 信号（连续数据，非分段）
② 零相位带通滤波：低频带（相位来源）+ 高频带（幅值来源），掐掉两端边缘
③ 希尔伯特变换：低频相位 ph、高频包络 am
④ 相位分箱 → 调制指数 MI（观测值）
⑤ 循环平移置换检验 → 零分布 → p 值
⑥ 仿真验证 → comodulogram → 组级统计

类比一下：这套流程就像考驾照——先在练习场（仿真数据）上确认车没问题， 再上真实道路（真实 EEG）。本课按顺序把每一环写出来， 每段代码都先回答"为什么"。

二、第一步：零相位带通滤波

相位是本课的主角，所以滤波方式是安全底线： 必须零相位。普通 filter 会给不同频率引入不同的相移， 直接把低频相位整体挪歪——相位分析当场作废。filtfilt 正反两次滤波，幅值特性平方、相位恰好抵消，是相位类分析的标准选择。

Matlab
% ---- 步骤 1：零相位带通滤波（需要 Signal Processing Toolbox）----
fs = 500;
x  = eeg_data(:);                         % 单通道（或 ROI 平均）连续信号
[b1, a1] = butter(2, [4 8]    / (fs/2));  % 低频带 theta：相位的来源
[b2, a2] = butter(2, [60 100] / (fs/2));  % 高频带 gamma：幅值的来源
lo = filtfilt(b1, a1, x);                 % filtfilt：正反两次 → 零相位
hi = filtfilt(b2, a2, x);                 % 若用 filter()，相移会污染相位结果

% 边缘效应：滤波在数据两端产生不可靠的过渡——分析前掐头去尾
cut = 2 * fs;                             % 保守起见掐掉两端各 2 秒
lo = lo(cut+1 : end-cut);
hi = hi(cut+1 : end-cut);
⚠️ 两个高频报错点

① 阶数别贪高：butter(2, [f1 f2]) 已是四阶带通，再陡会振铃；② 边缘必须掐：哪怕只掐 1 秒，也比不掐强——上一课的陷阱清单第 1 条在这里落地。

三、第二步：希尔伯特变换取相位与包络

Matlab
% ---- 步骤 2：希尔伯特变换——相位与包络 ----
ph = angle(hilbert(lo));                  % 低频瞬时相位 ∈ (-pi, pi]
am = abs(hilbert(hi));                    % 高频瞬时幅值（包络）
% 直觉：hilbert 把实信号补成复解析信号。angle 回答"此刻在圆周哪个角度"，
% abs 回答"此刻绕圈的半径多大"。两者都只对窄带信号有意义——所以必须先滤波

% 永远先目检这两个序列（取前 3 秒）：
figure;
subplot(2,1,1); plot(lo(1:3*fs)); hold on; plot(ph(1:3*fs));
legend('低频信号', '低频相位'); title('相位应随信号起伏而循环转动');
subplot(2,1,2); plot(hi(1:3*fs)); hold on; plot(am(1:3*fs));
legend('高频信号', '高频包络'); title('包络应贴着高频的峰顶走');

四、第三步：相位分箱与调制指数

把上一课的直觉代码封装成函数——科研代码的第一步就是把"能复用的东西"变成函数， 后面的置换检验和 comodulogram 都要反复调用它。

Matlab · modulation_index.m
function mi = modulation_index(ph, am, nBins)
% 调制指数 MI：高频幅值在低频相位上的分布偏离均匀的程度
% 直觉：把每箱平均幅值画成柱状图，MI 度量"这排柱子有多不像一条平线"
% 保存为 modulation_index.m（函数文件名 = 函数名）
bIdx = discretize(ph, linspace(-pi, pi, nBins + 1));  % 逐时刻分箱
ampBin = accumarray(bIdx, am, [nBins 1], @mean);      % 每箱平均幅值
p = ampBin / sum(ampBin);                             % 归一化成概率分布
mi = sum(p .* log(max(p, eps) * nBins));              % 与均匀分布的 KL 散度
end

% 调用：得到观测值
miObs = modulation_index(ph, am, 18);
fprintf('观测 MI = %.4f\n', miObs);

五、第四步：置换检验——给 MI 一个零分布

MI = 0.03 算大算小？没有参照系就无法回答。置换检验造一个参照系： 反复破坏相位与幅值的时间对应关系，重算 MI，得到"没有耦合时 MI 长什么样"的分布。 这里用循环平移而不是随机打乱——平移保留了信号自身的频谱与自相关结构， 得到的零分布更保守、更合理。

Matlab
% ---- 步骤 4：置换检验（循环平移包络）----
nPerm  = 1000;
miNull = zeros(nPerm, 1);
for p = 1:nPerm
    shift = randi(numel(am) - 1);
    miNull(p) = modulation_index(ph, circshift(am, shift), 18);
end
pval = (1 + sum(miNull >= miObs)) / (1 + nPerm);    % +1 避免 p = 0
zval = (miObs - mean(miNull)) / std(miNull);
fprintf('MI = %.4f, z = %.2f, p = %.4f\n', miObs, zval, pval);

figure; histogram(miNull, 40); hold on;
xline(miObs, 'r', 'LineWidth', 2);                  % 观测值应远在零分布右尾
xlabel('零分布 MI'); title(sprintf('置换检验 p = %.3f', pval));

六、先仿真，再真实数据

这一步是科研必备习惯：构造一个已知耦合的信号，看管线能不能把它挖出来。 连自己亲手埋的耦合都找不到，说明代码有 bug；挖出来的位置或强度不对， 说明参数有问题。修好之前，别碰真实数据。

Matlab
% ---- 仿真验证：把已知的耦合埋进信号，看管线能否还原 ----
fs = 1000; T = 60; t = (0:1/fs:T)';
fP = 6; fA = 80;                                  % 埋一个 (6Hz 相位, 80Hz 幅值) 耦合
theta = sin(2*pi*fP*t);                           % 低频载波
env   = 1 + 0.8 * sin(2*pi*fP*t - pi/2);          % 包络被 6Hz 调制，滞后四分之一周期
gamma = env .* sin(2*pi*fA*t);                    % 幅值被调制的高频载波
xSim  = theta + 2*gamma + 0.5*randn(size(t));     % 加噪声，模拟真实信噪比

% 用同一套管线处理 xSim（频带相应改为 5-7Hz 与 70-90Hz）
[b1, a1] = butter(2, [5 7]   / (fs/2));
[b2, a2] = butter(2, [70 90] / (fs/2));
loS = filtfilt(b1, a1, xSim);
hiS = filtfilt(b2, a2, xSim);
miSim = modulation_index(angle(hilbert(loS)), abs(hilbert(hiS)), 18);
fprintf('仿真数据 MI = %.4f（应明显高于噪声水平）\n', miSim);
% 再对纯噪声跑一遍同样管线，得到"无耦合基线"，两者应拉开明显差距
💡 仿真三连问

每次写完新管线问自己三个问题：① 已知耦合能被找到吗（灵敏度）？② 纯噪声会被误报吗（特异性）？③ 参数微调时结果稳定吗（稳健性）？三个都答"是"，才有资格上真实数据。

七、可视化：comodulogram

📖 概念：comodulogram（调制图）

把所有"低频相位频率 × 高频幅值频率"组合的 MI 画成一张热图。亮斑所在的位置，就是数据里存在的耦合频对——它是 PAC 结果的"全家福"。

Matlab
% ---- comodulogram：扫描所有频对的 MI ----
% modulation_index 与 bp_filt 保存为独立函数文件（或放在脚本末尾）
fPhList = 2:1:12;                     % 低频相位候选
fAmList = 40:10:150;                  % 高频幅值候选
M = zeros(numel(fPhList), numel(fAmList));
for i = 1:numel(fPhList)
    for j = 1:numel(fAmList)
        loB = bp_filt(x, fPhList(i) + [-1 1], fs);     % 低频 ±1Hz 窄带
        hiB = bp_filt(x, fAmList(j) + [0 20], fs);     % 高频 20Hz 宽带
        M(i, j) = modulation_index(angle(hilbert(loB)), ...
                                   abs(hilbert(hiB)), 18);
    end
end
figure; imagesc(fAmList, fPhList, M); axis xy; colorbar;
xlabel('高频幅值频率 (Hz)'); ylabel('低频相位频率 (Hz)');
title('Comodulogram：每个点是一个频对的 MI');

function y = bp_filt(x, fRange, fs)
% 教学用零相位带通；正式分析建议统一改用 FieldTrip / EEGLAB 的滤波器
[b, a] = butter(2, fRange / (fs/2));
y = filtfilt(b, a, x);
end

八、组级分析骨架与报告规范

组级的思路很朴素：给每个被试算一张 comodulogram（或先验频对上的单个 MI）， 然后逐频对做统计。骨架如下：

  1. 循环被试 → 每人得到一个 MI 矩阵 Msub（频对 × 或通道 ×）；
  2. 逐频对做单样本检验（MI 是否大于置换零分布）或组间 t / 置换检验；
  3. 全部 p 值做 FDR 或 cluster 校正（模块 5 的方法原样搬过来）；
  4. 报告先验频对的效应量，comodulogram 作为补充图。
📌 PAC 论文报告清单

滤波器类型与阶数、是否零相位、边缘剔除长度、低频/高频带定义、相位箱数、MI 的具体定义、置换次数与打乱方式（循环平移）、多重比较校正方法、仿真验证结果、低频波形正弦性检查。这十项写全，方法部分就立于不败之地。

📌 本节小结

  • 完整 PAC 管线六步：滤波 → 希尔伯特 → 分箱 → MI → 置换检验 → 仿真与可视化
  • 相位类分析的底线：filtfilt 零相位 + 掐掉两端边缘
  • 置换用循环平移而非随机打乱：保留信号自身的频谱与自相关结构
  • 仿真验证是上真实数据前的必经关卡：灵敏度、特异性、稳健性三连问
  • comodulogram 是 PAC 的全家福；组级统计沿用模块 5 的置换 + FDR 框架

✏️ 课后练习

  1. 对一段静息态数据（至少 60 秒）跑完整管线，报告先验频对的 MI、z 值与置换 p 值。
  2. 把仿真代码中的调制深度从 0.8 改成 0.2 与 0.05，观察 MI 与 z 值如何衰减，体会信噪比的影响。
  3. 对同一数据画 comodulogram，标记出最亮的频对，并检查该低频的波形是否接近正弦。