简介：面向静息态功能磁共振研究者的 rsHRF 工具箱，聚焦血氧动力学响应函数的估计与反卷积，并提供连通性分析支持。代码以 MATLAB 实现，既可作为独立工具运行，也能以 SPM 插件形式接入既有预处理流程，适合具备一定脑影像数据处理基础、希望从峰值统计转向逐体素 HRF 参数建模的研究生与科研人员，用于静息态 HRF 形态差异、延迟估计及后续连通性分析等场景。压缩包共 148 个文件、约 5.82MB，主体为 130 个 .m 函数与配置脚本，另含 .mat 示例数据、.txt 与 .md 说明文档、.pdf 手册、.gii 与 .nii 影像样例、.yml 环境配置及 LICENSE 许可文件，结构完整、便于二次开发。随包附带的演示数据集与 demo 脚本可直接跑通估计与反卷积流程，SPM 配置入口与数据组织示例也让插件式调用更易上手；对希望复现静息态 HRF 分析、排查参数设置问题的读者具备直接参考价值。目前已有 455 人学习。
1. 从一张“糊掉的 HRF”说起：rsHRF 到底解决什么问题
做完静息态预处理，跑 DCM 或者 Granger，结果忽好忽坏，回头一查，问题常常出在那个谁都没动过的默认 HRF 上。SPM 给的是双伽马标准曲线，全脑一套参数、所有人一套参数，可真实血流动力学响应在个体之间、脑区之间差得很远——峰值时间从 4 秒到 8 秒我都见过。rsHRF 就是冲着这件事来的：它用反卷积从静息态 BOLD 信号本身把 HRF 估出来，再把 BOLD 还原成神经事件序列，接着做有效连接分析。整套东西是系统开源的 MATLAB 实现，依赖 SPM12，适合已经跑通预处理、想往下做个体化 HRF 和连接分析的开发者。
2. 目录结构与运行环境：把 rsHRF 和 SPM12 稳妥地装到一起
2.1 源码包里都有什么，先按功能认，别按文件名背
拿到这个包，第一件事不是急着跑，而是把目录过一遍。rsHRF 的组织方式很典型的 MATLAB 工具箱风格：根目录扔入口函数和核心算法，辅助函数单独放，示例数据单独放。不同版本的文件名会差一两个字，但你按功能去认，一眼就能对上号。
| 位置 | 大致内容 | 你改动它的频率 |
|---|---|---|
根目录rsHRF_*.m | 批处理入口、HRF 估计主函数、参数初始化、连接分析 | 高 |
| 子函数目录 | AR 预白化、基函数生成、数据裁切、迭代求解 | 中 |
| 示例数据目录 | 小体量的 demo 数据，用来验证环境通不通 | 低 |
| 文档与注释 | 函数头注释、参数说明、方法出处 | 中 |
几个你一定会打交道的函数，职责大致是这样：rsHRF_main是批处理入口，顶部有一段配置区，改数据目录、TR、掩模路径就能整批跑；rsHRF_est是真正干活的那个，输入时间序列和参数，输出 HRF 参数、去卷积后的神经事件、HRF 曲线；rsHRF_paron负责生成默认参数结构体，所有默认值都在这里，改参数前先来这里看一眼比翻文档快；rsHRF_est_GS做的是网格搜索版本，参数化方法下用它找更优的 HRF 形状参数；rsHRF_conn系列则是在去卷积结果上做连接分析。
提示：这套代码是系统开源的纯
.m实现，没有 p-code 加密，函数头注释写得比多数工具箱详细。遇到不确定的字段，与其猜，不如edit rsHRF_paron直接看默认值，这条路比查任何二手资料都靠谱。
2.2 SPM12 依赖与 addpath 顺序，这里翻过车
rsHRF 不自己实现 NIfTI 读写、GLM 设计矩阵、AR 模型这些底层活，它直接调 SPM12 的函数。所以 SPM12 必须先装好、先上路径。MATLAB 的搜索路径是有顺序的，谁在前面谁优先，这一点在混装多个工具箱的机器上尤其致命。
% startup_rsHRF.m —— 我习惯每次开 MATLAB 先跑这一段 addpath('/home/user/work/matlab/spm12'); % SPM12 必须排在前面 spm('defaults','FMRI'); % 初始化 SPM 全局变量，不写这句后面全是空值 addpath(genpath('/home/user/work/matlab/rsHRF')); % 递归加，子函数目录一起进来 savepath; % 确认没问题再固化，出问题好回退第一行的顺序不能反。如果先把 rsHRF 递归加进去，而包里恰好带了某个旧版spm_*.m同名函数，它会遮蔽 SPM12 的正版实现，报出来的错往往指向一个你根本没调用过的函数，排查起来极其难受。第二行的spm('defaults','FMRI')也不能省，它设置的是 SPM 的全局结构体，省了之后读写头信息、取体素尺寸都可能拿到空值或者默认值。
验证环境通不通，不用上真实数据，直接跑包里的 demo。跑通了说明路径、SPM 依赖、数据读写这条链没问题；跑不通就老老实实从报错的第一行往上找，十有八九是路径顺序或者 SPM 没初始化。
2.3 数据要求：预处理到哪一步才能喂进去
这一点必须先说清楚，否则后面全是白忙。rsHRF 吃的是已经做完常规预处理的 4D NIfTI：时间层校正、头动校正、配准到标准空间、空间平滑，这些都得先做完。它内部会自己做高通滤波（由P.cutoff控制），所以你不需要提前滤一遍，重复滤波反而会把低频信息切掉两次。
TR 必须准确。这个参数不是随便填的，它直接决定了去卷积的时间栅格和 HRF 的采样密度。填错 TR，估计出来的峰值时间会整体偏移，而且偏得很隐蔽——曲线形状看着还挺像那么回事，只是峰值落在了 3 秒而不是 5 秒，你不细看根本发现不了。
掩模方面，体素级全脑跑是可行的，但代价很大：每个体素都要解一次反卷积，几万体素乘上时间点，单被试跑几个小时很正常。我一般先用 ROI 或者脑网络模板把数据降到几十上百个节点，节点内取平均时间序列再送进去。降维会损失空间特异性，但对于连接分析这个下游任务来说，节点级别的 HRF 估计已经够用了。
3. rsHRF_est 核心拆解：一次反卷积到底走了哪几步
3.1 参数结构 P：每个字段都在影响什么
所有行为都由一个结构体P控制，用rsHRF_paron生成默认值再覆盖字段，是改动最小、最不容易漏参的写法。
P = rsHRF_paron(struct( ... 'TR', 2, ... % 重复时间，秒，必须和采集一致 'T', (0:2:(T-1)*2)', ...% 时间点向量，长度要等于数据行数 'T0', 1, ... % 起始时间偏移，通常留 1 个 TR 避开初始瞬态 'dt', 0.2, ... % 去卷积的微时间分辨率，一般是 TR/10 'AR', 1, ... % 是否做 AR 预白化，1 开启 'estimation', 'canonical', ... % 估计方法：canonical / fir / sFIR 'TD_DD', 1, ... % 是否带时间导数与弥散导数 'lag', 4, ... % 基函数滞后阶数 'cutoff', 128, ... % 高通截止周期，秒 'local', 1, ... % 局部迭代优化开关 'mid', 1)); % 中间层优化开关几个字段值得单独说。dt是整个去卷积的时间分辨率，设得比 TR 细是因为 BOLD 采样太粗，直接在 TR 栅格上解卷积会把神经事件的时序钉死在 2 秒一格上，峰值位置误差可能到半个 TR。常见的做法是取TR/10，再细收益递减，计算量却线性上涨。
cutoff的单位是秒不是赫兹，128 秒对应大约 0.0078 Hz，这是静息态功能连接的常用下限。填成 0.01 这种赫兹值，滤波会把大部分信号切掉，出来的 HRF 平得像一条直线——这个坑我见过不止一次。
AR开启后会对残差做自回归预白化。fMRI 时间序列的自相关很强，不白化的话最小二乘的标准误会被低估，HRF 幅度估计偏大。代价是多一层计算，单被试耗时大概增加两三成，我认为这个代价值得付。
lag控制基函数向后延伸多少个时间单位。设太小，HRF 的负向回落（undershoot）截不完整，估计出来的曲线尾巴是断的；设太大，设计矩阵列数膨胀，容易过拟合噪声。默认值在多数 2 秒 TR 的数据上是合适的，TR 长于 3 秒时可以适当加。
3.2 三种估计方法：canonical、fir、sFIR 怎么选
这是整个工具箱最核心的分叉点，选错了后面全歪。
| 方法 | 基函数 | 优点 | 代价 | 适用场景 |
|---|---|---|---|---|
canonical | 标准双伽马 + 时间导数 + 弥散导数 | 参数少、极稳、可跨被试比较 | 形状被约束在标准曲线附近 | 组水平分析、样本量大 |
fir | 一组有限脉冲响应基 | 形状完全自由，能抓双峰、初始下凹 | 参数多，噪声大时容易抖 | 单被试深挖、事件相关设计 |
sFIR | FIR 加平滑约束 | 兼顾自由度与稳定性 | 平滑强度需要调 | 静息态个体化 HRF 的折中方案 |
canonical是我在组水平分析里的默认选择。它假设 HRF 形状接近标准双伽马，只允许在幅度、峰值时间偏移、弥散程度三个维度上变化，参数少所以抗噪，被试之间也能直接比这三个参数。缺点是如果真实 HRF 有明显双峰或者初始下凹，它拟合不出来，残差里会留下系统性结构。
fir走另一个极端，把 HRF 当成一条自由曲线来估。自由度上去了，什么形状都能拟合，但静息态数据本来信噪比就低，自由度高意味着噪声也被拟进去了，估计出来的 HRF 常常是锯齿状的，看着就不像生理信号。要用的话，得配合后续平滑或者足够长的数据。
sFIR是这两者之间的务实选择，在 FIR 基础上加平滑惩罚，允许形状偏离标准曲线但不允许剧烈抖动。我在单被试分析里更倾向用它，原因是静息态没有刺激时刻作为锚点，去卷积本身就是一个病态问题，多一点约束不是坏事。
3.3 调用与输出：三个返回值分别是什么
数据读进来、转成「时间点 × 体素」的矩阵之后，一行调用就能出结果。
% 读 4D 数据并整理成 T×N V = spm_vol('sub-01_space-MNI_desc-preproc_bold.nii'); Y = spm_read_vols(V); % [x, y, z, T] Y = reshape(Y, [], size(Y,4))'; % [T, Nvox] % ROI 内取平均，把 Nvox 降到节点数 mask = spm_read_vols(spm_vol('roi_node.nii')) > 0; y = mean(Y(:, mask), 2); % [T, 1] % 核心调用 [beta_hrf, event_bold, hrf] = rsHRF_est(y, P.T, P.T0, P, P.estimation);beta_hrf是拟合出来的 HRF 参数，行对应体素（或节点），列的含义随方法变化：canonical 下第一列是幅度，第二列是峰值时间偏移，第三列是弥散宽度，TD_DD关掉时列数会减少。event_bold是去卷积得到的神经事件序列，尺寸和输入y一样，这是后面做连接分析的原料。hrf是重建出来的 HRF 曲线本身，长度由P.T和dt决定，画出来就是那条你熟悉的波形。
这三个返回值里，hrf是最该先看的。不用做任何统计，直接plot(P.T(1:size(hrf,1)), hrf(:,1))画一条出来，峰值位置对不对、有没有负向回落、曲线是不是抖成锯齿，肉眼就能判断。我一般会先抽三五个节点画一遍，确认形状合理再往下跑批处理，这比跑完一百个被试再回来查问题省太多时间。
3.4 迭代与 AR 白化：为什么不是一次最小二乘就完事
初次解出来的 HRF 和神经事件是互相依赖的：要估 HRF 得知道事件，要估事件得知道 HRF。工具箱的处理方式是迭代——先用一个初始 HRF 反卷积得到事件估计，再用这些事件重新拟合 HRF，来回几轮直到收敛。local和mid这两个开关控制的就是迭代策略的不同层次。
这个过程收敛与否，直接影响结果的可信度。我习惯在批处理脚本里把迭代次数和最终残差都记下来，跑完扫一眼有没有异常大的残差。残差突然飙高的被试，往往是头动大、配准差或者 TR 填错，数据本身有问题，不是算法的问题。
注意：迭代收敛不代表结果正确。如果输入数据的高通滤波截止频率设得过于激进，或者 TR 填错，迭代照样会收敛到一个看起来很像 HRF、实际上峰值整体偏移的曲线上。参数检查必须在跑之前做，不能指望算法自己兜底。
4. 批处理与结果质检：从单被试脚本到一百个被试
4.1 批处理循环：把 try-catch 写进去
单被试跑通之后，剩下的就是把这段逻辑包成循环。这一步没什么玄学，但有两个细节决定了你第二天早上看到的是完整结果还是一堆中断。
subs = {'sub-01','sub-02','sub-03'}; % 被试列表 outdir = fullfile(pwd, 'hrf_out'); if ~exist(outdir, 'dir'), mkdir(outdir); end for i = 1:numel(subs) try f = fullfile(pwd, 'data', [subs{i} '_bold.nii']); V = spm_vol(f); Y = spm_read_vols(V); Y = reshape(Y, [], size(Y,4))'; y = mean(Y(:, mask), 2); [beta_hrf, event_bold, hrf] = rsHRF_est(y, P.T, P.T0, P, P.estimation); % 把 P 一起存下来，否则三个月后你不知道这批结果用的什么参数 save(fullfile(outdir, [subs{i} '_hrf.mat']), ... 'beta_hrf','event_bold','hrf','P','-v7.3'); fprintf('[ OK ] %s\n', subs{i}); catch ME fprintf('[FAIL] %s : %s\n', subs{i}, ME.message); end endtry-catch不是可有可无的装饰。一百个被试跑四五个小时，中间只要有一个数据文件损坏或者掩模维度对不上，没有容错的话整个循环就停在那里，前面跑完的也得重来。把错误信息打出来、继续下一个，早上起来看日志定位失败被试，补跑一遍就行。
save那句里把P一起存进去，是我踩过坑之后养成的习惯。半年前跑的一批结果，回头想复核，参数结构体没存，dt用的是 0.2 还是 0.25 完全想不起来，只能整批重跑。多存一个结构体，省的是未来的自己。
数据量大、机器内存够的话，把for换成parfor能明显缩短时间，但要注意两点：一是每个 worker 都要能访问到 SPM12 的路径，二是不要在循环里做savepath这类全局操作。SPM 的部分函数不是完全线程安全的，parfor跑出随机崩溃的话，退回for观察，别硬扛。
4.2 输出组织：文件名自己定，别依赖默认
工具箱自带的批处理入口会按它自己的规则命名输出文件，但在自己写的流水线里，文件名最好由你控制。我习惯把方法名和关键参数编进文件名，sub-01_canonical_dt02.mat这种，事后一眼能看出来这批结果是什么配置下跑的。
输出的三样东西，用途各不相同：beta_hrf进组水平统计，做被试间比较；event_bold进连接分析，是下一章的主角；hrf进质控，画图看形状。三者分开存、分开处理，不要混在一个大矩阵里。
4.3 质检：把峰值时间和半高宽画出来
批处理跑完，最忌讳直接拿去做统计。先做一轮形态质检，成本很低但能捞出一批脏数据。
peak = nan(numel(subs),1); t2p = nan(numel(subs),1); for i = 1:numel(subs) S = load(fullfile(outdir, [subs{i} '_hrf.mat'])); [peak(i), idx] = max(S.hrf(:,1)); % 峰值幅度与索引 t2p(i) = (idx-1) * S.P.dt; % 索引换算成秒 end % 一眼看分布，离群值直接标出来 figure; histogram(t2p, 20); xlabel('time to peak (s)'); ylabel('count'); outlier = subs(t2p > 8 | t2p < 3); disp(outlier);正常成人的 HRF 峰值时间大概落在 4 到 7 秒这个区间，中位数在 5 到 6 秒。如果一批数据里有被试的峰值时间小于 3 秒或者大于 9 秒，先别急着当发现，多半是数据问题：头动过大、TR 填错、掩模里混进了脑室或者白质。把这些被试挑出来单独看原始时间序列，能省掉后面统计阶段一堆解释不清的异常值。
半高宽也是个有用的指标，反映 HRF 的“宽窄”。太窄的曲线往往意味着去卷积过拟合了噪声，太宽则可能是平滑过度或者信号本身很弱。两个指标一起看，比只看峰值可靠。
5. 有效连接分析：去卷积之后的第二条流水线
5.1 为什么必须先反卷积再做连接
BOLD 信号是神经活动经过 HRF 卷积之后的结果，这个卷积操作在时间上做了平滑，会把原本存在的时序关系糊掉，甚至制造出方向相反的假象。同一个 TR 内发生的两个事件，在 BOLD 上可能表现为同时变化；而一个区域领先另一个区域几百毫秒的真实因果，在 BOLD 上可能被 HRF 的延迟吞掉，看都看不见。
所以做有效连接之前先反卷积，逻辑是站得住的：先把每个节点的 BOLD 还原成神经事件序列，再在这些事件序列上估计时序依赖。这一步做完，采样率也从 TR 提升到了dt的栅格，对 Granger 这类对时间分辨率敏感的模型是实打实的改善。
但要清醒：反卷积不是无损还原。HRF 估计本身有误差，噪声也被一并放大，去卷积之后的序列比原始 BOLD 更“毛”。所以连接分析的结论，得结合 HRF 拟合质量一起看——HRF 都估得不像样的节点，它的连接结果没什么参考价值。
5.2 一个最小可用的 MVAR-Granger 实现
工具箱里带了连接分析的函数，但接口和参数结构跟你自己的流水线未必对得上。我一般自己写一个几十行的 MVAR-Granger，好处是完全可控，坏处是要自己保证正确性。下面这段在多个数据集上验证过，可以直接抄。
function [F, pval] = mvar_gc(x, y, order) % x, y : 去卷积后的神经事件序列，列向量，等长 % order: MVAR 模型阶数，静息态常用 3~5 % F : 从 x 到 y 的 Granger F 值 % pval : 对应 F 检验的 p 值 x = x(:); y = y(:); T = numel(y); n = T - order; % 构造滞后设计矩阵：第 k 组两列 = [x 的 k 阶滞后, y 的 k 阶滞后] Z = zeros(n, 2*order); for k = 1:order Z(:, 2*k-1) = x(order+1-k : end-k); Z(:, 2*k ) = y(order+1-k : end-k); end target = y(order+1:end); % 受限模型：只用自己的滞后预测自己 Zr = Z(:, 2:2:end); br = Zr \ target; rss_r = sum((target - Zr*br).^2); % 全模型：加入 x 的滞后 bf = Z \ target; rss_f = sum((target - Z*bf).^2); % F 检验，分子自由度 order，分母自由度 n - 2*order - 1 F = ((rss_r - rss_f) / order) / (rss_f / (n - 2*order - 1)); pval = 1 - fcdf(F, order, n - 2*order - 1); end索引那几行是整段的核心。Z(:, 2*k-1)取的是 x 的第 k 阶滞后，Z(:, 2*k)取的是 y 的第 k 阶滞后，所以2:2:end正好是全模型里所有属于 y 自己的列，也就是受限模型的设计矩阵。这个排列方式一旦写错，受限模型和全模型用的就不是同一批样本，F 值直接失去意义。
自由度也别写错。分母自由度是样本量减去全模型参数个数2*order再减一，写少了 p 值会系统性偏小，一堆本该不显著的结果变成显著。
调用的时候，把每个节点的event_bold列拿出来两两配对：
order = 3; N = size(event_bold, 2); GC = zeros(N, N); for i = 1:N for j = 1:N if i ~= j [GC(i,j), ~] = mvar_gc(event_bold(:,i), event_bold(:,j), order); end end endGC(i,j)的含义是节点 i 到节点 j 的因果强度。注意方向，别读反了，这个错误在组会汇报里出现过不止一次。
5.3 去卷积信号进连接模型前的三个处理
直接拿event_bold喂进 MVAR 通常会得到一堆虚高的 F 值，原因是反卷积把噪声也放大了，高频成分尤其明显。有三个处理我基本每次都会做。
第一，平滑。对去卷积后的序列做一次轻微的高斯平滑，窗口大概 2 到 3 个dt，目的是压掉反卷积引入的高频抖动。平滑太狠会把真实的时序关系也抹掉，所以要克制。
第二，降采样。dt往往是 0.2 秒，而 MVAR 的阶数乘上采样间隔要能覆盖真实的生理延迟。在 0.2 秒栅格上用 3 阶模型，覆盖的时间跨度只有 0.6 秒，对 fMRI 这种以秒为单位的信号来说太短了。常见的做法是降采样到 1 秒左右，再用 3 到 5 阶。
第三，去均值、去趋势。MVAR 假设序列是平稳的，残留的线性漂移会让模型把趋势误判成因果。做一次线性去趋势，成本很低。
提示：这三步处理会显著改变连接结果的绝对数值，所以在报告结果时要把处理流程写清楚。跨研究比较 Granger 的绝对值基本没有意义，能比的是方向和相对强弱，以及在同一套流程下组间的差异。
6. 常见问题排查与二次开发：把它改成你流水线里的一块
6.1 三个高频翻车点
| 现象 | 原因 | 解决 |
|---|---|---|
报错指向一个没调用过的spm_*函数 | 路径顺序反了，旧版同名函数遮蔽了 SPM12 | 先addpathSPM12，再spm('defaults','FMRI')，最后加 rsHRF |
| HRF 曲线几乎是一条平线 | cutoff填成了赫兹值，滤波把信号切光了 | cutoff单位是秒，静息态常用 128 |
| 峰值时间集中在 2 到 3 秒且被试间几乎无差异 | TR 填错或P.T长度与数据行数不一致 | 核对 NIfTI 头里的 TR，检查numel(P.T) == size(y,1) |
第一条是最容易让人怀疑人生的，因为报错信息完全不指向真正的原因。第二条我见过两次，都是直接照搬别人的脚本没看单位。第三条最隐蔽，因为结果看起来“很一致”，反而容易被当成好现象。
6.2 二次开发的三个改造点
想把它接进自己的流水线，有三处值得动手。第一处是把rsHRF_est的输入输出包一层，让数据读写走你自己的 BIDS 组织逻辑，而不是依赖工具箱的文件扫描；这样换数据集不用改算法代码。第二处是把P的构造抽成配置文件，被试级的 TR 和节点数从表里读，避免在脚本里硬编码。第三处是在估计之后加一个自动质控钩子，峰值时间越界或者拟合残差超阈值的被试直接打标记进日志，批处理跑完先看日志再决定下一步。
改完之后一定要做回归验证：拿一小批已经跑通的数据，新旧两条路径各跑一遍，比对beta_hrf的差异。峰值时间的偏差应该小于dt，HRF 曲线的相关系数应该在 0.95 以上。达不到就说明改造过程中动了不该动的东西，老老实实回退到一个能对上的版本再逐个改动点排查。
从那以后我每次改这类工具箱的调用封装，都强制先跑一遍小样本回归对比，曲线对不上绝不往下走。希望帮到你。
本文还有配套的精品资源，点击获取
