企业官网建设流程全解析

近红外数据分析在认知神经科学、心理学和医学研究中的应用越来越广泛，但很多研究者和学生在入门时，常常被复杂的信号处理流程、晦涩的算法原理和繁琐的绘图步骤所困扰。网上资料要么过于理论化，要么代码片段零散不成体系，导致从原始数据到可发表图表之间，往往需要耗费大量时间摸索和排错。

本文旨在整合一套从零开始的近红外数据分析实战闭环方案。我们将以最常用的 MATLAB 环境为例，结合经典的 Homer2 和 NIRS-KIT 工具箱，系统性地讲解数据预处理、脑功能成像计算、统计分析和结果可视化的全流程。无论你是心理学、神经科学专业的学生，还是希望将近红外技术应用于临床或工程研究的开发者，都能通过本文的步骤和完整代码，快速上手，避开那些常见的“坑”。

1. 近红外数据分析核心概念与背景

在深入代码之前，我们必须理解我们在处理什么，以及为什么要这样处理。这能帮助你在后续步骤中做出正确的判断，而不是机械地执行操作。

1.1 什么是功能性近红外光谱技术？

功能性近红外光谱技术是一种利用近红外光（波长通常为650-950nm）穿透生物组织（如头皮和颅骨）来检测大脑皮层血红蛋白浓度变化的光学成像技术。它主要基于“神经血管耦合”原理：当大脑某个区域神经元活动增强时，该区域的局部血流量和血氧水平会随之发生变化。

fNIRS设备通过发射器向头皮发射近红外光，并由探测器接收经组织散射后的光信号。通过测量不同波长光强的衰减，我们可以推算出氧合血红蛋白和脱氧血红蛋白浓度的相对变化，从而间接反映神经活动。

与其它脑成像技术的对比：

• 功能磁共振成像：空间分辨率高，但时间分辨率较低（秒级），设备昂贵、笨重，对被试活动限制大。
• 脑电图/事件相关电位：时间分辨率极高（毫秒级），但空间分辨率很差，难以精确定位活动脑区。
• fNIRS：在空间分辨率（厘米级）和时间分辨率（可达0.1秒）之间取得了较好的平衡，设备相对便携，允许被试进行一定范围的自然活动，非常适合研究自然情境下的认知过程、儿童发育、康复评估等。

1.2 fNIRS 数据分析的基本流程

一个标准的 fNIRS数据分析流程可以概括为以下步骤，这也是本文实战部分的结构框架：

1. 数据采集与格式转换：从设备导出原始光强数据，并转换为分析软件（如Homer2）可识别的格式。
2. 数据预处理：这是最关键且最易出错的环节，目的是去除噪声，提取有效的血红蛋白信号。包括：
   • 检测并标记运动伪迹。
   • 将原始光强转换为光密度变化。
   • 进行带通滤波，去除心跳、呼吸等生理噪声和低频漂移。
   • 利用修正的比尔-朗伯定律，将光密度变化转换为血红蛋白浓度变化。
3. 个体水平分析：对单个被试的数据进行建模，计算任务条件与基线条件相比的脑激活。
   • 定义实验范式（事件 onset 和 duration）。
   • 构建一般线性模型，估计每个通道的β值（激活强度）。
4. 组水平分析：将多个被试的数据进行整合，进行群体统计推断。
   • 将个体β值配准到标准脑空间（如MNI空间）。
   • 进行单样本t检验、配对t检验或方差分析等。
5. 结果可视化：生成可发表的统计地图和图表。
   • 绘制脑激活地形图。
   • 绘制时间序列曲线。
   • 绘制条形图和误差线。

理解这个流程后，我们就能明白每一步代码的目的，从而在出现问题时能够有效排查。

2. 环境准备与工具箱安装

工欲善其事，必先利其器。我们将搭建一个基于 MATLAB 的核心分析环境。

2.1 软件与工具箱版本说明

• MATLAB: 推荐使用 R2018b 或更高版本。本文示例基于 R2021a，但核心函数在较新版本中兼容性良好。
• Homer2: 这是最经典、使用最广泛的 fNIRS 预处理工具箱。我们将使用其稳定版本。
• NIRS-KIT: 一个功能强大的国产工具箱，集成了预处理、个体/组分析、统计和绘图的完整流程，对中文用户非常友好，且与 Homer2 数据格式兼容。
• SPM12: 一个通用的神经影像统计分析包，NIRS-KIT 在组分析时会调用其部分功能。

重要提示：不同工具箱的版本可能存在函数接口差异。本文提供的代码和配置思路是通用的，但如果你遇到函数未定义的错误，请首先检查工具箱的文档和函数名。不要盲目复制代码。

2.2 安装步骤详解

假设你的 MATLAB 安装路径为C:\MATLAB，我们建议在用户目录下（如D:\fNIRS_Analysis）创建工作文件夹。

步骤1：下载工具箱

1. 访问 Homer2 官网或 GitHub 仓库，下载最新稳定版，解压到D:\fNIRS_Analysis\Toolboxes\homer2。
2. 访问 NIRS-KIT 的 GitHub 仓库，下载完整工具箱，解压到D:\fNIRS_Analysis\Toolboxes\nirs_kit。
3. 访问 SPM 官网，下载 SPM12，解压到D:\fNIRS_Analysis\Toolboxes\spm12。

步骤2：设置 MATLAB 路径这是最关键的一步，路径设置错误会导致工具箱函数无法调用。 打开 MATLAB，在“主页”选项卡点击“设置路径”。在弹出的对话框中：

1. 点击“添加并包含子文件夹”。
2. 依次添加以下三个文件夹：
   • D:\fNIRS_Analysis\Toolboxes\homer2
   • D:\fNIRS_Analysis\Toolboxes\nirs_kit
   • D:\fNIRS_Analysis\Toolboxes\spm12
3. 点击“保存”，然后关闭对话框。务必保存，否则下次启动 MATLAB 需要重新添加。

步骤3：验证安装在 MATLAB 命令窗口中分别输入以下命令，不报错即表示路径设置成功。

% 验证 Homer2 which hmrR_Intensity2OD % 应返回类似：D:\fNIRS_Analysis\Toolboxes\homer2\hmrR_Intensity2OD.m % 验证 NIRS-KIT (以其中一个函数为例) which NIRS_KIT_version % 应返回类似：D:\fNIRS_Analysis\Toolboxes\nirs_kit\NIRS_KIT_version.m

3. 数据预处理实战：从原始光强到干净信号

我们将使用一个模拟的示例数据来演示完整流程。假设你的原始数据文件为sub-01_task-motor_raw.nirs，这是一个包含光强、刺激标记、采样频率等信息的标准.nirs文件。

3.1 数据加载与初步检查

首先，我们需要将数据读入 MATLAB 工作空间，并查看其基本结构。

% 文件路径：preprocess_pipeline.m clear; clc; close all; % 1. 加载数据 data_path = 'D:\fNIRS_Analysis\Data\sub-01\'; filename = 'sub-01_task-motor_raw.nirs'; load(fullfile(data_path, filename), '-mat'); % 加载后变量名为 'raw_data' % 2. 查看数据结构 disp('=== 数据结构信息 ==='); whos raw_data % 通常包含以下字段： % d: 原始光强数据矩阵 [时间点 x 通道] % s: 刺激标记矩阵 [时间点 x 条件数] % t: 时间轴向量 [时间点 x 1] % aux: 辅助信号（如加速度计数据） % SD: 探头结构体，包含光源、探测器的位置和配对信息 % 3. 查看关键参数 fs = 1 / mean(diff(raw_data.t)); % 计算采样频率 disp(['采样频率: ', num2str(fs), ' Hz']); disp(['数据长度: ', num2str(length(raw_data.t)), ' 个时间点']); disp(['通道数量: ', num2str(size(raw_data.d, 2))]);

运行这段代码，你可以确认数据是否被正确加载，并了解其基本维度，这是后续所有处理的基础。

3.2 核心预处理步骤（基于Homer2函数）

预处理的目标是去除噪声，保留与任务相关的血红蛋白信号变化。我们按照标准流程进行。

% 文件路径：preprocess_pipeline.m (续) % 4. 将原始光强转换为光密度变化 % 这是应用比尔-朗伯定律的第一步 od_data = hmrR_Intensity2OD(raw_data.d); % 5. 检测运动伪迹 (使用经典的 tMotion 和 tMask 方法) % tMotion: 判断信号变化是否超过阈值的标准差倍数 % tMask: 确定需要插值的时段长度 tMotion = 1.0; % 推荐值 0.5-1.0 tMask = 1.0; % 推荐值 0.5-1.5 (秒) [od_data_corrected, ~] = hmrR_MotionCorrectPCArecurse(od_data, raw_data.t, raw_data.SD, tMotion, tMask); % 6. 带通滤波 % 去除高频生理噪声（如心跳~1Hz）和低频漂移（如 Mayer波~0.1Hz） hpf = 0.01; % 高通滤波截止频率，去除低频漂移 (单位: Hz) lpf = 0.5; % 低通滤波截止频率，去除高频噪声 (单位: Hz) od_data_filtered = hmrR_BandpassFilt(od_data_corrected, raw_data.t, hpf, lpf); % 7. 将光密度转换为血红蛋白浓度变化 % 使用修正的比尔-朗伯定律，需要消光系数和微分路径因子 ppf = [6.0, 6.0]; % 微分路径因子，通常对两个波长设为相同值 hb_data = hmrR_OD2Conc(od_data_filtered, raw_data.SD, ppf); % hb_data 是一个结构体，通常包含： % hb_data.HbO: 氧合血红蛋白浓度变化 [时间点 x 通道] % hb_data.HbR: 脱氧血红蛋白浓度变化 [时间点 x 通道] % hb_data.HbT: 总血红蛋白浓度变化 (HbO+HbR) disp('预处理完成！');

关键参数解释与避坑指南：

• tMotion 和 tMask：这两个参数对运动伪迹校正效果影响巨大。值设得太小，可能无法检测到真实的运动；值设得太大，可能将生理信号误判为运动而过度校正。务必通过绘制原始和校正后的信号来肉眼检查效果。
• 滤波频率：hpf=0.01意味着保留周期低于100秒（1/0.01）的信号变化，这通常能有效去除缓慢的基线漂移。lpf=0.5意味着去除频率高于0.5Hz的信号，这可以滤除大部分心跳噪声（~1Hz）。根据你的任务设计（如事件间隔）调整这些值。
• 微分路径因子：这是一个经验值，表示光在头皮和大脑之间传播的实际路径长度与光源-探测器距离的比值。成人常用值为6.0。对于婴儿或特殊人群，需要查阅文献使用特定值。

3.3 预处理结果可视化检查

在继续分析前，必须可视化检查预处理效果。这是避免“垃圾进，垃圾出”的关键。

% 文件路径：visualize_preprocess.m % 选择一个代表性通道进行可视化检查 ch_to_plot = 10; % 假设检查第10通道 time = raw_data.t; figure('Position', [100, 100, 1200, 800]); % 子图1：原始光强 vs 光密度 subplot(4,1,1); plot(time, raw_data.d(:, ch_to_plot), 'k'); title(['通道 ', num2str(ch_to_plot), ' - 原始光强']); xlabel('时间 (秒)'); ylabel('光强 (a.u.)'); grid on; subplot(4,1,2); plot(time, od_data(:, ch_to_plot), 'b'); title('转换为光密度变化'); xlabel('时间 (秒)'); ylabel('OD'); grid on; % 子图2：运动校正前后对比 subplot(4,1,3); plot(time, od_data(:, ch_to_plot), 'b'); hold on; plot(time, od_data_corrected(:, ch_to_plot), 'r--', 'LineWidth', 1.5); title('运动校正前后对比 (蓝色:原始, 红色:校正后)'); xlabel('时间 (秒)'); ylabel('OD'); legend('原始', '校正后'); grid on; % 子图3：最终血红蛋白浓度变化 subplot(4,1,4); plot(time, hb_data.HbO(:, ch_to_plot), 'r', 'LineWidth', 1.5); hold on; plot(time, hb_data.HbR(:, ch_to_plot), 'b', 'LineWidth', 1.5); plot(time, hb_data.HbT(:, ch_to_plot), 'g', 'LineWidth', 1.5); title('血红蛋白浓度变化'); xlabel('时间 (秒)'); ylabel('\Delta\mu M'); legend('HbO', 'HbR', 'HbT'); grid on; % 标记刺激事件 % 假设第一个刺激标记在 s 矩阵的第一列 stim_times = find(raw_data.s(:,1) > 0); for i = 1:length(stim_times) x_line = time(stim_times(i)); subplot(4,1,4); line([x_line x_line], ylim, 'Color', 'k', 'LineStyle', '--', 'LineWidth', 0.5); end

通过这个多面板图，你可以清晰地看到数据在每个处理阶段的变化，确认运动伪迹是否被有效抑制，以及血红蛋白信号是否清晰。如果某个通道信号质量极差（如全程饱和或噪声过大），应考虑在后续分析中将其排除。

4. 个体水平分析：构建一般线性模型

预处理后，我们得到了干净的 HbO/HbR 时间序列。接下来，我们要量化每个通道在任务期间的激活程度。最常用的方法是一般线性模型。

4.1 使用 NIRS-KIT 进行便捷的个体分析

NIRS-KIT 封装了 GLM 分析的完整流程，使用起来比手动编写设计矩阵更加方便。

% 文件路径：individual_glm_nirskit.m % 假设你已经完成了预处理，并将结果保存为 hb_data 和 raw_data % 1. 准备输入结构体 IndividualData = struct(); IndividualData.oxy = hb_data.HbO; % 氧合血红蛋白数据 IndividualData.dxy = hb_data.HbR; % 脱氧血红蛋白数据 IndividualData.name = {'sub-01'}; % 被试ID % 注意：NIRS-KIT期望数据维度为 [时间点 x 通道 x 被试] % 目前我们只有一个被试，需要增加一个维度 IndividualData.oxy = permute(IndividualData.oxy, [1, 2, 3]); % 变成 [时间点 x 通道 x 1] IndividualData.dxy = permute(IndividualData.dxy, [1, 2, 3]); % 2. 设置分析参数 para = struct(); para.T1 = 0; % 刺激开始时间 (相对于扫描开始，通常为0) para.T2 = 20; % 刺激结束时间 (根据你的实验范式设定，单位秒) para.interval = [0, 15]; % 分析的时间窗口，例如分析刺激后0-15秒的血流响应 para.base = 5; % 基线时间长度 (秒)，用于计算基线均值 para.fs = fs; % 采样频率 % 3. 定义实验范式 % 假设是一个简单的区块设计，有3个条件（如任务A，任务B，休息） % raw_data.s 矩阵的每一列代表一个条件，值为1表示该时间点该条件发生 condition_names = {'TaskA', 'TaskB', 'Rest'}; % 与 raw_data.s 的列对应 para.condName = condition_names; % 4. 运行个体水平 GLM 分析 % 这个函数会为每个通道、每个条件、每个血红蛋白类型计算β值（激活强度）和t值 [beta, tval, ~, ~] = NIRS_GLM(IndividualData, raw_data.s, para); disp('个体水平GLM分析完成！'); % beta 是一个结构体：beta.oxy 和 beta.dxy，维度为 [通道数 x 条件数 x 被试数] % tval 同理，是统计检验的t值

运行后，beta.oxy(:, 1, 1)就代表了第一个被试在所有通道上，对于“TaskA”条件的 HbO 激活强度估计值。这些 β 值是进行组水平统计的基础。

5. 组水平分析与统计推断

单个被试的结果受个体差异影响很大。我们需要将一组被试的数据放在一起，进行统计检验，判断哪些脑区的激活在群体水平上是显著的。

5.1 准备组水平数据

假设我们已经处理了10个被试的数据，并将每个被试的 β 值保存了下来。现在需要将它们整合到一个数据结构中。

% 文件路径：group_analysis_prep.m % 假设有10个被试，数据保存在一个cell数组或结构数组中 num_subjects = 10; num_channels = size(beta.oxy, 1); % 假设所有被试通道数相同 num_conditions = size(beta.oxy, 2); GroupData_HbO = zeros(num_channels, num_conditions, num_subjects); GroupData_HbR = zeros(num_channels, num_conditions, num_subjects); for sub = 1:num_subjects % 这里需要根据你实际的数据加载方式来编写 % 例如：load(['beta_sub', num2str(sub, '%02d'), '.mat'], 'beta'); % 然后赋值： % GroupData_HbO(:, :, sub) = beta.oxy; % GroupData_HbR(:, :, sub) = beta.dxy; end % 我们以模拟数据为例 % 假设通道5在条件1下，所有被试都有较强的激活（正β值） GroupData_HbO(5, 1, :) = randn(num_subjects, 1) + 2.0; % 均值2，标准差1的正态分布 GroupData_HbR(5, 1, :) = randn(num_subjects, 1) - 0.5; % HbR可能轻微负激活 % 其他通道和条件设为随机噪声 for ch = 1:num_channels for cond = 1:num_conditions if ~(ch==5 && cond==1) GroupData_HbO(ch, cond, :) = randn(num_subjects, 1) * 0.3; GroupData_HbR(ch, cond, :) = randn(num_subjects, 1) * 0.3; end end end

5.2 执行单样本t检验（与0比较）

这是最常用的检验，用于判断某个条件下，某个通道的激活是否显著不为零（即是否存在显著激活）。

% 文件路径：group_ttest.m % 对 HbO，条件1，进行单样本t检验 cond_idx = 1; alpha = 0.05; % 显著性水平 tail = 'both'; % 双尾检验 tvals_HbO = zeros(num_channels, 1); pvals_HbO = zeros(num_channels, 1); for ch = 1:num_channels data = squeeze(GroupData_HbO(ch, cond_idx, :)); % 提取该通道所有被试的β值 [h, p, ci, stats] = ttest(data, 0, 'Alpha', alpha, 'Tail', tail); tvals_HbO(ch) = stats.tstat; pvals_HbO(ch) = p; end % 进行多重比较校正（例如FDR校正） [p_fdr_HbO, ~] = mafdr(pvals_HbO, 'BHFDR', true); % 需要生物信息学工具箱 % 或者使用简单的 Bonferroni 校正 p_corr_HbO = pvals_HbO * num_channels; % Bonferroni校正 p_corr_HbO(p_corr_HbO > 1) = 1; % p值不能大于1 % 找出经过校正后显著的通道 sig_channels_HbO = find(p_fdr_HbO < alpha); disp(['HbO 条件 ', num2str(cond_idx), ' 经过FDR校正后显著的通道有: ', num2str(sig_channels_HbO')]);

对于配对t检验（比较两个条件）或方差分析（比较多个条件），可以使用 MATLAB 的ttest或anova1函数，逻辑类似，但需要组织不同的数据输入格式。

6. 结果可视化：绘制专业图表

得到统计结果后，我们需要将其直观地呈现出来。常见的图表包括脑激活地形图和条件对比条形图。

6.1 绘制脑激活地形图

我们需要探头位置信息（raw_data.SD）来将通道的统计值映射到头皮空间。

% 文件路径：plot_topography.m % 假设我们使用 t 值来绘制地形图 stat_values = tvals_HbO; % 使用上一步计算出的t值 % 或者使用 beta 值的组平均 % stat_values = mean(GroupData_HbO(:, cond_idx, :), 3); % 提取光源和探测器的2D或3D坐标 % 这里假设 SD 结构体中包含 'SrcPos' 和 'DetPos' src_pos = raw_data.SD.SrcPos; det_pos = raw_data.SD.DetPos; % 计算每个通道的中点位置（通常作为通道的坐标） ch_pos = zeros(length(raw_data.SD.MeasList)/2, 3); for ch = 1:size(ch_pos, 1) src_idx = raw_data.SD.MeasList(ch, 1); det_idx = raw_data.SD.MeasList(ch, 2); ch_pos(ch, :) = (src_pos(src_idx, :) + det_pos(det_idx, :)) / 2; end % 绘制2D地形图 (假设坐标已经是2D或我们只取前两维) figure; scatter(ch_pos(:,1), ch_pos(:,2), 200, stat_values, 'filled'); colorbar; colormap(jet); % 可以使用其他colormap，如 parula, hot, coolwarm title(['HbO 条件 ', num2str(cond_idx), ' 激活地形图 (t值)']); xlabel('X (mm or a.u.)'); ylabel('Y (mm or a.u.)'); axis equal; grid on; % 标记显著通道 hold on; if exist('sig_channels_HbO', 'var') scatter(ch_pos(sig_channels_HbO,1), ch_pos(sig_channels_HbO,2), 250, 'k', 'x', 'LineWidth', 2); legend('通道t值', 'FDR显著通道', 'Location', 'best'); end

为了得到更平滑、更美观的“热图”式地形图，通常需要进行空间插值。可以使用griddata或scatteredInterpolant函数。

6.2 绘制时间序列响应曲线

这对于展示血红蛋白信号在任务期间随时间变化的模式非常有用。

% 文件路径：plot_timecourse.m % 选择一个感兴趣的通道和条件 ch_of_interest = 5; cond_of_interest = 1; % 提取该通道所有被试在任务期间的时间序列（需要原始预处理后的数据） % 假设我们有一个cell数组 all_hb_data，包含了所有被试预处理后的 hb_data time_window = para.interval; % 例如 [0, 15] 秒 time_idx = find(raw_data.t >= time_window(1) & raw_data.t <= time_window(2)); time_vector = raw_data.t(time_idx); % 初始化矩阵存储所有被试的时间序列 all_HbO_timeseries = zeros(length(time_vector), num_subjects); all_HbR_timeseries = zeros(length(time_vector), num_subjects); for sub = 1:num_subjects % 加载或访问第 sub 个被试的 hb_data % hb_data_sub = all_hb_data{sub}; % 这里用模拟数据代替 % 模拟一个典型的HRF（血流动力学响应函数）形状 hrf = gampdf(time_vector-2, 6, 0.8); % 峰值在约4-5秒的Gamma函数 hrf = hrf / max(hrf) * 3; % 缩放振幅 all_HbO_timeseries(:, sub) = hrf + randn(size(hrf))*0.2; % 加上噪声 all_HbR_timeseries(:, sub) = -0.3*hrf + randn(size(hrf))*0.1; end % 计算组平均和标准误 mean_HbO = mean(all_HbO_timeseries, 2); sem_HbO = std(all_HbO_timeseries, 0, 2) / sqrt(num_subjects); mean_HbR = mean(all_HbR_timeseries, 2); sem_HbR = std(all_HbR_timeseries, 0, 2) / sqrt(num_subjects); % 绘制带有阴影误差带的曲线 figure; hold on; % HbO fill([time_vector; flipud(time_vector)], ... [mean_HbO - sem_HbO; flipud(mean_HbO + sem_HbO)], ... [1, 0.8, 0.8], 'EdgeColor', 'none', 'FaceAlpha', 0.5); plot(time_vector, mean_HbO, 'r-', 'LineWidth', 2); % HbR fill([time_vector; flipud(time_vector)], ... [mean_HbR - sem_HbR; flipud(mean_HbR + sem_HbR)], ... [0.8, 0.8, 1], 'EdgeColor', 'none', 'FaceAlpha', 0.5); plot(time_vector, mean_HbR, 'b-', 'LineWidth', 2); xlabel('时间 (秒)'); ylabel('\Delta\mu M'); title(['通道 ', num2str(ch_of_interest), ' - 条件 ', num2str(cond_of_interest), ' 的血流响应']); legend('HbO ± SEM', 'HbO 均值', 'HbR ± SEM', 'HbR 均值'); grid on; xlim(time_window);

这张图能清晰展示出典型的 HbO 上升、HbR 下降的血流动力学响应模式，是论文中非常有力的证据。

7. 常见问题与排查思路

在实际操作中，你几乎一定会遇到下面这些问题。这里提供一份排查清单。

8. 最佳实践与工程化建议

将分析流程工程化、规范化，能极大提升研究效率和结果的可重复性。

1. 建立标准化的文件夹结构

   Project_Root/ ├── Data/ │ ├── Raw/ # 存放原始 .nirs/.snirf 文件 │ ├── Preprocessed/ # 存放预处理后的 .mat 文件 │ └── Demographics/ # 被试信息表 ├── Code/ │ ├── 01_Preprocess/ │ ├── 02_Individual_Analysis/ │ ├── 03_Group_Analysis/ │ ├── 04_Visualization/ │ └── utils/ # 自定义函数 ├── Results/ │ ├── Figures/ # 生成的图表 │ └── Stats/ # 统计结果表格 └── Docs/ # 实验协议、分析笔记

   使用脚本自动按照此结构组织数据和结果。

2. 编写可复用的分析脚本而非交互式命令

   • 为每个主要步骤（预处理、个体分析、组分析、绘图）编写独立的.m脚本或函数。
   • 在脚本开头用clear; clc; close all;清空环境，确保结果可重复。
   • 使用明确的变量名，并添加充足的注释。
   • 将关键参数（如滤波频率、运动校正阈值）定义为脚本开头的变量，方便集中修改。

3. 实现批处理与自动化

   % 示例：批处理预处理 subject_list = {'sub-01', 'sub-02', 'sub-03'}; for i = 1:length(subject_list) sub_id = subject_list{i}; raw_file = fullfile('Data', 'Raw', [sub_id, '_task-motor_raw.nirs']); output_file = fullfile('Data', 'Preprocessed', [sub_id, '_preprocessed.mat']); % 调用你的预处理函数 preprocess_pipeline(raw_file, output_file); fprintf('已完成被试 %s 的预处理\n', sub_id); end

4. 数据与代码版本管理

   • 对分析代码使用 Git 进行版本控制。
   • 原始数据永远只读，任何处理都生成新文件。
   • 在结果文件和图表文件名中包含关键参数和版本信息（如topo_TaskA_HbO_tval_fdr_20240410.png）。

5. 结果验证与敏感性分析

   • 改变关键预处理参数（如滤波带宽），检查统计结果是否稳定。
   • 尝试不同的运动校正算法，比较结果。
   • 使用不同的多重比较校正方法，报告其结果。

6. 生产环境（如实验室共享分析管道）注意事项

   • 编写详细的README.md，说明环境依赖、安装步骤、数据格式要求和运行示例。
   • 考虑将成熟的流程封装成带有图形用户界面的小工具，供不熟悉编程的同事使用。
   • 定期备份完整的分析环境（包括 MATLAB 版本和工具箱版本），以确保长期可重复性。

掌握从数据到图表的完整链条，理解每一步背后的原理和潜在陷阱，你就能从容应对实际研究中的大部分分析任务。接下来，可以进一步探索更高级的主题，如功能连接分析、图论指标计算、机器学习解码等，这些都可以在稳固的基础之上进行拓展。