简介modulation_index-master是一份面向神经科学和脑电信号分析的轻量级MATLAB工具包围绕相位-幅度耦合PAC问题提供计算调制指数MI的函数实现。PAC用于刻画低频神经振荡相位对高频振荡幅度的调控MI是量化这种关联强度的常用统计量尤其适合海马区学习记忆、空间导航等场景的神经信号分析。压缩包共4个文件包含MATLAB函数源码、Markdown说明文档、开源许可证和Git属性配置整体仅3KB结构紧凑计算时通常对低频信号提取相位、对高频信号提取幅度再检验二者的相关性工具已将这一流程封装为可直接调用的函数。目前已有414人浏览学习适合神经科学方向的学生、科研新手以及需要快速上手PAC计算的实验室场景。通过这份资源使用者可免去从零搭建算法框架的时间直接获得封装好的调制指数计算函数和调用说明并理解PAC分析的关键步骤把更多精力放在结果解读与统计检验上。1. modulation_index 计算 PAC一个被低估的神经耦合度量拿到一段局部场电位或脑电数据想回答“高频活动是不是被低频节律牵着走”最常见的技术路径就是计算相位-幅度耦合Phase-Amplitude Coupling, PAC而 PAC 的定量核心就是调制指数Modulation Index, MI。一个容易被忽略的事实MI 的绝对值通常很小0.010.05 的数值在短数据里完全可能是噪声产物直接拿它下结论十有八九要翻车。modulation_index-master这类命名常见于从代码仓库下载的算法包目录它把计算 PAC 的最基础函数封装成了可直接调用的模块。这篇内容会沿着“理论拆解 → 最小实现 → 参数调优 → 统计验证”的顺序把这条链路完整走一遍适合正在处理电生理数据、想搞懂 MI 到底算了什么以及准备把它作为特征做条件比较的工程师和研究者。2. 调制指数计算 PAC 的完整链路从滤波到熵2.1 PAC 的生理信号含义与调制指数的两种计算范式相位-幅度耦合描述的是一个低频振荡比如 theta 412 Hz 或 delta 14 Hz的相位状态对另一个高频振荡gamma 30150 Hz的幅度大小产生调制。在神经科学里海马 theta-gamma 耦合被认为是工作记忆编码的重要机制在皮层信号里慢波相位对 spindle 或 ripples 幅度的调制也常被用来刻画睡眠深度、麻醉深度甚至意识状态。计算 MI 有两大流派。第一派是 Tort 等人 2008 年提出的基于香农熵的方法把相位分箱后对每个相位箱里的平均幅度做归一化再计算该分布相对均匀分布的偏离程度。第二派是 Canolty 等人 2006 年提出的复数向量法把瞬时相位作为角度、瞬时幅度作为模长构建一个复合向量向量的平均模长就是 MI。两者的数学动机不同但信号预处理环节几乎完全一样。Tort 方法对幅度分布更稳健Canolty 方法对短数据更敏感。modulation_index类仓库里最常见的实现是 Tort 变体因为它需要的假设最少且天然落在一个固定区间 [0, 1] 内。2.2 计算 MI 前的信号预处理带通滤波与希尔伯特变换无论选哪种范式第一步都是把原始 LFP 拆成“低频相位信号”和“高频幅度信号”。对于目标相位频率 f_phase需要做一个带通滤波常见带宽是 ±12 Hz对于目标幅度频率 f_amp带宽通常是 ±510 Hz 或相对带宽 20%30%。滤波完成后对低频信号做希尔伯特变换求瞬时相位对高频信号做希尔伯特变换求瞬时幅度包络。import numpy as np from scipy.signal import butter, filtfilt, hilbert def extract_phase_amp(lfp, fs, phase_freq, amp_freq, phase_bw2.0, amp_bw20.0, order3): # 提取低频相位 lo_ph phase_freq - phase_bw / 2 hi_ph phase_freq phase_bw / 2 b_ph, a_ph butter(order, [lo_ph / (fs / 2), hi_ph / (fs / 2)], btypeband) phase_signal filtfilt(b_ph, a_ph, lfp) phase np.angle(hilbert(phase_signal)) # 提取高频幅度包络 lo_amp amp_freq - amp_bw / 2 hi_amp amp_freq amp_bw / 2 b_amp, a_amp butter(order, [lo_amp / (fs / 2), hi_amp / (fs / 2)], btypeband) amp_signal filtfilt(b_amp, a_amp, lfp) amp_env np.abs(hilbert(amp_signal)) return phase, amp_envfiltfilt采用零相位滤波先正向后反向消除了滤波器引入的相位延迟这对相位计算是必须的如果改用lfilter相位会整体偏移MI 值虽然可能不变但后面画相-幅图时峰值位置就是错的。滤波器阶数默认用 3 阶巴特沃斯常规做法是阶数越低越稳高阶级数在锐利谱峰附近容易产生吉布斯振铃。2.3 相位分箱、平均幅度与香农熵的数学表达相位范围 [-π, π] 被均匀分成 N 个箱常见 N18即每个箱 20°。将所有时间点的幅度按瞬时相位落入的箱号分组计算每个箱的平均幅度。设平均幅度序列为 (j)j 1…N归一化为 P(j) (j) / ∑ 。如果相位与幅度完全无关P(j) 应该接近 1/N 的均匀分布如果强耦合P(j) 会呈现一个明显的峰值。MI 的公式是MI (ln(N) ∑ P(j) * ln(P(j))) / ln(N)等价的写法是 MI (ln(N) - H(P)) / ln(N)其中 H(P) 是香农熵最大值是 ln(N)此时 MI 为 0分布越集中H(P) 越小MI 越接近 1。这个归一化让 MI 不受 N 取值的影响N18 和 N36 计算出的 MI 是可以相互比较的这一点比直接用熵值本身做跨实验比较更稳妥。2.3.1 为什么不能跳过相位箱数量的一致性检查不同包默认的箱数不一致有的用 18 有的用 36。如果你把箱数不同的两组 MI 结果直接做 t 检验会发现低频部分因为箱数少而平均幅度统计量更平稳噪声引起的 MI 偏移也更小导致组间差异被掩盖或放大。我一般在处理多条件对比前会先确认所有数据的 N 完全一致这个检查写在分析的参数配置里不写在主流程里但它是结果可信的前提。3. 用 Python 复现 modulation_index 的最小可运行流程3.1 构造已知耦合的合成信号作为验证基准直接拿真实数据调参很难判断代码算得对不对正确顺序是先构造一个“答案已知”的信号。设计一个 10 秒的合成 LFP采样率 1000 Hztheta 频率 6 Hzgamma 频率 60 Hzgamma 的幅度包络被 theta 相位以 90° 相位偏移调制fs 1000 duration 10 t np.arange(0, duration, 1/fs) f_theta 6.0 f_gamma 60.0 # 低频 theta 成分 theta np.sin(2 * np.pi * f_theta * t) # 高频 gamma 的包络被 theta 调制峰值在 theta 的上升沿 env 1.0 0.8 * np.sin(2 * np.pi * f_theta * t np.pi / 2) # 合成 gamma gamma env * np.sin(2 * np.pi * f_gamma * t) # 混合少量噪声 rng np.random.default_rng(42) lfp theta * 0.3 gamma * 0.7 rng.standard_normal(len(t)) * 0.05这段合成信号里gamma 总量与 theta 分量在同一量级噪声标准差只占总方差的很小比例信噪比足够高MI 应该稳定落在 0.20.4 之间。如果你算出来的 MI 远低于 0.1说明预处理链路某个环节有问题而不是数据没耦合。3.2 核心 MI 计算函数的逐行实现在 synth_data 基础上写一个独立函数计算 Tort MI这一步不依赖任何现成的 PAC 库方便你自己控制每一个变量。def tort_mi(phase, amp_env, n_bins18): # 构造分箱边界 bin_edges np.linspace(-np.pi, np.pi, n_bins 1) bin_idx np.digitize(phase, bin_edges) - 1 bin_idx[bin_idx n_bins] 0 # 把刚好落在 π 上的点归入最后一个箱 # 计算每个相位箱的平均幅度 mean_amp np.zeros(n_bins) for b in range(n_bins): mask bin_idx b if mask.sum() 0: mean_amp[b] amp_env[mask].mean() # 归一化为概率分布 mean_amp mean_amp / mean_amp.sum() # 香农熵 H -(mean_amp * np.log(mean_amp 1e-12)).sum() # 归一化到 [0, 1] mi (np.log(n_bins) - H) / np.log(n_bins) return mi, mean_amp, bin_edgesnp.digitize返回的是“元素落在哪个 bin 区间”的索引左侧闭区间。边界点 π 会被分到最后一个箱之外所以要把索引 n_bins 强制改回 0。mean_amp加1e-12是防止某个箱没有样本时出现 log(0)如果你的相位是均匀分布的基本不会出现空箱但碰到极端滤波或者数据段包含 NaN 时这是必要的保护。3.3 运行验证合成数据的 MI 输出与相-幅图数据phase, amp_env extract_phase_amp( lfp, fsfs, phase_freqf_theta, amp_freqf_gamma, phase_bw2.0, amp_bw20.0 ) mi, mean_amp, bin_edges tort_mi(phase, amp_env) print(fMI {mi:.4f}) # 检查峰值相位箱用于确认耦合相位位置 centers (bin_edges[:-1] bin_edges[1:]) / 2 peak_bin int(np.argmax(mean_amp)) print(f峰值相位箱中心 {centers[peak_bin]:.3f} rad)对应 6 Hz theta 调制 60 Hz gammaMI 输出通常在 0.3 附近峰值相位箱在 π/2 即 1.57 rad 左右。峰值位置由合成信号的调制相位决定而 MI 值只反映“耦合强不强”不反映“耦合发生在哪个相位”。如果你后续要报告“gamma 在 theta 下降沿增强”不能只写 MI必须配合平均幅度分布的峰值位置或方向统计量。3.3.1 常见失败模式MI 接近 0 且峰值相位随机合成信号里噪声权重过高或者滤波器中心频率落在信号频带边缘都会让 MI 退化为接近 0。排查顺序是先画瞬时相位的直方图确认相位均匀覆盖 02π再画幅度包络的频谱确认 60 Hz 附近有能量最后检查滤波后的信号方差是否明显小于原始信号。这个流程比盲目调分箱数量更有效。4. 关键参数怎么设频率对、滤波器阶数与分箱数量4.1 相位频率带宽窄带滤波什么时候失效相位频率的带宽设置直接影响瞬时相位的质量。带宽太窄比如 ±0.5 Hz滤波后的信号能量太小希尔伯特变换提取的相位会被噪声主导带宽太宽比如 ±4 Hz相位信号里混入相邻频率成分相位的时间演化不再干净。经验做法是对 18 Hz 的低频用 ±12 Hz 固定带宽对 830 Hz 的频率用中心频率的 20%30% 作为相对带宽。一个更细致的方式是先用多窗谱估计峰值宽度再设定滤波范围但常规流水线里频率步进 1 Hz、固定带宽的默认做法已经够用。def build_freq_pairs(phase_range(2, 12, 1), amp_range(30, 150, 10), phase_bw2.0, amp_bw20.0): pairs [] for p in range(phase_range[0], phase_range[1] 1, phase_range[2]): for a in range(amp_range[0], amp_range[1] 1, amp_range[2]): if a p * 2: # 幅度频率至少是相位频率的2倍 pairs.append((p, a, phase_bw, amp_bw)) return pairs这个函数生成的是 comodogram 的扫描网格。限制a p * 2是因为相位与幅度太近时滤波器在频域上会重叠算出的 MI 是伪影。幅度频率的步进设为 10 Hz 时gamma 频段的耦合图案比较粗糙但也够看如果要精细定位峰值的中心频率步进要缩到 2 Hz代价是计算时间成倍增加。4.2 滤波器阶次与滤波方向的取舍巴特沃斯滤波器的阶数在高阶时过渡带更陡但通带内的相位响应也更剧烈。filtfilt配合 3 阶是多数电生理分析中默认的配置它对幅度的保真度足够又不会在通带边界产生明显振铃。如果数据长度很短小于 5 个振荡周期高阶滤波器会让边缘数据失真严重此时可以改用 2 阶或者只取滤波后中间 80% 的数据段进行计算。边缘失真在filtfilt下表现为首尾 100200 ms 的剧烈波动我在处理 1 秒 trial 时通常会先拼连续数据再整段滤波然后裁剪掉每段首尾各 0.2 秒比逐 trial 滤波稳定得多。4.3 分箱数量 N 与刺激相位直方图分箱数量的默认值 18 对应每个箱 20°这个分辨率对观察 PAC 的相位调制形状已经足够。把 N 提到 36 可以更精细地定位峰值相位但每个箱内的样本数减半平均幅度的方差变大对含有少量伪迹的数据更不友好。一个折中方式是计算时用 18 箱在确定峰值相位时用第二次插值到 36 箱。后面这个方法对相位精度要求高的场景非常实用。参数默认值取值范围调参方向phase_bw2 Hz13 Hz数据短时取窄信噪比低时取宽amp_bw20 Hz1030 Hz频率越低带宽越窄滤波器阶数325数据短取低阶n_bins181236样本充足时可加大采样率要求≥ 500 Hz越高越好低于 400 Hz 不建议做 100 Hz 的 PAC4.4 数据长度下限与 MI 的偏差问题MI 是有偏估计量数据越短相位箱内样本越少平均幅度的波动越大归一化后的分布越不均匀MI 的期望值被系统性抬高。这是不依赖信号本身是否真的有耦合的非零偏移。为了压制这个偏差需要做到两点一是对每个 trial 单独计算 MI 后取平均不能把所有 trial 拼接成一段长信号再算那样会抹掉 trial 间的相位信息二是做基线校正典型做法是取刺激前的无耦合时段计算 MI把它作为基线电平从刺激时段的 MI 里减掉或除以基线标准差做 z-score。5. 验证与进阶置换检验、comodogram 绘制与条件比较5.1 用置换检验给 MI 一个可靠的 p 值算出一个 MI0.05到底是真耦合还是随机波动最直接的办法是破坏相位-幅度关系构造零分布。把幅度包络在时间轴上按随机偏移循环移位移位量大于低频周期的 3 倍这样破坏耦合但保留幅度本身的统计特性。重复 2001000 次得到随机 MI 分布真实 MI 位于该分布第 95 百分位以上才能算有效。def permutation_test(phase, amp_env, n_perm200, n_bins18, fs1000, shift_range(100, 900)): mi_real, _, _ tort_mi(phase, amp_env, n_bins) shift_samples np.arange(int(shift_range[0] * fs / 1000), int(shift_range[1] * fs / 1000) 1) rng np.random.default_rng(0) mi_null np.zeros(n_perm) n len(amp_env) for i in range(n_perm): shift int(rng.choice(shift_samples)) amp_shifted np.concatenate([amp_env[shift:], amp_env[:shift]]) mi_null[i], _, _ tort_mi(phase, amp_shifted, n_bins) p_value (np.sum(mi_null mi_real) 1) / (n_perm 1) return mi_real, p_value, mi_null移位量的下界必须大于低频周期的数倍否则移位前后的幅度包络仍有部分相位重叠零分布被高估p 值会偏保守。上界不能超过数据长度的 90%避免移位后与原序列重叠太少导致方差异常。p 值计算加 1 是为了把真实 MI 也纳入分布尾部的统计避免出现 p0 的假象。5.2 全频段 comodogram 的绘制从单对到二维图谱单对频率的 MI 只回答一个窄问题全频段扫描才看得出耦合的整体格局。把 4.1 中的build_freq_pairs跑完后把结果整理成二维图像用pcolormesh展示x 轴是相位频率y 轴是幅度频率颜色是 MI 值。正则化时不要用全局最大值做归一化而要用每个幅度频率段内的基线分布做 z-score这样低频耦合强、值大的区域不会把 gamma 段的小值全部压成同一种颜色。5.3 跨条件的 MI 统计比较技巧同一批动物做清醒与麻醉两种状态对比建议先算 z-score 再进入统计模型zMI (MI_real - mean(MI_null)) / std(MI_null)。这一变换消除了基线偏差和方差差异让高耦合状态的 zMI 直接反映超出随机水平的倍数。另一种做法是把 MI 直接作为因变量做混合效应模型但前提是所有条件的 MI 都基于完全相同的相位箱数量与滤波器参数否则组间差异会被参数不一致污染。还有一种小技巧是只对先验频段如 theta 412 Hz × gamma 3080 Hz内的 MI 峰值做条件间比较而不是对全图逐像素做多重比较这样既保留了探索性又控制了假阳性率。本文还有配套的精品资源点击获取 相位-幅度耦合与调制指数的Python计算与统计验证
简介modulation_index-master是一份面向神经科学和脑电信号分析的轻量级MATLAB工具包围绕相位-幅度耦合PAC问题提供计算调制指数MI的函数实现。PAC用于刻画低频神经振荡相位对高频振荡幅度的调控MI是量化这种关联强度的常用统计量尤其适合海马区学习记忆、空间导航等场景的神经信号分析。压缩包共4个文件包含MATLAB函数源码、Markdown说明文档、开源许可证和Git属性配置整体仅3KB结构紧凑计算时通常对低频信号提取相位、对高频信号提取幅度再检验二者的相关性工具已将这一流程封装为可直接调用的函数。目前已有414人浏览学习适合神经科学方向的学生、科研新手以及需要快速上手PAC计算的实验室场景。通过这份资源使用者可免去从零搭建算法框架的时间直接获得封装好的调制指数计算函数和调用说明并理解PAC分析的关键步骤把更多精力放在结果解读与统计检验上。1. modulation_index 计算 PAC一个被低估的神经耦合度量拿到一段局部场电位或脑电数据想回答“高频活动是不是被低频节律牵着走”最常见的技术路径就是计算相位-幅度耦合Phase-Amplitude Coupling, PAC而 PAC 的定量核心就是调制指数Modulation Index, MI。一个容易被忽略的事实MI 的绝对值通常很小0.010.05 的数值在短数据里完全可能是噪声产物直接拿它下结论十有八九要翻车。modulation_index-master这类命名常见于从代码仓库下载的算法包目录它把计算 PAC 的最基础函数封装成了可直接调用的模块。这篇内容会沿着“理论拆解 → 最小实现 → 参数调优 → 统计验证”的顺序把这条链路完整走一遍适合正在处理电生理数据、想搞懂 MI 到底算了什么以及准备把它作为特征做条件比较的工程师和研究者。2. 调制指数计算 PAC 的完整链路从滤波到熵2.1 PAC 的生理信号含义与调制指数的两种计算范式相位-幅度耦合描述的是一个低频振荡比如 theta 412 Hz 或 delta 14 Hz的相位状态对另一个高频振荡gamma 30150 Hz的幅度大小产生调制。在神经科学里海马 theta-gamma 耦合被认为是工作记忆编码的重要机制在皮层信号里慢波相位对 spindle 或 ripples 幅度的调制也常被用来刻画睡眠深度、麻醉深度甚至意识状态。计算 MI 有两大流派。第一派是 Tort 等人 2008 年提出的基于香农熵的方法把相位分箱后对每个相位箱里的平均幅度做归一化再计算该分布相对均匀分布的偏离程度。第二派是 Canolty 等人 2006 年提出的复数向量法把瞬时相位作为角度、瞬时幅度作为模长构建一个复合向量向量的平均模长就是 MI。两者的数学动机不同但信号预处理环节几乎完全一样。Tort 方法对幅度分布更稳健Canolty 方法对短数据更敏感。modulation_index类仓库里最常见的实现是 Tort 变体因为它需要的假设最少且天然落在一个固定区间 [0, 1] 内。2.2 计算 MI 前的信号预处理带通滤波与希尔伯特变换无论选哪种范式第一步都是把原始 LFP 拆成“低频相位信号”和“高频幅度信号”。对于目标相位频率 f_phase需要做一个带通滤波常见带宽是 ±12 Hz对于目标幅度频率 f_amp带宽通常是 ±510 Hz 或相对带宽 20%30%。滤波完成后对低频信号做希尔伯特变换求瞬时相位对高频信号做希尔伯特变换求瞬时幅度包络。import numpy as np from scipy.signal import butter, filtfilt, hilbert def extract_phase_amp(lfp, fs, phase_freq, amp_freq, phase_bw2.0, amp_bw20.0, order3): # 提取低频相位 lo_ph phase_freq - phase_bw / 2 hi_ph phase_freq phase_bw / 2 b_ph, a_ph butter(order, [lo_ph / (fs / 2), hi_ph / (fs / 2)], btypeband) phase_signal filtfilt(b_ph, a_ph, lfp) phase np.angle(hilbert(phase_signal)) # 提取高频幅度包络 lo_amp amp_freq - amp_bw / 2 hi_amp amp_freq amp_bw / 2 b_amp, a_amp butter(order, [lo_amp / (fs / 2), hi_amp / (fs / 2)], btypeband) amp_signal filtfilt(b_amp, a_amp, lfp) amp_env np.abs(hilbert(amp_signal)) return phase, amp_envfiltfilt采用零相位滤波先正向后反向消除了滤波器引入的相位延迟这对相位计算是必须的如果改用lfilter相位会整体偏移MI 值虽然可能不变但后面画相-幅图时峰值位置就是错的。滤波器阶数默认用 3 阶巴特沃斯常规做法是阶数越低越稳高阶级数在锐利谱峰附近容易产生吉布斯振铃。2.3 相位分箱、平均幅度与香农熵的数学表达相位范围 [-π, π] 被均匀分成 N 个箱常见 N18即每个箱 20°。将所有时间点的幅度按瞬时相位落入的箱号分组计算每个箱的平均幅度。设平均幅度序列为 (j)j 1…N归一化为 P(j) (j) / ∑ 。如果相位与幅度完全无关P(j) 应该接近 1/N 的均匀分布如果强耦合P(j) 会呈现一个明显的峰值。MI 的公式是MI (ln(N) ∑ P(j) * ln(P(j))) / ln(N)等价的写法是 MI (ln(N) - H(P)) / ln(N)其中 H(P) 是香农熵最大值是 ln(N)此时 MI 为 0分布越集中H(P) 越小MI 越接近 1。这个归一化让 MI 不受 N 取值的影响N18 和 N36 计算出的 MI 是可以相互比较的这一点比直接用熵值本身做跨实验比较更稳妥。2.3.1 为什么不能跳过相位箱数量的一致性检查不同包默认的箱数不一致有的用 18 有的用 36。如果你把箱数不同的两组 MI 结果直接做 t 检验会发现低频部分因为箱数少而平均幅度统计量更平稳噪声引起的 MI 偏移也更小导致组间差异被掩盖或放大。我一般在处理多条件对比前会先确认所有数据的 N 完全一致这个检查写在分析的参数配置里不写在主流程里但它是结果可信的前提。3. 用 Python 复现 modulation_index 的最小可运行流程3.1 构造已知耦合的合成信号作为验证基准直接拿真实数据调参很难判断代码算得对不对正确顺序是先构造一个“答案已知”的信号。设计一个 10 秒的合成 LFP采样率 1000 Hztheta 频率 6 Hzgamma 频率 60 Hzgamma 的幅度包络被 theta 相位以 90° 相位偏移调制fs 1000 duration 10 t np.arange(0, duration, 1/fs) f_theta 6.0 f_gamma 60.0 # 低频 theta 成分 theta np.sin(2 * np.pi * f_theta * t) # 高频 gamma 的包络被 theta 调制峰值在 theta 的上升沿 env 1.0 0.8 * np.sin(2 * np.pi * f_theta * t np.pi / 2) # 合成 gamma gamma env * np.sin(2 * np.pi * f_gamma * t) # 混合少量噪声 rng np.random.default_rng(42) lfp theta * 0.3 gamma * 0.7 rng.standard_normal(len(t)) * 0.05这段合成信号里gamma 总量与 theta 分量在同一量级噪声标准差只占总方差的很小比例信噪比足够高MI 应该稳定落在 0.20.4 之间。如果你算出来的 MI 远低于 0.1说明预处理链路某个环节有问题而不是数据没耦合。3.2 核心 MI 计算函数的逐行实现在 synth_data 基础上写一个独立函数计算 Tort MI这一步不依赖任何现成的 PAC 库方便你自己控制每一个变量。def tort_mi(phase, amp_env, n_bins18): # 构造分箱边界 bin_edges np.linspace(-np.pi, np.pi, n_bins 1) bin_idx np.digitize(phase, bin_edges) - 1 bin_idx[bin_idx n_bins] 0 # 把刚好落在 π 上的点归入最后一个箱 # 计算每个相位箱的平均幅度 mean_amp np.zeros(n_bins) for b in range(n_bins): mask bin_idx b if mask.sum() 0: mean_amp[b] amp_env[mask].mean() # 归一化为概率分布 mean_amp mean_amp / mean_amp.sum() # 香农熵 H -(mean_amp * np.log(mean_amp 1e-12)).sum() # 归一化到 [0, 1] mi (np.log(n_bins) - H) / np.log(n_bins) return mi, mean_amp, bin_edgesnp.digitize返回的是“元素落在哪个 bin 区间”的索引左侧闭区间。边界点 π 会被分到最后一个箱之外所以要把索引 n_bins 强制改回 0。mean_amp加1e-12是防止某个箱没有样本时出现 log(0)如果你的相位是均匀分布的基本不会出现空箱但碰到极端滤波或者数据段包含 NaN 时这是必要的保护。3.3 运行验证合成数据的 MI 输出与相-幅图数据phase, amp_env extract_phase_amp( lfp, fsfs, phase_freqf_theta, amp_freqf_gamma, phase_bw2.0, amp_bw20.0 ) mi, mean_amp, bin_edges tort_mi(phase, amp_env) print(fMI {mi:.4f}) # 检查峰值相位箱用于确认耦合相位位置 centers (bin_edges[:-1] bin_edges[1:]) / 2 peak_bin int(np.argmax(mean_amp)) print(f峰值相位箱中心 {centers[peak_bin]:.3f} rad)对应 6 Hz theta 调制 60 Hz gammaMI 输出通常在 0.3 附近峰值相位箱在 π/2 即 1.57 rad 左右。峰值位置由合成信号的调制相位决定而 MI 值只反映“耦合强不强”不反映“耦合发生在哪个相位”。如果你后续要报告“gamma 在 theta 下降沿增强”不能只写 MI必须配合平均幅度分布的峰值位置或方向统计量。3.3.1 常见失败模式MI 接近 0 且峰值相位随机合成信号里噪声权重过高或者滤波器中心频率落在信号频带边缘都会让 MI 退化为接近 0。排查顺序是先画瞬时相位的直方图确认相位均匀覆盖 02π再画幅度包络的频谱确认 60 Hz 附近有能量最后检查滤波后的信号方差是否明显小于原始信号。这个流程比盲目调分箱数量更有效。4. 关键参数怎么设频率对、滤波器阶数与分箱数量4.1 相位频率带宽窄带滤波什么时候失效相位频率的带宽设置直接影响瞬时相位的质量。带宽太窄比如 ±0.5 Hz滤波后的信号能量太小希尔伯特变换提取的相位会被噪声主导带宽太宽比如 ±4 Hz相位信号里混入相邻频率成分相位的时间演化不再干净。经验做法是对 18 Hz 的低频用 ±12 Hz 固定带宽对 830 Hz 的频率用中心频率的 20%30% 作为相对带宽。一个更细致的方式是先用多窗谱估计峰值宽度再设定滤波范围但常规流水线里频率步进 1 Hz、固定带宽的默认做法已经够用。def build_freq_pairs(phase_range(2, 12, 1), amp_range(30, 150, 10), phase_bw2.0, amp_bw20.0): pairs [] for p in range(phase_range[0], phase_range[1] 1, phase_range[2]): for a in range(amp_range[0], amp_range[1] 1, amp_range[2]): if a p * 2: # 幅度频率至少是相位频率的2倍 pairs.append((p, a, phase_bw, amp_bw)) return pairs这个函数生成的是 comodogram 的扫描网格。限制a p * 2是因为相位与幅度太近时滤波器在频域上会重叠算出的 MI 是伪影。幅度频率的步进设为 10 Hz 时gamma 频段的耦合图案比较粗糙但也够看如果要精细定位峰值的中心频率步进要缩到 2 Hz代价是计算时间成倍增加。4.2 滤波器阶次与滤波方向的取舍巴特沃斯滤波器的阶数在高阶时过渡带更陡但通带内的相位响应也更剧烈。filtfilt配合 3 阶是多数电生理分析中默认的配置它对幅度的保真度足够又不会在通带边界产生明显振铃。如果数据长度很短小于 5 个振荡周期高阶滤波器会让边缘数据失真严重此时可以改用 2 阶或者只取滤波后中间 80% 的数据段进行计算。边缘失真在filtfilt下表现为首尾 100200 ms 的剧烈波动我在处理 1 秒 trial 时通常会先拼连续数据再整段滤波然后裁剪掉每段首尾各 0.2 秒比逐 trial 滤波稳定得多。4.3 分箱数量 N 与刺激相位直方图分箱数量的默认值 18 对应每个箱 20°这个分辨率对观察 PAC 的相位调制形状已经足够。把 N 提到 36 可以更精细地定位峰值相位但每个箱内的样本数减半平均幅度的方差变大对含有少量伪迹的数据更不友好。一个折中方式是计算时用 18 箱在确定峰值相位时用第二次插值到 36 箱。后面这个方法对相位精度要求高的场景非常实用。参数默认值取值范围调参方向phase_bw2 Hz13 Hz数据短时取窄信噪比低时取宽amp_bw20 Hz1030 Hz频率越低带宽越窄滤波器阶数325数据短取低阶n_bins181236样本充足时可加大采样率要求≥ 500 Hz越高越好低于 400 Hz 不建议做 100 Hz 的 PAC4.4 数据长度下限与 MI 的偏差问题MI 是有偏估计量数据越短相位箱内样本越少平均幅度的波动越大归一化后的分布越不均匀MI 的期望值被系统性抬高。这是不依赖信号本身是否真的有耦合的非零偏移。为了压制这个偏差需要做到两点一是对每个 trial 单独计算 MI 后取平均不能把所有 trial 拼接成一段长信号再算那样会抹掉 trial 间的相位信息二是做基线校正典型做法是取刺激前的无耦合时段计算 MI把它作为基线电平从刺激时段的 MI 里减掉或除以基线标准差做 z-score。5. 验证与进阶置换检验、comodogram 绘制与条件比较5.1 用置换检验给 MI 一个可靠的 p 值算出一个 MI0.05到底是真耦合还是随机波动最直接的办法是破坏相位-幅度关系构造零分布。把幅度包络在时间轴上按随机偏移循环移位移位量大于低频周期的 3 倍这样破坏耦合但保留幅度本身的统计特性。重复 2001000 次得到随机 MI 分布真实 MI 位于该分布第 95 百分位以上才能算有效。def permutation_test(phase, amp_env, n_perm200, n_bins18, fs1000, shift_range(100, 900)): mi_real, _, _ tort_mi(phase, amp_env, n_bins) shift_samples np.arange(int(shift_range[0] * fs / 1000), int(shift_range[1] * fs / 1000) 1) rng np.random.default_rng(0) mi_null np.zeros(n_perm) n len(amp_env) for i in range(n_perm): shift int(rng.choice(shift_samples)) amp_shifted np.concatenate([amp_env[shift:], amp_env[:shift]]) mi_null[i], _, _ tort_mi(phase, amp_shifted, n_bins) p_value (np.sum(mi_null mi_real) 1) / (n_perm 1) return mi_real, p_value, mi_null移位量的下界必须大于低频周期的数倍否则移位前后的幅度包络仍有部分相位重叠零分布被高估p 值会偏保守。上界不能超过数据长度的 90%避免移位后与原序列重叠太少导致方差异常。p 值计算加 1 是为了把真实 MI 也纳入分布尾部的统计避免出现 p0 的假象。5.2 全频段 comodogram 的绘制从单对到二维图谱单对频率的 MI 只回答一个窄问题全频段扫描才看得出耦合的整体格局。把 4.1 中的build_freq_pairs跑完后把结果整理成二维图像用pcolormesh展示x 轴是相位频率y 轴是幅度频率颜色是 MI 值。正则化时不要用全局最大值做归一化而要用每个幅度频率段内的基线分布做 z-score这样低频耦合强、值大的区域不会把 gamma 段的小值全部压成同一种颜色。5.3 跨条件的 MI 统计比较技巧同一批动物做清醒与麻醉两种状态对比建议先算 z-score 再进入统计模型zMI (MI_real - mean(MI_null)) / std(MI_null)。这一变换消除了基线偏差和方差差异让高耦合状态的 zMI 直接反映超出随机水平的倍数。另一种做法是把 MI 直接作为因变量做混合效应模型但前提是所有条件的 MI 都基于完全相同的相位箱数量与滤波器参数否则组间差异会被参数不一致污染。还有一种小技巧是只对先验频段如 theta 412 Hz × gamma 3080 Hz内的 MI 峰值做条件间比较而不是对全图逐像素做多重比较这样既保留了探索性又控制了假阳性率。本文还有配套的精品资源点击获取