简介：RQA对离散时间序列进行递归图分析的MATLAB代码包，面向时间序列分析与信号处理方向的学生、研究人员和工程技术人员，提供完整的数据预处理、递归图构建与量化分析流程，可用于心电图、脑电、金融波动和机械故障诊断等实际场景。整个zip压缩包仅1.3MB，包含2个文件：一个MATLAB脚本与一份PDF说明文档，脚本实现了从数据读入到提取复发率、最长对角线长度、平均长度线和熵等核心指标的处理链路，文档则针对每个量化指标给出计算说明与解读建议，并解释距离阈值等关键参数的选择方法。目前已有160人学习下载，适合作为快速入门RQA的轻量示例。通过这套代码，读者可直接将其用于自己的离散序列，修改参数观察指标变化，快速验证非线性动态特征，同时结合PDF说明理解递归图背后的理论，为后续深入分析打下基础。

1. 递归量化分析（RQA）对离散时间序列到底在算什么：一张递归图和七个统计量

很多从业者手里有一段离散时间序列：传感器振动值、设备温度、订单量、心电图R-R间期。用均值、方差和功率谱看，信号平稳得像什么都没发生，可系统实际上已经换了状态。递归量化分析（Recurrence Quantification Analysis）的思路跟线性工具完全不同：先把一维序列嵌入高维相位空间，画一张递归图（Recurrence Plot），再用递归率、确定度、最长对角线等指标把图里的纹理压缩成数字。这份资源是一套MATLAB代码，输入一列时间序列即可完成相位重构、递归图生成和RQA指标计算；适合做故障诊断、状态识别或非线性信号分析的从业者，也适合刚接触非线性时间序列、想拿一套能跑通的代码做对比实验的人。

2. 为什么要看轨迹重访：从递归图到RQA指标的构造逻辑

2.1 相位空间重构和多维向量的距离

一维时间序列看不出系统状态，是因为单个变量只记录了系统在某一方向上的投影。Takens延迟嵌入定理告诉我们：用延迟坐标构造的向量可以恢复系统的相空间拓扑结构。对序列 u(1), u(2), …, u(N)，嵌入维数为 m，延迟时间为 τ，得到的第 i 个状态向量为：

x(i) = [u(i), u(i+τ), u(i+2τ), …, u(i+(m-1)τ)]

注意这里的 i 取值范围是 1 到 N - (m-1)τ。这个有效点数 Nrec 是整个分析里最容易出错的地方，后面避坑章会专门讲。构造好状态向量后，递归图的本质就是一个阈值化的距离矩阵。对任意两个向量 x(i) 和 x(j)，计算两者之间的距离，如果距离小于等于临界半径 ε，就认为系统在 i 时刻和 j 时刻的轨迹发生了“重访”，在递归图上记 1，否则记 0。

距离的取法直接影响递归图纹理。最常用的是欧氏距离，即 L2 范数；也有人用切比雪夫距离（L∞ 范数），它取的是各分量差值的最大值。对同一段数据，L∞ 通常会给出更大的递归区域，因为它的判定条件更宽松：只要有一个分量接近就算递归。我一般会把默认算法写成欧氏距离，切比雪夫作为备用选项，方便在不同任务里快速切换。

2.2 递归矩阵与七个核心指标的物理解释

递归矩阵 R 是 Nrec × Nrec 的 0/1 矩阵，R(i,j) 表示 x(i) 与 x(j) 是否足够接近。对所有递归点做统计分析，就得到RQA指标。最核心的七个指标分别描述递归图的不同侧面。

递归率 RR 是图中递归点占所有点对的比例，反映系统轨迹在相位空间中的“密度”。确定性 DET 是落在长度不小于 Lmin 的对角线段上的递归点数量占递归点总数的比例，反映系统的确定性和周期性。平均对角线长度 L 是所有长度不小于 Lmin 的对角线段的平均长度；最长对角线 Lmax 则是主对角线之外最长的连续对角线段。这两个指标对混沌系统中的轨道分离速度很敏感。对角线熵 ENTR 是对角线长度分布的香农熵，衡量递归图纹理的复杂程度，周期信号熵低，混沌信号熵高。

剩下两个指标面向垂直方向。层流率 LAM 是落在垂直方向连续线段上的递归点占比，捕获系统“停留在某个状态附近”的时间；平均垂直长度 TT 也叫湍流时间，衡量系统在局部状态停留的平均时间。这七个指标组合起来，相当于把一张二维纹理图压缩成一组可以跨信号比较的标量。

2.3 离散序列和连续序列在递归图上的本质差异

连续时间序列的距离分布是连续的，可以选一个合适的 ε 让递归率落在合理范围；但离散序列，尤其是整数值序列，距离只取少数几个值。比如取值为 {0, 1, 2, 3} 的符号序列，任意两个状态向量之间的距离只有有限种可能，同一个 ε 下会出现大片连续递归区域，递归图看起来像一块块色斑，而不是连续信号那种细碎纹理。

这意味着离散序列用固定绝对阈值很容易翻车：ε 设小，递归点太少；ε 加一，递归点猛增。更稳的做法是用固定百分比邻域，也就是对每个点只保留距离最近的固定比例邻域作为递归点，或者用第 4 章介绍的目标递归率校准法。还有一个常见处理是先对离散序列做差分或标准化，再进行相位空间重构，这样能减小取值离散程度对递归图纹理的扭曲。

3. MATLAB实现：三段式代码把时间序列变成可比较的指标

3.1 相位空间重构函数：有效点数的来源

代码包里第一个核心函数是相位空间重构。输入一维列向量 u，返回状态向量矩阵 Y，每一行对应一个相位空间点。

function [Y, Nrec] = phaseSpaceReconstruct(u, m, tau) % 时间延迟嵌入，将一维序列映射到 m 维相位空间 % u : 一维时间序列，列向量 % m : 嵌入维数，正整数 % tau : 时间延迟，正整数，单位是采样间隔倍数 % Y : Nrec x m 的状态向量矩阵 % Nrec: 有效点数，等于 N - (m-1)*tau u = u(:); N = length(u); Nrec = N - (m-1)*tau; if Nrec < 1 error('序列长度不够，请减小 m 或 tau'); end Y = zeros(Nrec, m); for j = 1:m % 第 j 列对应延迟 (j-1)*tau Y(:, j) = u((j-1)*tau + (1:Nrec)); end end

这段代码的关键是起点和终点。每个状态向量最后一个分量是 u(i+(m-1)τ)，所以 i 最大只能取到 N-(m-1)τ。这里输出的 Nrec 就是后续递归矩阵的尺寸。m 和 τ 越大，Nrec 越小，数据浪费越多；对采样点只有几百的短序列，m 不建议超过 5。还要注意 u 被强制转为列向量，避免行向量输入时矩阵形状出错。

3.2 递归矩阵生成：距离计算方式与阈值化

有了状态向量矩阵，下一步是构造递归矩阵。下面这个函数默认用欧氏距离，同时支持切比雪夫距离。

function R = recurrenceMatrix(u, m, tau, eps, normType) % 生成递归矩阵 R，逻辑矩阵，1 表示递归点 % u : 一维时间序列 % m, tau : 嵌入维数与延迟 % eps : 临界半径，标量 % normType: 'euclidean' 或 'chebychev' Y = phaseSpaceReconstruct(u, m, tau); Nrec = size(Y, 1); R = false(Nrec, Nrec); for i = 1:Nrec diffMat = Y - Y(i, :); switch lower(normType) case 'euclidean' d = sqrt(sum(diffMat.^2, 2)); case 'chebychev' d = max(abs(diffMat), [], 2); otherwise error('不支持的范数类型'); end R(i, :) = (d <= eps).'; end end

这段实现用循环逐行计算距离，逻辑清晰但速度一般。对 Nrec 在 1000 以下运行没问题；Nrec 到 5000 时计算时间会明显上升。如果机器上装了 MATLAB 自带的统计与机器学习工具箱，可以把内层循环换成 pdist2，能快一个数量级。R 矩阵用逻辑类型存储，每个元素只占 1 字节，这是控制内存的关键习惯，不要用 double 存 0/1。

3.3 RQA统计：对角线和垂直线段的运行长度编码

计算RQA指标的核心，是对递归矩阵的每个对角线偏移和每个垂直方向做运行长度统计。下面给出一个能直接使用的指标函数：

function rqa = computeRQA(R, Lmin) % R : 逻辑递归矩阵，已去除主对角线 % Lmin : 参与统计的最小线段长度，默认取 2 if nargin < 2, Lmin = 2; end N = size(R, 1); % 递归率：去掉主对角线后的点对总数 % 主对角线全部置为 0，避免自递归点干扰 R = R .* ~eye(N); NR = sum(R(:)); totalPairs = N * (N - 1); rqa.RR = NR / totalPairs; % 对角线方向统计 diagFreq = zeros(1, N); sumDiagLength = 0; maxDiag = 0; for d = 1:N-1 seg = R(1:N-d, 1+d:N); % 这一列依次按行展开后，连续为 1 的段就是一段对角线段 lineVals = seg(:).'; idx = find(diff([0, lineVals, 0]) ~= 0); for k = 1:2:length(idx)-1 len = idx(k+1) - idx(k); if len >= Lmin diagFreq(len) = diagFreq(len) + 1; sumDiagLength = sumDiagLength + len; maxDiag = max(maxDiag, len); end end end rqa.L = sumDiagLength / max(1, sum(diagFreq(Lmin:end))); rqa.Lmax = maxDiag; rqa.DET = sumDiagLength / max(1, NR); % 垂直方向统计 vertFreq = zeros(1, N); sumVertLength = 0; for col = 1:N colVals = R(:, col).'; idx = find(diff([0, colVals, 0]) ~= 0); for k = 1:2:length(idx)-1 len = idx(k+1) - idx(k); if len >= Lmin vertFreq(len) = vertFreq(len) + 1; sumVertLength = sumVertLength + len; end end end rqa.LAM = sumVertLength / max(1, NR); rqa.TT = sumVertLength / max(1, sum(vertFreq(Lmin:end))); % 对角线熵：基于对角线段长度分布 p = diagFreq(Lmin:end) / max(1, sum(diagFreq(Lmin:end))); p = p(p > 0); rqa.ENTR = -sum(p .* log(p)); end

这段代码把每个偏移方向上的序列拉直后做运行长度统计，原理上等价于沿对角线扫描。重点逻辑是 idx 变量定位连续 1 的起止位置，从而实现任意长度的线段统计。DET 的分母用的是 NR，也就是去除主对角线后的总递归点数；如果某个实现把主对角线也算进递归点，DET 会虚高，跨工具对比时要注意。Lmin 默认取 2 是多数文献的惯例，但如果你想复现一篇具体论文，务必先确认论文中的 Lmin 取值。

3.4 一次完整运行：从加载数据到输出指标表和递归图

把三段函数串起来，主脚本只需要几行：

load demoSignal.mat; % 假设数据里有一个变量 sig，是列向量 m = 3; tau = 1; eps = 0.5; Lmin = 2; Y = phaseSpaceReconstruct(sig, m, tau); R = recurrenceMatrix(sig, m, tau, eps, 'euclidean'); R = R .* ~eye(size(R, 1)); % 主动剔除主对角线 rqa = computeRQA(R, Lmin); imagesc(R); title('Recurrence Plot');

运行后输出 rqa.RR、rqa.DET、rqa.L、rqa.Lmax、rqa.ENTR、rqa.LAM、rqa.TT。我习惯把 R 矩阵直接显示成黑白图，先肉眼看纹理再做数值统计，因为很多异常信号在递归图上会有明显特征：突变表现为横向或纵向的断续线，漂移表现为递归图边缘变暗。建议每一步都打印 Nrec，确认有效点数和预期一致。

4. 三个关键参数怎么设：嵌入维数、时间延迟与临界半径的实操

4.1 嵌入维数m：用Cao方法代替人工试错

嵌入维数 m 决定状态向量的长度。m 太小，相空间没有充分展开，距离计算掩盖系统结构；m 太大，每个状态向量需要的延迟点数变多，Nrec 变小，统计噪声变大。实际工作里很少直接试 m，我一般用 Cao 方法，它不需要先验地估计噪声水平，适合短序列。

Cao 方法的核心是计算 E1(d) = E(d+1) / E(d)，其中 E(d) 是基于最近邻距离定义的均值。当 E1(d) 随着 d 增大开始逼近 1 并保持平稳时，对应的 d 就是合适的嵌入维数。下面的函数计算 E1 随维数的变化：

function E1 = caoE1(u, maxDim, tau) % 简化版 Cao 方法，只输出 E1 指标 u = u(:); for d = 1:maxDim Y = phaseSpaceReconstruct(u, d, tau); Nrec = size(Y, 1); % 对每个点找最近邻距离 nnDist = zeros(Nrec, 1); for i = 1:Nrec dst = sqrt(sum((Y - Y(i,:)).^2, 2)); dst(i) = inf; % 排除自身 nnDist(i) = min(dst); end a = zeros(Nrec, 1); % 计算加入下一维后的距离变化 Y1 = phaseSpaceReconstruct(u, d+1, tau); if size(Y1, 1) < Nrec Nrec = size(Y1, 1); Y = Y(1:Nrec, :); nnDist = nnDist(1:Nrec); end for i = 1:Nrec dstPlus = sqrt(sum((Y1 - Y1(i,:)).^2, 2)); dstPlus(i) = inf; nnPlus = min(dstPlus); a(i) = dstPlus / max(nnDist(i), eps); end E(d) = mean(a); end E1 = E(2:end) ./ E(1:end-1); plot(2:maxDim, E1, '-o');

这段代码没有做文献里的分段平均改进，但对大多数中等长度序列足够用。实际操作时，我会把 E1 从 2 画到 8，观察曲线第一次进入“平台区”的转折点。例如 E1 在 d=4 之后稳定在 1 附近且没有继续明显上升，就选 m=4。Cao 方法输出的不是一个唯一确定值，而是一个合理区间，所以不必纠结分数差异。

4.2 时间延迟τ：互信息极小值与自相关法

τ 决定每个状态向量两个相邻分量的间隔。τ 太大，相邻分量几乎不相关，重构轨迹被拉散；τ 太小，分量之间冗余，递归图纹理被压缩。最常用的两套选择法：自相关法和互信息法。

自相关法简单直接，计算序列自相关函数，找到它首次下降到初始值的 1-1/e 约 36.8% 时的滞后，作为 τ。互信息法更全面，它计算原始序列和延迟序列之间的互信息 I(τ)，取第一个极小值对应的 τ。互信息能捕捉非线性相关性，所以对混沌信号更可靠，但计算量稍大。一个简化版的互信息计算如下：

function I = mutualInfo(u, tau) % 简化版互信息：用直方图估计联合分布 u = u(:); N = length(u); x = u(1:N-tau); y = u(1+tau:N); nbins = max(5, round(sqrt(numel(x)))); [c, edges] = histcounts(x, nbins); px = c / sum(c); [~, edgesY] = histcounts(y, nbins); py = histcounts(y, edgesY) / numel(y); pJoint = histcounts2(x, y, edges, edgesY) / numel(x); logTerm = pJoint ./ (px(:) * py(:).'); logTerm(pJoint == 0) = 1; % 避免 log(0) I = nansum(pJoint .* log(abs(logTerm))); end

这个实现没有做互信息的均值校正，但找极小值的趋势是可靠的。我的习惯是把 τ 从 1 到 10 都算一遍，画出 I(τ) 曲线，取第一个明确极小值。如果曲线太平、没有清晰极小值，就用 τ=1 作为默认值，因为对采样率较高的离散信号，过大的 τ 会丢失时间结构。

4.3 临界半径ε：固定阈值、百分比阈值与目标递归率

ε 是 RQA 里最敏感的参数，它直接决定递归点是稀疏还是密集，也影响 DET 和 Lmax。经验上，ε 太小，递归点太少，对角线统计噪声大；ε 太大，递归点太多，几乎全部点都变成递归点，DET 趋近 1 失去区分度。说直白点，ε 的选取在有些论文里写得像玄学，实质是尺度选择问题。

三种常用方案各有边界。固定绝对阈值适合所有被比较的序列具有相同量纲和幅值的场景，比如同一传感器的多段振动数据；百分比阈值是把每个点到其所有近邻的距离排序，取某个分位数作为该点的自适应阈值，适合幅值漂移信号；目标递归率校准是先设定期望的 RR，比如 0.05，再用二分法搜索 ε，使实际 RR 逼近目标值。第三种方案对跨工况对比最友好，我会在最后一章给出具体代码。

4.4 常用参数组合表与默认建议

下面这个表是跨多个使用场景的经验区间，不是硬性标准，用于给第一次跑这套代码的人一个参考起点：

拿到新数据时，我会先用 m=4、tau=1、目标RR=0.05 作为默认粗跑一遍，再根据递归图纹理和具体任务调整。不要盲目追求某个“最优”参数，RQA 指标在 m 和 τ 的合理邻域内变化通常不大，真正会导致指标剧烈变化的是 ε，所以参数优先级是 ε > τ > m。

5. 避坑指南：在离散序列上跑RQA的六个翻车现场

5.1 主对角线剔除与否影响RR和DET

现象：同一段数据用两套代码算 RR 和 DET，结果差别很大，一套 RR=0.04，另一套 RR=0.06。

原因：主对角线上 x(i) 和它自身的距离恒为 0，只要 ε>0 就是递归点。如果实现把主对角线计入递归点总数，RR 和 DET 的分母分子同时变化，指标虚高；不同实现是否剔除主对角线不一致，导致跨工具对比直接失真。

解决：在统计前强制把主对角线置 0，也就是代码中的R = R .* ~eye(N)，并且把分母固定为 N(N-1)。这个决定需要在实验记录里写明，避免自己和读者误判。

5.2 序列长度不同，指标能不能比

现象：A 序列 1000 点、B 序列 500 点，B 的 DET 明显更高，看起来 B 更“确定”。

原因：短序列的有效点数少，边界效应影响大，递归图边缘部分不完整；同时 Lmax 天然受序列长度限制，短序列的 Lmax 相对噪声更大。

解决：比较前把各序列截取成相同长度，或采用滑动窗口对全序列提取指标曲线。如果必须在不同长度下比较，优先使用 RR、DET 和 LAM 这些比值型指标，Lmax 只能作为辅助，且最好用 Lmax/Nrec 归一化。

5.3 大N递归矩阵的内存瓶颈

现象：Nrec 超过 5000 时计算变慢，Nrec 超过 8000 时 MATLAB 内存占用巨大，甚至直接无响应。

原因：递归矩阵是 Nrec × Nrec 的逻辑矩阵，Nrec=8000 时也有 64M 个元素，如果中途用 double 保存距离矩阵，内存占用立刻翻 8 倍。

解决：只保存 R 的逻辑矩阵，不保留完整距离矩阵；距离计算用分块或逐行处理。如果只看指标不需要绘图，可以逐条偏移计算对角线统计，跳过整矩阵存储。代码包里专门提供了一个内存友好版本，用分块循环替代全矩阵构造。

5.4 离散整数序列的“块状递归图”

现象：取值为 {0, 1, 2} 的序列画递归图，大范围区域整块变黑，几乎看不出纹理。

原因：离散取值少，状态向量重复率高；用绝对距离阈值时，很多不同状态的距离也小于 ε，递归点大面积连成片。

解决：先做序列标准化或一阶差分得到连续分布的数据，再用百分比阈值。对符号序列，要明确说明结果属于“符号化递归”，与连续信号的 RQA 指标不在同一尺度上，不能直接混比。

5.5 DET高不等于混沌：周期与混沌的区分

现象：纯正弦信号的 DET 接近 1，Logistic 混沌映射的 DET 也能到 0.9 以上，单看 DET 完全分不开。

原因：DET 只反映轨迹返回的确定性结构，周期性数据和混沌数据都能产生明显的对角结构。它们真正的差异在 Lmax 和 ENT R：周期信号最长对角线接近 Nrec，熵低；混沌信号最长对角线较短，对角线长度分布更随机，熵较高。

解决：比较时同时看 Lmax/Nrec 和 ENTR。周期信号两者呈“长而低熵”组合，混沌信号呈“短而高熵”组合。如果只报 DET，结论说服力不够。

5.6 工具箱依赖导致代码跑不起来

现象：代码下载后直接运行，报错提示缺少某个函数或工具箱。

原因：部分公开实现依赖统计工具箱里的 pdist2、图像处理工具箱里的连通域函数，而基础版 MATLAB 不带这些模块。

解决：资源内的实现全部基于基础 MATLAB，不依赖任何附加工具箱。如果要用 pdist2 加速，运行前先判断license('test', 'Statistics_Toolbox')，不可用时回退到循环版本。这样既保证可复现，又能照顾有工具箱的用户。

6. 进阶技巧：滑动窗口RQA与目标递归率校准

6.1 滑动窗口RQA的快速实现

RQA 指标的静态值只适合平稳段分析，对非平稳信号和状态切换检测，需要把序列切成窗口滑动计算。窗口长度建议取 200 到 1000 点，步长取窗口的 1/4 到 1/2。核心思路是对每个窗口重新做相位空间重构和 ε 校准，得到一组随时间变化的指标序列。

function rqaSeries = slidingRQA(u, W, S, m, tau, targetRR) tAxis = W:S:length(u); rqaSeries = zeros(numel(tAxis), 7); k = 0; for t = tAxis k = k + 1; seg = u(t-W+1:t); epsOpt = calibrateEpsilon(seg, m, tau, targetRR); R = recurrenceMatrix(seg, m, tau, epsOpt, 'euclidean'); R = R .* ~eye(size(R, 1)); rqa = computeRQA(R, 2); rqaSeries(k, :) = [rqa.RR, rqa.DET, rqa.L, rqa.Lmax, rqa.ENTR, rqa.LAM, rqa.TT]; end end

6.2 校准ε让不同信号在同样递归密度下比较

对每个窗口重新校准 ε，是为了保证滑动过程中 RR 始终接近目标值，这样才能把 DET、Lmax 的变化归因于动力学结构变化，而不是递归密度变化。二分法搜索在 ε 与 RR 单调递增的前提下可靠，迭代 20 次足够收敛。

function epsOpt = calibrateEpsilon(u, m, tau, targetRR) Y = phaseSpaceReconstruct(u, m, tau); Nrec = size(Y, 1); % 计算距离上界，用于设置二分区间 maxDist = 0; for i = 1:Nrec dst = sqrt(sum((Y - Y(i,:)).^2, 2)); maxDist = max(maxDist, max(dst)); end lo = 0; hi = maxDist; for iter = 1:20 mid = (lo + hi) / 2; R = recurrenceMatrix(u, m, tau, mid, 'euclidean'); RR = sum(R(:) & ~eye(Nrec)) / (Nrec * (Nrec - 1)); if RR < targetRR lo = mid; else hi = mid; end end epsOpt = (lo + hi) / 2; end

这段代码会随递归矩阵计算重复运行，代价不小，但胜在稳定。我之前用固定 ε 对比两段不同工况的振动数据，结果一个 RR 是 0.002，另一个是 0.15，DET 完全不可比；后来写成目标递归率校准，才发现两段信号在动力学复杂度上其实很接近。从那以后，凡是需要跨工况或跨样本比较 RQA 指标，我都强制走一遍“定窗口长度—校准 RR—检查递归图边缘是否有截断”的流程，希望帮到你。

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