简介这是一套面向静息态功能磁共振研究者的MATLAB工具箱聚焦HRF血流动力学响应函数估计、反卷积与连通性分析支持独立运行及作为SPM插件使用适合具备一定MATLAB与神经影像基础、希望深入探究静息态HRF动态特性的科研人员。资源包共148个文件以130个.m脚本为核心辅以.mat数据、.txt与.md说明文档、.zip依赖包、.pdf教程及.nii、.gii影像示例等压缩包约5.82MB目录结构清晰便于按模块查阅与二次开发。工具箱内置演示数据集与demo脚本覆盖从数据准备到反卷积、连通性分析的完整流程并附有分步说明与Python独立版本、BIDS应用的相关指引可帮助读者快速上手静息态HRF分析、理解算法实现细节并迁移到自身研究场景。目前已有455人学习下载适合需要系统掌握HRF反卷积方法的进阶用户参考。1. 静态 HRF 反卷积到底在解决什么问题从 BOLD 信号里把神经活动“抠”出来如果你做过静息态 fMRI 分析大概率遇到过这个场景数据跑完预处理准备做功能连接结果发现不同被试、不同脑区的 BOLD 信号延迟差异很大有的区域响应快有的慢半拍直接算相关就把这些差异全混进去了。更麻烦的是静息态没有刺激时序你根本不知道神经事件什么时候发生传统 GLM 那套“已知刺激函数卷积 HRF”的思路直接失效。二阶梯度 matlab 代码 rsHRF 这个工具箱干的就是这件事在静息态下用盲反卷积把每个体素的 BOLD 信号拆成“神经活动事件”和“血流动力学响应函数”两部分同时估计出 HRF 的形状参数峰值时间、响应幅度、延迟等。它把 HRF 当成一个带二阶梯度约束的模型来拟合避免解出物理上不合理的振荡解。适合谁用做静息态功能连接、有效连接、或者想校正 HRF 变异性对连接结果影响的人。如果你只是做任务态 GLM这个工具箱不是给你用的。2. rsHRF 的模型底子二阶梯度约束和 HRF 参数化怎么选2.1 为什么是二阶梯度而不是直接解卷积盲反卷积本身是个欠定问题给定一个 BOLD 时间序列理论上可以有无穷多组“神经事件 HRF”的组合。常见做法是加正则项比如 L1 稀疏约束神经事件或者假设 HRF 属于某个参数族。rsHRF 走的是后一条路但它对 HRF 的拟合加了二阶梯度惩罚。直白说二阶梯度约束的是 HRF 曲线的“平滑程度”。一阶梯度管斜率二阶梯度管曲率。如果只约束一阶HRF 可能出现尖刺约束到二阶相当于要求 HRF 的曲率变化不能太剧烈解出来的曲线更接近真实血流动力学响应的钟形。这个约束在频域上等价于一个低通滤波把高频的数值噪声压掉。我一般会先跑默认参数如果发现某些体素的 HRF 峰值时间估计得离谱比如小于 2 秒或大于 10 秒再回头调二阶梯度的权重。2.2 HRF 基函数的选择双伽马 vs 逆对数正态rsHRF 支持几种 HRF 参数化方式最常用的是双伽马函数和逆对数正态。双伽马就是 SPM 默认那个两个伽马函数相减参数是峰值延迟和 undershoot 延迟。逆对数正态在拟合早期响应时更灵活但参数解释起来麻烦一点。选哪个如果你后续要做组间比较建议用双伽马因为参数物理意义明确峰值时间直接对应 BOLD 响应的达峰时间。如果只是做单被试的 HRF 估计逆对数正态拟合残差通常更小。我自己的习惯是先跑双伽马看拟合优度R²低于 0.6 的体素比例如果超过 20%换逆对数正态再跑一遍对比。2.3 工具箱的目录结构和依赖解压后你会看到几个核心文件夹rsHRF主函数目录、HRF_est估计模块、connectivity连通性分析模块、utils工具函数。依赖方面它需要 MATLAB 的统计工具箱和信号处理工具箱。如果你用的是 Octave部分函数会报错因为用到了fit和nlinfit这类 MATLAB 特有的拟合函数。提示跑之前先执行addpath(genpath(rsHRF_root))把整个目录树加进路径否则hrf_estimation会找不到子函数。3. 从 NIfTI 到 HRF 参数图完整跑一遍的代码和参数说明3.1 数据准备和预处理边界rsHRF 的输入是预处理后的 NIfTI 文件。预处理至少要做完 slice timing 和 motion correction是否做空间平滑取决于你的目的如果要做体素级 HRF 估计平滑核不要超过 6mm否则 HRF 参数会被空间模糊掉如果只做 ROI 级分析8mm 平滑问题不大。去线性漂移和带通滤波要注意rsHRF 内部有自己的去趋势步骤如果你在预处理里已经做了 0.01-0.1Hz 带通建议把带通改成 0.008-0.15Hz给 HRF 估计留一点低频信息。我见过有人用 0.01-0.08Hz 的窄带数据跑 rsHRF结果 HRF 峰值时间全被压到 4 秒附近这就是带通把 HRF 的低频成分滤掉了。3.2 单被试 HRF 估计的完整调用下面这段代码是我常用的单被试估计流程输入是预处理后的 4D NIfTI输出是 HRF 参数图和神经活动时间序列。% 添加工具箱路径 addpath(genpath(/path/to/rsHRF)); % 输入输出设置 input_nii sub-01_task-rest_bold.nii; output_dir sub-01_rsHRF_out; mask_nii sub-01_mask.nii; % 核心参数 cfg struct(); cfg.TR 2; % 重复时间单位秒必须和数据一致 cfg.T 200; % 时间点数留空则自动从数据读取 cfg.hrf canonical; % HRF 基函数canonical 双伽马IL 逆对数正态 cfg.estimation wiener; % 估计方法wiener 或 irls cfg.lag 0; % 允许的 HRF 延迟搜索范围0 表示不搜索 cfg.threshold 0.1; % 神经事件幅度阈值低于此值视为无事件 cfg.penalty 0.1; % 二阶梯度惩罚权重越大曲线越平滑 % 运行估计 hrf_estimation(input_nii, output_dir, mask_nii, cfg);逻辑说明hrf_estimation是主入口函数它先读取 4D 数据对每个体素的时间序列做盲反卷积。cfg.TR和cfg.T决定时间轴如果 TR 设错HRF 峰值时间会整体缩放。cfg.penalty是二阶梯度约束的强度默认 0.1调到 0.5 以上 HRF 会变得非常平滑但可能欠拟合调到 0.01 以下可能出现振荡解。参数说明cfg.estimation选wiener用的是维纳滤波反卷积速度快但假设噪声平稳选irls是迭代重加权最小二乘慢但对非平稳噪声更稳。我一般先用wiener跑一遍看整体质量如果个别体素残差大再单独用irls重跑。3.3 输出文件解读和参数图提取跑完之后output_dir里会生成几个文件HRF_parameters.nii包含峰值时间、响应幅度、延迟等参数neural_events.nii是估计出的神经活动时间序列fit_quality.nii是拟合优度 R²。提取峰值时间图的代码% 读取 HRF 参数图 hrf_param_nii fullfile(output_dir, HRF_parameters.nii); V spm_vol(hrf_param_nii); Y spm_read_vols(V); % 假设第一个 volume 是峰值时间 peak_time Y(:,:,:,1); % 保存为新的 NIfTI V_out V(1); V_out.fname fullfile(output_dir, peak_time.nii); spm_write_vol(V_out, peak_time);逻辑说明spm_vol和spm_read_vols是 SPM 的函数如果你没装 SPM可以用niftiread替代。HRF_parameters.nii的第四个维度顺序是固定的峰值时间、响应幅度、延迟、宽度。别搞错顺序否则后续统计全错。3.4 组水平统计的输入准备单被试跑完后组水平分析需要把每个被试的峰值时间图配准到标准空间然后做体素级单样本 t 检验。常见做法是用 DARTEL 或 FNIRT 做非线性配准配准时用峰值时间图作为输入不要用原始 BOLD 图因为原始 BOLD 的对比度在配准时不如参数图稳定。配准后用 SPM 的spm_glm或fitglm做组分析。我一般会额外把fit_quality.nii作为协变量放进模型把拟合差的体素权重降低。这一步不是必须的但能减少噪声体素对组结果的污染。4. 连通性分析用 HRF 参数做加权还是做校正4.1 两种连通性分析思路的区别rsHRF 的连通性模块提供两种模式一种是把 HRF 参数比如峰值时间作为加权因子对功能连接矩阵做逐边加权另一种是用估计出的神经活动时间序列替代原始 BOLD 做连接相当于把 HRF 变异性校正掉。第一种思路适合研究“HRF 变异性本身是否与连接强度相关”的问题。比如你假设某脑区的 HRF 延迟越长它和另一个脑区的连接越弱那就用加权模式。第二种思路适合研究“校正 HRF 后连接结果是否更稳定”的问题。我自己的经验是如果被试间 HRF 峰值时间标准差超过 1 秒校正后的连接矩阵和原始矩阵差异会很明显这时候用校正模式更合理。4.2 连通性分析的代码实现% 读取神经活动时间序列 neural_nii fullfile(output_dir, neural_events.nii); V spm_vol(neural_nii); Y spm_read_vols(V); % 提取 ROI 时间序列假设 ROI 掩膜已准备好 roi_mask spm_read_vols(spm_vol(roi_mask.nii)); n_roi max(roi_mask(:)); ts zeros(size(Y,4), n_roi); for r 1:n_roi idx find(roi_mask r); ts(:,r) squeeze(mean(Y(repmat(idx,1,size(Y,4)) ... (0:size(Y,4)-1)*numel(roi_mask)), 1)); end % 计算 Pearson 相关矩阵 conn_mat corr(ts); % 如果要做 HRF 加权读取峰值时间图 peak_time spm_read_vols(spm_vol(fullfile(output_dir, peak_time.nii))); roi_peak zeros(n_roi,1); for r 1:n_roi idx find(roi_mask r); roi_peak(r) mean(peak_time(idx)); end % 加权峰值时间差异越大权重越低 weight_mat 1 ./ (1 abs(roi_peak - roi_peak)); conn_weighted conn_mat .* weight_mat;逻辑说明neural_events.nii是 4D 文件第四个维度是时间。提取 ROI 时间序列时用均值也可以用第一主成分但均值更稳。加权矩阵weight_mat的设计不唯一这里用的是反比函数你也可以用高斯核。关键是权重要有下界否则峰值时间差异大的边权重会趋近于零连接矩阵变得太稀疏。参数说明roi_mask.nii的标签要从 1 开始连续编号如果有跳号max(roi_mask(:))会多出空 ROI。跑之前用unique(roi_mask(:))检查一下。4.3 结果可视化和统计检验连接矩阵算完后组水平比较用 Fisher z 变换把相关系数转成近似正态再做双样本 t 检验。如果用了加权注意加权后的连接值不再是相关系数不能直接做 Fisher z要先做标准化。可视化我一般用imagesc画矩阵配合colorbar和caxis固定色标范围方便跨被试对比。别用默认的自动色标否则被试间的视觉差异会被色标缩放掩盖。5. 避坑与排查rsHRF 跑崩的五个常见原因5.1 报错 “Undefined function nlinfit”现象运行hrf_estimation时 MATLAB 报错找不到nlinfit。原因没装统计工具箱或者工具箱路径没加全。解决ver命令查看已安装工具箱确认 Statistics and Machine Learning Toolbox 在列表里。如果装了但还报错执行restoredefaultpath后重新addpath。5.2 HRF 峰值时间全是 0 或 NaN现象输出的peak_time.nii里大部分体素是 0 或 NaN。原因掩膜文件和数据没对齐或者cfg.threshold设得太高所有神经事件都被阈值滤掉了。解决先用spm_check_registration看掩膜和数据是否重叠然后把cfg.threshold降到 0.01 再跑一遍。如果还是全 0检查cfg.TR是否设成了 0。5.3 拟合优度 R² 普遍低于 0.3现象fit_quality.nii均值只有 0.2 左右。原因数据预处理带通太窄把 HRF 的低频成分滤掉了或者cfg.penalty太大HRF 被过度平滑。解决把带通改成 0.008-0.15Hzcfg.penalty降到 0.05。如果还不行换cfg.hrf IL试试。5.4 跑得特别慢一个被试要半小时以上现象hrf_estimation卡在某个体素不动。原因cfg.estimation irls在噪声大的体素上迭代次数会爆炸或者数据没做掩膜全脑几万个体素都在跑。解决先用wiener快速跑一遍然后用fit_quality做掩膜只对 R² 0.5 的体素用irls重跑。另外确认掩膜文件是二值化的别把概率掩膜直接丢进去。5.5 组分析时配准效果差现象峰值时间图配准到标准空间后脑区边界模糊。原因用原始 BOLD 图做的配准对比度不适合参数图。解决用峰值时间图本身做配准输入或者在配准前对峰值时间图做一次轻度平滑4mm提高配准的鲁棒性。配准后记得检查一下几个关键脑区的峰值时间是否合理比如初级视觉皮层应该在 5-6 秒附近。6. 进阶技巧用 HRF 参数做被试筛选和模型诊断跑完一批被试后别急着做组统计。我习惯先做一轮被试级诊断把 HRF 估计质量差的被试标出来。具体做法是计算每个被试全脑fit_quality的中位数如果某个被试的中位数低于所有被试中位数的 2 个标准差就标记为异常。这一步能筛掉那些头动大或者预处理有问题的被试比单纯看头动参数更直接。另一个技巧是用 HRF 峰值时间的空间分布做模型诊断。如果某个被试的峰值时间图在灰质和白质之间没有对比度说明反卷积没收敛。正常情况灰质的峰值时间应该在 4-7 秒白质因为血流量低峰值时间估计会偏大或者拟合优度很低。如果灰质白质一个样回去检查掩膜和预处理。还有一个我踩过的坑用 rsHRF 估计出的神经活动时间序列做功能连接时别忘了这个序列是反卷积出来的它的自相关结构和原始 BOLD 不一样。做统计检验时自由度不能直接用时间点数减 2要用有效自由度。常见做法是用 AR(1) 模型估计有效自由度或者用置换检验。我一般用置换检验跑 5000 次虽然慢但省心。从那以后我每次跑 rsHRF 都强制走一遍诊断流程先看fit_quality的空间分布再看峰值时间的灰质白质对比最后检查神经活动序列的自相关。这三步走完后面的组统计才敢往下做。希望帮到你。本文还有配套的精品资源点击获取