简介：本资源是面向静息态功能磁共振研究者的rsHRF工具箱，基于MATLAB实现静态HRF反卷积与连通性分析，支持独立运行及作为SPM插件使用，适合具备一定MATLAB与神经影像基础、需要估计和去除静息态HRF变异影响的研究人员。压缩包共148个文件，约5.82MB，以130个.m函数脚本为核心，辅以.mat数据、.txt与.md说明文档、.zip依赖包、.pdf教程及.nii、.gii影像示例等，覆盖配置、算法与演示数据。资源内含demo_codes演示脚本与示例数据集，并附分步教程，可帮助读者快速跑通HRF估计、反卷积与后续连通性分析流程，理解工具箱的目录结构与调用方式。目前已有455人学习下载，适合作为静息态HRF分析入门的实操参考。

1. 静态HRF反卷积到底在解决什么问题：从BOLD信号到神经活动的反推

做过静息态fMRI分析的人大概都遇到过这个尴尬：你手里有一堆BOLD时间序列，做功能连接、做脑网络、算图论指标，结果审稿人一句“你考虑过HRF变异性对连接结果的影响吗”就能把你问住。BOLD信号本身不是神经活动，它是神经活动经过血流动力学响应函数（HRF）卷积之后的结果。不同脑区、不同被试的HRF在峰值时间和宽度上都有差异，直接拿BOLD信号算相关，等于把HRF的差异混进了连接强度里。rsHRF这个MATLAB工具箱要干的事情，就是把这个卷积过程反过来——从观测到的BOLD信号中估计出每个体素或每个脑区的HRF，然后做反卷积，还原出潜在的神经活动时间序列，再基于还原后的信号做连通性分析。它适合两类人：一是做静息态连接分析、对HRF变异性敏感的研究者；二是想从BOLD信号里提取HRF参数（峰值时间、峰值高度、宽度）做个体差异分析的从业者。工具箱本身是开源的，代码结构清晰，核心函数都有注释，二次开发的成本不算高。

2. rsHRF工具箱的核心方法：参数化与非参数化反卷积怎么选

2.1 两种反卷积策略的数学逻辑

rsHRF提供了两条技术路线。第一条是参数化方法，假设HRF可以用一个双伽马函数（或它的变体）来描述，通过非线性最小二乘拟合来估计每个体素的HRF参数。双伽马函数的典型形式是两个伽马函数的差：

h(t) = A * (t^(a1-1) * exp(-t/b1) / (b1^a1 * gamma(a1)) - c * t^(a2-1) * exp(-t/b2) / (b2^a2 * gamma(a2)))

其中A控制幅度，a1、b1决定峰值形态，a2、b2决定 undershoot 的形态，c是 undershoot 相对于峰值的比例。参数化方法的优势是输出可解释的HRF参数，每个体素都能拿到峰值时间、峰值高度、FWHM（半高全宽），方便做组间统计。缺点是如果真实HRF偏离双伽马假设（比如被试有血管病变、或者脑区靠近大血管），拟合会强行套一个不合适的形状，反卷积出来的神经信号反而失真。

第二条是非参数化方法，不预设HRF形状，直接用正则化的最小二乘或维纳滤波来估计HRF。rsHRF里实现的是基于样条基函数的非参数估计，通过控制平滑度参数来平衡拟合优度和HRF的平滑性。非参数方法的优势是灵活，能捕捉到参数化方法漏掉的HRF形态变异；代价是计算量大，而且平滑度参数的选择对结果影响明显——太小平滑不够，HRF会出现振荡；太大则过度平滑，峰值被抹掉。

2.2 参数选择与实操建议

我一般会先跑参数化方法，把每个体素的HRF参数导出来看一眼分布。如果峰值时间的分布集中在5-7秒、FWHM在4-6秒，说明数据质量不错，参数化结果可信。如果峰值时间出现大量小于3秒或大于10秒的体素，要么是数据质量有问题（头动、大血管），要么是这些体素的HRF确实偏离标准形状，这时候再考虑非参数化方法。

参数化方法的关键参数在rsHRF_estimation.m里：

% 参数化HRF估计的核心参数设置 options.TR = 2; % 重复时间，单位秒，必须和你的扫描参数一致 options.T = 200; % 时间点数，根据实际数据长度设置 options.T0 = 1; % 起始时间偏移，通常设为1 options.TD = 1; % 时间导数标志，1表示同时估计HRF的时间导数 options.AR = 1; % 自回归模型阶数，用于预白化，通常设为1 options.estimation = 'parametric'; % 选择参数化方法 options.hrf = 'canonical'; % 初始HRF形状，可选'canonical'或'boxcar'

TR必须和实际扫描参数一致，这个参数错了后面全错。TD设为1时会多估计一个时间导数，相当于允许HRF在时间轴上轻微平移，对HRF延迟变异大的数据有帮助，但也会增加计算量。AR阶数用于处理fMRI时间序列的自相关，标准做法是设1或2，设太高会过度校正。

非参数化方法的调用方式不同，需要指定基函数数量和平滑参数：

% 非参数化HRF估计 options.estimation = 'nonparametric'; options.nbasis = 10; % 样条基函数数量，通常8-12 options.lambda = 0.1; % 平滑参数，越大越平滑 options.dt = 0.1; % 时间分辨率，用于HRF采样

nbasis太少（比如小于6）会导致HRF估计过于粗糙，峰值可能被抹平；太多（大于15）则容易过拟合，HRF出现不合理的振荡。lambda的选择更玄学，我一般从0.01、0.1、1三个值里试，看哪个跑出来的HRF形态最合理——所谓合理，就是峰值在5-7秒、有明确的undershoot、没有明显的负向振荡。

2.3 反卷积与连通性分析的衔接

估计出HRF之后，反卷积步骤在rsHRF_deconvolution.m里完成。核心思路是用估计的HRF构造卷积矩阵，然后通过正则化最小二乘求解神经活动：

% 反卷积得到神经活动时间序列 [beta, resid] = rsHRF_deconvolution(Y, HRF, options); % Y是观测的BOLD信号，HRF是估计的血流动力学响应函数 % beta就是还原后的神经活动时间序列

拿到beta之后，连通性分析就和常规流程一样了——算相关、算偏相关、做Fisher z变换。但要注意，反卷积后的信号自相关结构会发生变化，做统计检验时自由度需要重新估计。rsHRF提供了一个rsHRF_connectivity.m函数来处理这个问题，内部用有效自由度来校正：

% 基于反卷积信号的连通性分析 conn = rsHRF_connectivity(beta, options); % conn.r是相关系数矩阵，conn.p是校正后的p值 % options.dof_adjust = 1 开启有效自由度校正

有效自由度校正的逻辑是：反卷积相当于对信号做了滤波，滤波后的信号自相关增强，实际独立样本数减少。如果不校正，p值会偏小，假阳性率升高。这个坑我在早期项目里踩过，审稿人直接指出连接显著性被高估，后来老老实实把dof_adjust打开才通过。

3. 从BOLD到神经活动：完整实操流程与参数配置

3.1 数据准备与预处理要求

rsHRF对输入数据的要求不算苛刻，但有几个硬性条件。第一，数据必须是预处理后的，至少完成时间层校正、头动校正、空间标准化、平滑。第二，需要去趋势和滤波，但滤波方式有讲究——rsHRF官方建议不要做带通滤波，因为滤波会改变HRF的形态，影响估计精度。如果一定要滤波，建议在反卷积之后对神经活动信号做滤波，而不是在BOLD信号上做。

第三，需要回归掉头动参数、白质信号、脑脊液信号。这些 nuisance regressor 的回归要在HRF估计之前完成，否则头动相关的信号波动会被误认为是神经活动。我一般用spm_regress或者自己写脚本做这一步：

% 加载预处理后的BOLD数据 load('preprocessed_bold.mat'); % 假设变量名为Y，维度是T x V % Y的每一列是一个体素的时间序列，T是时间点数 % 回归 nuisance regressor nuisance = [head_motion, white_matter, csf]; % 头动6参数+WM+CSF Y_regressed = Y - nuisance * (nuisance \ Y); % 去线性趋势 Y_detrended = detrend(Y_regressed); % 注意：不要在这里做带通滤波

头动参数用6个（3个平移+3个旋转）还是24个（加上导数）取决于头动大小。如果被试头动小于0.5mm，6个够用；如果头动在0.5-1mm之间，建议用24参数模型。白质和脑脊液信号从各自的mask里提取，取平均时间序列。

3.2 逐体素HRF估计与反卷积

数据准备好之后，就可以跑rsHRF的主流程了。我一般写一个循环，逐体素处理，这样内存占用小，也方便并行：

% 逐体素HRF估计与反卷积 nVoxels = size(Y_detrended, 2); HRF_params = zeros(nVoxels, 4); % 存储峰值时间、峰值高度、FWHM、拟合优度 neural_signal = zeros(size(Y_detrended)); for v = 1:nVoxels y = Y_detrended(:, v); % 跳过方差过小的体素（通常是背景或脑室） if std(y) < 0.1 * mean(std(Y_detrended)) continue; end % 参数化HRF估计 options.TR = 2; options.T = length(y); options.estimation = 'parametric'; options.hrf = 'canonical'; options.TD = 1; options.AR = 1; [hrf, params] = rsHRF_estimation(y, options); % 保存HRF参数 HRF_params(v, :) = [params.peak_time, params.peak_height, ... params.fwhm, params.r2]; % 反卷积 [beta, ~] = rsHRF_deconvolution(y, hrf, options); neural_signal(:, v) = beta; end % 保存结果 save('hrf_results.mat', 'HRF_params', 'neural_signal', '-v7.3');

这段代码有几个细节值得说。第一，跳过方差过小的体素是为了省时间，背景体素的HRF估计没有意义，强行跑还会拖慢整体速度。第二，options.T必须和实际时间序列长度一致，不一致会报错。第三，params.r2是拟合优度，我一般会设一个阈值（比如0.3），低于这个值的体素标记为HRF估计不可靠，后续分析时剔除。

并行化可以用parfor替换for，但要注意每个worker的内存占用。如果体素数量超过10万，建议分块处理，每块5万个体素，跑完一块保存一块。

3.3 连通性计算与统计检验

神经信号拿到之后，连通性计算就是标准流程了。但前面提过，反卷积后的信号自相关结构变了，需要做有效自由度校正：

% 定义ROI并提取平均神经信号 roi_indices = load('roi_mask.mat'); % ROI体素索引 nROI = length(roi_indices); roi_signals = zeros(size(neural_signal, 1), nROI); for r = 1:nROI roi_signals(:, r) = mean(neural_signal(:, roi_indices{r}), 2); end % 计算连通性矩阵 conn_matrix = corr(roi_signals); % 有效自由度校正 T = size(roi_signals, 1); TR = 2; % 反卷积后的信号自相关更强，有效自由度约为T/2到T/3 % 具体值取决于HRF估计的平滑程度 dof_adjust = 2.5; % 保守估计，实际可用rsHRF_connectivity自动计算 t_values = conn_matrix .* sqrt((T/dof_adjust - 2) ./ (1 - conn_matrix.^2)); p_values = 2 * tcdf(-abs(t_values), T/dof_adjust - 2); % FDR校正 [p_fdr, ~] = fdr_bh(p_values(:), 0.05); p_fdr_matrix = reshape(p_fdr, size(p_values));

dof_adjust这个值我一般设2.5，意思是有效自由度是原始时间点数的1/2.5。这个值不是拍脑袋来的——反卷积相当于对信号做了低通滤波，滤波后的信号相邻时间点相关性增强，独立样本数减少。rsHRF的rsHRF_connectivity.m内部用AR模型估计自相关长度，然后算有效自由度，比手动设一个固定值更靠谱。如果不想自己算，直接用那个函数就行。

4. 避坑与排查：HRF估计翻车的五种典型情况

4.1 峰值时间出现大量异常值

现象：跑完参数化估计后，峰值时间的直方图出现双峰或者长尾，一部分体素的峰值时间小于3秒或大于10秒。

原因：小于3秒的体素通常位于大血管附近，BOLD信号受血流影响大，HRF形态偏离标准形状；大于10秒的体素可能是白质或者脑室，这些区域本来就没有明显的神经活动，拟合出来的HRF没有生理意义。另外，如果头动校正没做好，头动引起的信号跳变会被误认为是HRF峰值。

解决：第一，检查头动参数，如果某个被试的最大头动超过2mm，建议剔除或者用更严格的校正方法。第二，用灰质mask限制分析范围，白质和脑室的体素直接排除。第三，对于峰值时间异常的体素，可以尝试非参数化方法，看是否能得到更合理的HRF形态。如果非参数化也救不回来，标记为不可靠体素，后续分析剔除。

4.2 反卷积后的神经信号出现负值

现象：反卷积得到的神经活动时间序列有大量负值，甚至负值的幅度和正值差不多。

原因：反卷积是一个病态逆问题，观测信号里的噪声会被放大。如果HRF估计不准，或者BOLD信号的信噪比太低，反卷积结果就会出现不合理的负值。另外，如果nuisance regressor回归不充分，残留的头动或生理噪声也会导致负值。

解决：第一，检查HRF估计的拟合优度，r2低于0.3的体素反卷积结果不可信。第二，增加正则化强度，在rsHRF_deconvolution.m里调大lambda参数。第三，对反卷积后的信号做非负约束，rsHRF提供了一个options.nonneg = 1的选项，强制神经活动非负。但要注意，非负约束会引入偏差，只在负值特别严重时使用。

4.3 连通性结果与BOLD连通性差异过大

现象：基于反卷积神经信号的连通性矩阵和基于原始BOLD信号的连通性矩阵差异很大，甚至出现符号相反的情况。

原因：反卷积改变了信号的频谱特性，如果HRF估计不准，反卷积相当于对信号做了错误的滤波。另外，如果不同脑区的HRF差异很大，反卷积后神经信号的幅度和相位都会变化，连通性结果自然不同。

解决：第一，对比HRF参数的空间分布，如果某些脑区的峰值时间明显偏离全脑均值，检查这些脑区的数据质量。第二，做一步验证：用模拟数据（已知神经活动+已知HRF）跑一遍rsHRF，看反卷积能否还原出已知的神经活动。如果模拟数据都还原不了，说明参数设置有问题。第三，如果差异确实存在且无法消除，在论文里报告两种结果，说明HRF变异性对连通性的影响。

4.4 计算速度过慢

现象：逐体素跑完整个大脑需要几个小时甚至过夜。

原因：参数化估计的每次拟合都需要迭代优化，体素数量多的时候计算量很大。如果开了TD（时间导数估计），计算量翻倍。另外，如果用的是非参数化方法，基函数数量多的时候矩阵运算也很耗时。

解决：第一，用parfor并行，开4-8个worker，速度能提升3-5倍。第二，先跑一个子集（比如随机选5000个体素），确认参数设置没问题再跑全脑。第三，如果只是想做连通性分析，不需要每个体素的HRF参数，可以先用标准HRF做反卷积，跳过HRF估计步骤，速度会快很多。第四，非参数化方法如果只是为了得到神经信号，可以把nbasis降到6-8，lambda适当调大，牺牲一点灵活性换速度。

4.5 不同被试的HRF参数不可比

现象：组间比较时，A组的峰值时间显著大于B组，但看原始BOLD信号，两组的差异并不明显。

原因：HRF参数估计受数据质量影响很大。如果A组的头动比B组大，或者A组的TR比B组长，估计出来的HRF参数就会有系统偏差。另外，如果两组的扫描参数不同（比如TR不同），HRF估计的基线也不一样。

解决：第一，确保所有被试的扫描参数一致，TR不同的话需要重采样到同一时间分辨率。第二，在组间比较时，把头动参数作为协变量回归掉。第三，用非参数化方法交叉验证，如果参数化方法显示组间差异但非参数化方法没有，说明差异可能来自HRF形状假设的偏差，而不是真实的生理差异。

5. 进阶技巧：用模拟数据验证反卷积精度与HRF参数的可重复性

5.1 构造模拟数据集

验证rsHRF的反卷积精度，最直接的办法是构造模拟数据。先生成一个已知的神经活动时间序列，用一个已知的HRF去卷积它，得到模拟的BOLD信号，然后加噪声，再跑rsHRF看能不能还原出原始的神经活动和HRF。

% 构造模拟数据验证rsHRF TR = 2; T = 200; t = (0:T-1) * TR; % 生成已知的神经活动（稀疏的脉冲序列） neural_true = zeros(T, 1); neural_true(20:25) = 1; % 第一个事件 neural_true(80:85) = 0.8; % 第二个事件 neural_true(140:145) = 1.2; % 第三个事件 % 生成已知的HRF（双伽马） hrf_true = spm_hrf(TR, [6 16 1 1 6 0 32]); % SPM标准HRF % 卷积得到模拟BOLD信号 bold_clean = conv(neural_true, hrf_true); bold_clean = bold_clean(1:T); % 加噪声（SNR约为5） noise = randn(T, 1) * std(bold_clean) / 5; bold_noisy = bold_clean + noise; % 跑rsHRF估计 options.TR = TR; options.T = T; options.estimation = 'parametric'; options.hrf = 'canonical'; options.TD = 1; options.AR = 1; [hrf_est, params] = rsHRF_estimation(bold_noisy, options); [neural_est, ~] = rsHRF_deconvolution(bold_noisy, hrf_est, options); % 计算还原精度 corr_hrf = corr(hrf_true(1:30), hrf_est(1:30)); corr_neural = corr(neural_true, neural_est); fprintf('HRF还原相关系数: %.3f\n', corr_hrf); fprintf('神经活动还原相关系数: %.3f\n', corr_neural);

这段代码的关键在于spm_hrf函数，它来自SPM工具箱，生成标准的双伽马HRF。如果没装SPM，可以自己写一个双伽马函数。噪声水平用SNR控制，SNR=5相当于中等噪声水平，接近真实fMRI数据的信噪比。跑完之后看两个相关系数：HRF还原相关系数应该在0.9以上，神经活动还原相关系数在0.7以上算合格。如果低于这个值，检查参数设置，特别是AR阶数和TD标志。

5.2 参数可重复性评估

模拟数据只能验证方法在理想条件下的精度，真实数据的可重复性需要用分半信度来评估。把时间序列分成前后两半，分别跑rsHRF，比较两次估计的HRF参数：

% 分半信度评估 half1 = Y_detrended(1:round(T/2), :); half2 = Y_detrended(round(T/2)+1:end, :); % 分别估计 params1 = zeros(size(half1, 2), 4); params2 = zeros(size(half2, 2), 4); for v = 1:size(half1, 2) [~, p1] = rsHRF_estimation(half1(:, v), options); [~, p2] = rsHRF_estimation(half2(:, v), options); params1(v, :) = [p1.peak_time, p1.peak_height, p1.fwhm, p1.r2]; params2(v, :) = [p2.peak_time, p2.peak_height, p2.fwhm, p2.r2]; end % 计算组内相关系数（ICC） icc_peak_time = ICC([params1(:,1), params2(:,1)], 'A-1'); icc_fwhm = ICC([params1(:,3), params2(:,3)], 'A-1'); fprintf('峰值时间ICC: %.3f\n', icc_peak_time); fprintf('FWHM ICC: %.3f\n', icc_fwhm);

ICC的计算可以用MATLAB的ICC函数（需要统计工具箱），或者自己写。峰值时间的ICC一般要求在0.6以上，FWHM的ICC在0.5以上。如果ICC太低，说明HRF估计的稳定性不够，可能是时间序列太短（少于150个时间点），或者数据质量太差。我一般会在正式分析前跑一遍分半信度，ICC不达标的话，要么增加扫描时间，要么换非参数化方法试试。

5.3 一个容易忽略的细节：HRF估计的时间窗口

rsHRF默认的HRF估计窗口是32秒，这个窗口覆盖了HRF的峰值和大部分undershoot。但如果TR比较大（比如3秒），32秒只有10个时间点，HRF估计的精度会下降。这时候需要把窗口拉长到40秒甚至48秒。在rsHRF_estimation.m里，窗口长度由options.T和options.TR共同决定，但实际估计时用的是options.hrf_length参数：

% 调整HRF估计窗口 options.hrf_length = 40; % 单位秒，默认32

窗口太短，undershoot估计不准；窗口太长，计算量增加，而且如果时间序列本身不长，窗口太长会截断有效数据。我一般设40秒，兼顾精度和计算量。如果TR小于1.5秒，32秒够用；如果TR大于2.5秒，建议设40-48秒。

从那以后我每次跑rsHRF之前，都会先用模拟数据验证一遍参数设置，再跑分半信度确认ICC达标，最后才跑全脑。这个习惯帮我省了很多返工的时间——有一次没做模拟验证，直接跑全脑，结果反卷积出来的神经信号全是噪声，查了半天才发现是AR阶数设成了0，预白化没做，自相关把HRF估计带偏了。希望这些经验能帮到你，少走点弯路。

本文还有配套的精品资源，点击获取