简介：msdanalyzer 是一套面向生物物理、化学、材料与工程领域科研人员的 MATLAB 均方位移分析类，用于研究分子动力学、细胞迁移、颗粒扩散等场景中粒子的运动模式。它把轨迹导入、MSD 计算、扩散系数拟合与可视化封装为可复用的类方法，适合具备一定 MATLAB 基础、希望快速搭建分析流程的研究者。资源包共 24 个文件，以 22 个 .m 源码文件为主，另含 1 个 README.md 说明文档与 1 个 .gitignore，整体约 24KB，体积轻巧便于直接嵌入现有工程。源码覆盖 MSD 计算、均值与对数拟合、速度自相关、漂移校正、轨迹与曲线绘制及误差阴影等模块，并附带示例数据与使用指南，可帮助读者理解延迟时间、均方根位移、扩散系数等关键概念，按需调整时间窗口与步长参数，并保存结果用于共享。目前已有 921 人学习下载。
1. 单粒子追踪数据到手之后，msdanalyzer 到底替你省掉了什么
做单粒子追踪的人都有个共同体验：显微镜跑完一夜，拿到一堆轨迹坐标，真正的折磨才刚开始。均方位移（MSD）曲线怎么算、滞后时间怎么取、几条轨迹能不能平均、拟合出来的扩散系数到底可不可信，每一步都藏着主观选择。msdanalyzer 就是冲这件事来的——它是一个 MATLAB 类，把 MSD 计算、轨迹筛选、曲线拟合、误差估计这一整套流程封装成对象方法，让你用几行代码就能从原始轨迹走到可发表的扩散系数。它适合两类人：刚接手追踪数据、不想自己从零写 MSD 循环的新手，以及已经有一堆脚本、但想找个统一框架管理多组对比实验的熟手。这一章先把它的定位和边界讲清楚，后面几章再拆开揉碎讲怎么用、参数怎么调、哪里会翻车。
2. 把 msdanalyzer 对象建起来：从轨迹矩阵到第一条 MSD 曲线
2.1 为什么用类而不是一堆散装函数
散装脚本最大的问题是状态满天飞。你算完 MSD，想换个滞后时间范围重新拟合，得回头改参数、重跑、再手动记录结果。轨迹一多，哪个结果对应哪组参数就全靠文件名硬记，时间一长自己都不信。msdanalyzer 用面向对象的方式把这些状态收进一个对象里：轨迹、MSD 曲线、拟合结果、拟合参数都挂在同一个实例上，你调用方法改的是这个对象，查的也是这个对象。好处是实验可复现——同一个对象序列化存下来，别人加载后能看到你当时用的每一条轨迹和每一个拟合区间。
常见做法是每个实验条件建一个 msdanalyzer 实例，比如对照组一个、处理组一个，最后把两者的 MSD 曲线画在同一张图上对比。这样不会出现「我到底把哪组数据平均进去了」这种血泪问题。
2.2 轨迹数据的组织格式与导入
msdanalyzer 对输入轨迹的格式有明确要求：它期望的是轨迹的单元数组，每个单元是一条轨迹，每条轨迹是一个 N×3 的矩阵，三列分别是时间、x 坐标、y 坐标（如果是三维追踪就是 N×4，多一列 z）。时间列必须是单调递增的，单位统一，通常用秒。
下面这段代码演示怎么把一份常见的 CSV 轨迹文件读进来、整理成 msdanalyzer 能吃的格式。假设你的 CSV 每行是track_id, frame, time, x, y。
% 读取原始轨迹 CSV，整理成 msdanalyzer 需要的单元数组格式 raw = readmatrix('trajectories.csv'); % 列: id, frame, time, x, y ids = unique(raw(:,1)); tracks = cell(numel(ids), 1); for k = 1:numel(ids) idx = raw(:,1) == ids(k); seg = raw(idx, :); seg = sortrows(seg, 3); % 按时间列排序，保证单调递增 tracks{k} = [seg(:,3), seg(:,4), seg(:,5)]; % [time, x, y] end % 过滤掉点数太少的轨迹，少于 10 个点算 MSD 没意义 minPoints = 10; tracks = tracks(cellfun(@(t) size(t,1) >= minPoints, tracks)); fprintf('有效轨迹数: %d\n', numel(tracks));逻辑说明：readmatrix把 CSV 读成数值矩阵，unique拿到所有轨迹编号。循环里按编号切片、按时间排序，这一步不能省——追踪软件导出的数据偶尔会有帧顺序错乱，时间列不单调会让后面 MSD 计算出负的滞后时间。最后用cellfun过滤短轨迹，因为 MSD 在滞后时间接近轨迹长度时统计量极少，噪声极大，留着只会污染平均曲线。
参数说明：minPoints是最小点数阈值，我一般设 10 到 20。设太小，短轨迹的 MSD 尾部会剧烈上翘；设太大，长轨迹本来就不多，样本量又不够。这个值要根据你的平均轨迹长度来定，经验是取平均长度的十分之一左右。
2.3 创建对象并计算 MSD
数据整理好之后，创建 msdanalyzer 对象并计算 MSD 就是两三行的事。
% 创建 msdanalyzer 对象，指定时间与空间单位 ma = msdanalyzer(2, 's', 'um'); % 2 表示二维追踪 % 把轨迹塞进对象 ma = ma.addAll(tracks); % 计算 MSD，指定最大滞后时间为轨迹平均长度的四分之一 ma = ma.computeMSD; % 查看对象里现在有什么 ma逻辑说明：构造函数第一个参数是维度，二维追踪传 2，三维传 3。后面两个参数是时间单位和长度单位，纯粹为了画图时坐标轴标签好看，不影响计算。addAll把单元数组一次性加进去，也可以逐条addTrack。computeMSD是核心方法，它会对每条轨迹单独算 MSD，再按滞后时间点做平均，得到整体 MSD 曲线。
参数说明：computeMSD默认会用所有可用的滞后时间，但你可以传一个最大滞后比例参数来控制。常见做法是限制在轨迹长度的四分之一以内，因为超过这个范围，能参与平均的轨迹对数量骤减，曲线尾部会变成玄学。这个比例没有绝对标准，1/4 到 1/3 都有人用，关键是在论文里写清楚你用了多少。
2.4 画出第一条 MSD 曲线并读懂它
% 画出平均 MSD 曲线，带误差棒 figure; ma.plotMSD; xlabel('滞后时间 \tau (s)'); ylabel('MSD (\mum^2)'); title('平均 MSD 曲线');画出来的曲线通常分两段：前一小段接近线性，斜率对应短时扩散系数；后面逐渐变平或上翘，反映受限扩散或主动运输。如果曲线从头到尾都是直线，说明是纯布朗扩散；如果先线性后平台，说明粒子被限制在某个区域里；如果上翘，可能有定向流。读懂这三类形状，是后面选拟合区间的前提。
3. 拟合与参数提取：扩散系数不是随便取一段斜率就完事
3.1 MSD 曲线的三种典型形状与拟合区间选择
拟合 MSD 最容易被忽视的一步是选区间。很多人直接把整条曲线丢进线性拟合，结果扩散系数偏小，因为尾部平台或上翘把斜率拉平了。正确做法是先看曲线形状，再决定拟合哪一段。
纯布朗扩散的 MSD 曲线在全部滞后时间上都是线性的，理论上拟合整条没问题，但实际数据尾部噪声大，我一般只拟合前 1/4 到 1/3。受限扩散的曲线前段线性、后段平台，拟合必须只取线性段，平台段要单独分析受限尺寸。主动运输的曲线整体上翘，线性拟合没有物理意义，得改用二次项拟合或者分段处理。
判断线性段结束的位置，一个实用技巧是看 log-log 图的斜率：纯扩散斜率是 1，受限扩散斜率会降到 0，主动运输斜率大于 1。在斜率开始偏离 1 的地方切断，就是线性段的上界。
3.2 用 fitMSD 做线性拟合与参数提取
msdanalyzer 提供了fitMSD方法，可以对指定滞后时间区间做线性拟合，提取扩散系数。
% 对滞后时间 0.1 到 1.0 秒的区间做线性拟合 fitRange = [0.1, 1.0]; ma = ma.fitMSD(fitRange); % 查看拟合结果 ma.fitSummary逻辑说明：fitMSD接收一个两元素向量，表示拟合的滞后时间上下界。它内部会对每个滞后时间点的 MSD 值做加权线性拟合，权重通常取该点参与平均的轨迹对数量的倒数，这样尾部噪声大的点权重低，不会带偏斜率。拟合完成后，fitSummary会打印出扩散系数、截距和拟合优度。
参数说明：fitRange的选取直接决定结果。下界不能取 0，因为滞后时间趋近 0 时 MSD 也趋近 0，这个点对斜率没有信息量，反而可能因为定位误差导致截距偏移。上界按 3.1 节的方法确定。如果拟合优度低于 0.95，说明区间选得不对，要么混进了非线性段，要么区间太短噪声主导。
3.3 多组实验的批量拟合与结果对比
实际实验往往有多个条件，每个条件一个 msdanalyzer 对象。批量拟合可以写成循环，把结果汇总成表格。
% 假设 conditions 是多个 msdanalyzer 对象的单元数组 conditions = {maControl, maTreat1, maTreat2}; names = {'对照', '处理1', '处理2'}; fitRange = [0.1, 1.0]; results = table(); for k = 1:numel(conditions) obj = conditions{k}; obj = obj.fitMSD(fitRange); D = obj.getDiffusionCoefficient; % 提取扩散系数 results.Name{k} = names{k}; results.D(k) = D; results.N(k) = numel(obj.tracks); end disp(results);逻辑说明：循环里对每个对象调用fitMSD，然后用getDiffusionCoefficient把拟合得到的扩散系数取出来。汇总成 table 方便后续导出或画柱状图。注意每个对象的拟合区间可以不同——如果处理组的轨迹长度和对照组差异大，强行用同一个区间反而不合理。
参数说明：getDiffusionCoefficient返回的是拟合斜率的四分之一（二维情况），因为 MSD 在二维下的理论表达式是 4Dτ。如果你做的是三维追踪，系数是六分之一。这个换算 msdanalyzer 内部会处理，但你要清楚自己拿到的数是什么含义，别在论文里写错公式。
4. 避坑与排查：那些让 MSD 结果不可信的细节
4.1 轨迹拼接错误导致 MSD 整体偏移
现象：MSD 曲线在滞后时间趋近 0 时截距明显不为 0，甚至出现负值。
原因：追踪软件在粒子暂时消失又出现时，可能把两条不同粒子的轨迹拼成了一条，或者把同一条轨迹断成了两段。拼接错误会让坐标序列出现跳变，MSD 在短滞后时间就被抬高。
解决：导入轨迹后先画几条原始轨迹的 x-y 图，肉眼检查有没有突然的跳跃。也可以用速度阈值自动筛：相邻帧位移超过某个物理上限的轨迹段标记出来。msdanalyzer 本身不做轨迹拼接校正，这一步得在导入前完成。
4.2 滞后时间点数不足导致拟合不稳定
现象：拟合优度忽高忽低，换个区间扩散系数变化很大。
原因：拟合区间内只有三四个滞后时间点，线性拟合的自由度太低，一个点的波动就能改变斜率。
解决：拟合区间至少包含 5 到 8 个滞后时间点。如果数据本身帧率低、轨迹短，凑不够点数，那就得考虑提高采集帧率或延长追踪时间，而不是硬拟合。
4.3 定位误差被误当成扩散
现象：所有条件的扩散系数都偏大，且不同条件之间差异不明显。
原因：定位误差会给 MSD 曲线在短滞后时间处贡献一个常数项，表现为截距不为零。如果拟合区间下界取得太低，这个常数项会被部分吸收进斜率，导致扩散系数虚高。
解决：把拟合下界提高到定位误差影响可忽略的滞后时间。一个粗略判断方法是看 MSD 曲线在最短几个滞后时间点是否偏离线性——如果前几个点明显高于线性外推，说明定位误差在起作用。也可以单独用固定粒子测定位误差，算出对应的 MSD 偏移量，从数据中扣除。
4.4 单位不统一导致数量级错误
现象：扩散系数比文献值差几个数量级。
原因：时间列用了毫秒但当成秒，或者坐标用了像素但当成微米。
解决：导入前统一单位。msdanalyzer 构造函数里的单位参数只是标签，不会做换算。像素到微米的换算系数要用标定数据，别用估计值。我一般会在导入脚本里加一行断言，检查时间列的最大值是否在合理范围内，比如追踪 100 秒的数据，时间列最大值应该在 100 左右，如果看到 100000 就是毫秒没换算。
4.5 轨迹数量太少导致平均曲线噪声大
现象：MSD 曲线毛刺严重，误差棒比均值还大。
原因：每个滞后时间点参与平均的轨迹对数量不够。MSD 平均的统计误差大致按参与对数的平方根倒数下降。
解决：要么增加追踪的粒子数，要么放宽轨迹筛选阈值保留更多短轨迹。但放宽阈值会引入短轨迹尾部上翘的问题，需要在样本量和曲线质量之间权衡。一个折中做法是分层平均：长轨迹和短轨迹分别算 MSD，只在短滞后时间区间合并。
5. 进阶技巧：用 msdanalyzer 做变滞后时间分析与结果验证
前面讲的都是固定区间拟合，但有些体系用固定区间并不合适，比如扩散系数随时间变化的非稳态体系，或者 MSD 曲线没有明显线性段的复杂体系。这时候可以用变滞后时间分析：不拟合整段，而是逐点计算局部斜率，得到扩散系数随滞后时间的变化曲线。
% 变滞后时间分析：逐点计算局部斜率 msdData = ma.msd; % 取出平均 MSD 数据 tau = msdData(:,1); msdVal = msdData(:,2); % 用相邻点差分算局部斜率，再换算成扩散系数 localD = diff(msdVal) ./ diff(tau) / 4; % 二维除以 4 tauMid = (tau(1:end-1) + tau(2:end)) / 2; figure; loglog(tauMid, localD, 'o-'); xlabel('滞后时间 \tau (s)'); ylabel('局部扩散系数 D (\mum^2/s)'); title('扩散系数随滞后时间的变化');逻辑说明：diff做相邻点差分，得到每个滞后时间区间的局部斜率，除以 4 换算成扩散系数。横坐标取区间中点。如果局部扩散系数在短滞后时间高、长滞后时间低，说明存在受限扩散或测量噪声；如果基本恒定，说明纯布朗扩散，固定区间拟合的结果可信。
参数说明：差分对噪声敏感，如果曲线毛刺多，可以先对 MSD 数据做轻度平滑再差分。平滑窗口不要太大，否则会抹掉真实的短时行为。这个分析适合用来验证固定区间拟合的结果是否合理，而不是替代它。
另一个实用技巧是结果验证：把拟合得到的扩散系数代回 MSD 公式，画出理论曲线叠在实验曲线上，看两者在拟合区间内是否重合。如果理论线在区间两端明显偏离，说明区间选得不对。这个验证步骤花不了几分钟，但能避免很多审稿意见。
我自己养成的习惯是：每次拟合完，先不看扩散系数的数值，而是先画叠加图。图对了，数才可信。图不对，数再漂亮也是自欺欺人。这个习惯帮我拦下过好几次区间选错的问题，希望帮到你。
本文还有配套的精品资源，点击获取
与LALR：语法分析器的精度-代价权衡本质)
