简介面向时间序列分析与信号处理研究人员这份MATLAB代码包实现了循环量化分析RQA可直接对离散时间序列构建递归图并量化系统动态特征。压缩包共2个文件包括一个.m脚本核心算法与一份pdf说明文档总计约1.3MB轻量易用。已有160人学习下载。通过调整距离阈值等参数可计算复发率、最大/平均对角线长度及熵等指标适用于生物医学、金融、机械等领域的非线性动态分析。借由配套文档使用者可快速掌握递归图构建与结果解读并迁移至自身数据场景中是开展复杂性研究的实用工具。1. RQA不是画图是给“看不见的周期”算一笔账很多人第一次看到递归图都会觉得它只是把时间序列变成了一幅黑白花纹图。真正有价值的是隐藏在那幅图背后的 RQA递归量化分析Recurrence Quantification Analysis它把一维离散时间序列重建到相空间再把相空间中状态的往返程度压缩成一组可比较的数值。这套方法对短序列、带噪声序列、非平稳序列都相当能打尤其适合那些“肉眼看不到周期性、FFT也糊成一片”的场景。我见过不少做振动监测、生理信号分析、气象序列建模的人拿着一段看起来完全随机的数据用 RQA 一算DET 和 Lmax 两个指标立刻把潜在动力学结构挑了出来。这篇文章会把 RQA 对离散时间序列做递归图分析的完整链路讲清楚相空间重构怎么做、递归图和 RQA 指标怎么读、MATLAB 代码怎么组织、参数怎么定、哪些地方最容易翻车。适合想在 MATLAB 里把递归图分析真正跑进自己工程流程的工程师和研究者。2. 递归图与RQA指标先搞清楚图里的线在说什么2.1 相空间重构为什么一维序列要先升维再分析递归图的分析对象不是原始序列本身而是它的相空间重构结果。系统状态往往需要多个变量才能完整描述但我们手上通常只有一维观测值。Takens 嵌入定理保证了在合适的嵌入维数 m 和时间延迟 τ 下用延迟坐标构造的向量可以在拓扑意义上还原原系统的动力学。延迟坐标嵌入的构造方式很直接对离散时间序列 x(1), x(2), …, x(N)取 m 个延迟分量为一个状态向量X(i) [x(i), x(iτ), x(i2τ), …, x(i(m-1)τ)]于是原始序列被转换成 N - (m-1)τ 个状态向量这些向量构成相空间中的一条轨迹。对离散时间序列来说这一步有天然优势采样点本身就是等间隔的不需要像连续方程那样先做数值积分或插值拿来就能用。2.2 从距离矩阵到递归图阈值一加矩阵变图案相空间轨迹构造好后计算任意两个状态向量之间的距离得到一个距离矩阵。然后选定一个阈值 ε距离小于等于 ε 的两个状态被认为是“递归”的在图上记为 1否则记为 0。递归矩阵的数学表达是 R(i,j) Θ(ε - ‖X(i) - X(j)‖)Θ 是阶跃函数‖·‖ 常用欧氏距离。这个 0-1 矩阵可视化之后就是递归图。读图有非常直观的规律规则排列的平行长对角线说明系统具有周期性或准周期性短对角线碎块为主说明系统确定但混沌轨道在不断发散和折叠均匀散落无结构的孤立点说明序列接近随机过程出现明显的水平或垂直条带往往对应层流状态或锁频现象图上一角发黑、远离主对角线处递归点大量消失则提示数据可能带着趋势或发生了非平稳漂移。在 MATLAB 里画递归图我一般用 imagesc 加 gray 色彩映射几行就能出图figure; imagesc(rmat); colormap(gray); axis square; xlabel(i); ylabel(j); title(Recurrence Plot of Discrete Time Series);这里注意 imagesc 默认会翻转 y 轴方向看起来像是图被上下颠倒。想保持矩阵坐标和视觉方向一致可以在后面加一行set(gca, YDir, normal)但这只是一个阅读习惯问题不影响 RQA 数值计算。真正影响数值的是下一节里的统计口径。2.3 RQA核心指标RR、DET、Lmax、ENTR、LAM怎么读递归图本身是定性的RQA 的价值在于把它变成一组可量化、可对比的指标。以下是六个最常用指标及其工程含义。指标计算含义物理含义典型用途RR 递归率递归点占总矩阵元素比例系统状态在相空间中的“回头率”反映状态密度判断信号是否接近随机、状态是否发生整体变化DET 确定率长度不小于 lmin 的对角线段上的递归点占比系统确定性和可预测性的估计区分周期/混沌/噪声的关键指标Lmax 最长对角线长度递归图中最长对角线段长度轨道发散程度的倒数与最大李雅普诺夫指数相关估计混沌强度Lmax 越短说明发散越快L 平均对角线长度所有有效对角线的加权平均长度轨道保持相近状态的持续时间配合 DET 判断系统“稳定程度”ENTR 对角线熵对角线长度分布的香农熵系统结构复杂度分布越多样化熵越高区分规则运动和复杂动力学LAM 层流率垂直结构上的递归点占比系统停留在某个状态附近的时间比例检测间歇性、层流状态、状态切换TT 平均层流时间垂直结构长度的加权平均值层流状态的平均持续长度与 LAM 配合描述间歇过程的持久性这里有两个口径问题后面避坑章节会详细展开lmin 默认取 2主对角线算不算进统计会直接影响 DET 和 L 的绝对值。RQA 没有统一官方的“标准实现”不同论文、不同工具箱的默认值可能不同所以做对比之前必须先确认这些细节。3. 用MATLAB从零搭一套RQA管线代码骨架与最小可运行demo3.1 项目结构一个demo脚本加三个函数模块拿到一个 RQA 相关的 MATLAB 代码包第一步不是急着读源码而是先理解它的模块划分。常见做法是把整套流程拆成四个文件一个主脚本负责加载数据和展示结果两到三个函数模块分别负责相空间重构、递归矩阵生成和指标计算。这样拆的好处是后续换数据只需要改主脚本调参数时也不用在几百行代码里翻找。我建议的最小结构是main_demo.m 作为入口rp_reconstruct.m 做相空间重构rp_matrix.m 根据相空间矩阵和阈值生成递归图rqa_metrics.m 统计对角线分布并计算 RQA 指标。这样的划分也方便你在自己的项目里把指标计算替换成其他 RQA 变体。3.2 相空间重构与递归图生成30行以内能跑通先写相空间重构函数。输入是一维离散时间序列输出是状态向量矩阵 Y。需要注意边界长度N 个采样点重构后只剩下 N - (m-1)τ 个状态向量尾部数据会被截掉。function Y rp_reconstruct(x, m, tau) % 相空间重构: 延迟坐标嵌入 % 输入: % x - 一维离散时间序列, N×1 列向量 % m - 嵌入维数, 正整数 % tau - 时间延迟, 以采样点为单位, 正整数 % 输出: % Y - Np×m 相空间状态矩阵, Np N - (m-1)*tau N length(x); Np N - (m-1)*tau; if Np 0 error(序列长度太短, 请减小 m 或 tau); end Y zeros(Np, m); for j 1:m Y(:, j) x((j-1)*tau (1:Np).); end end这个函数的核心逻辑是用循环把延迟分量逐个填入矩阵。(j-1)*tau (1:Np)是当前延迟分量在原序列中的索引j 从 1 到 m 对应从当前时刻到 m-1 个延迟之后。用循环而不是直接调 hankel 或 toeplitz是为了让边界逻辑对新手更透明也方便后续修改延迟步长。接下来是递归矩阵生成。这里最关键的是内存控制如果先算一个 Np 乘 Np 的完整 double 距离矩阵序列稍长就会内存溢出。我一般按块计算每次只处理一个行块距离矩阵用完即弃。function rmat rp_matrix(Y, eps_th) % 根据相空间矩阵和阈值生成二值递归矩阵 % 输入: % Y - Np×m 相空间矩阵, 来自 rp_reconstruct % eps_th- 递归阈值, 状态向量欧氏距离 eps_th 记为递归点 % 输出: % rmat - Np×Np 的逻辑矩阵, 1 表示递归 Np size(Y, 1); blk 500; % 分块大小, 控制内存峰值 rmat false(Np, Np); % 逻辑矩阵每个元素只占1字节 for i 1:blk:Np idx i:min(iblk-1, Np); D pdist2(Y(idx, :), Y, euclidean); rmat(idx, :) (D eps_th); end % 主对角线状态与自身距离为0, 必然满足阈值, 自动为1 end分块参数 blk 决定内存峰值的上限。block 行乘 Np 列的 double 矩阵Np5000 时500×5000×8 字节约 20MB完全可以接受。如果 Np 再大把 blk 调小到 200 即可。pdist2 来自统计与机器学习工具箱如果你的环境没有这个工具箱可以用内层循环配合sqrt(sum((Y - Y(i,:)).^2, 2))替代但速度会慢一些生产环境建议优先装工具箱。3.3 RQA指标计算对角线分布统计是核心RQA 指标的计算核心是对角线段和垂直线段的长度分布统计。递归图上一条连续由 1 组成的对角线代表两个状态轨迹在相空间中并行运动了 l 个时间单位。统计所有偏移量上的连续 1 段长度得到分布 P(l)。function P line_hist(rmat, lmin) % 统计递归矩阵对角线长度分布 % 输入: % rmat - Np×Np 逻辑递归矩阵 % lmin - 有效线段最小长度, 默认2 % 输出: % P - Np×1向量, P(l)为长度为l的对角线段数量 if nargin 2, lmin 2; end N size(rmat, 1); P zeros(N, 1); for k -(N-1):(N-1) d diag(rmat, k); % 偏移k的对角线 d [0; d(:); 0]; % 前后补0, 便于检测连续段边界 starts find(d(2:end) 1 d(1:end-1) 0); ends find(d(1:end-1) 1 d(2:end) 0); for j 1:numel(starts) L ends(j) - starts(j); if L lmin P(L) P(L) 1; end end end end这个函数逐条提取不同偏移量的对角线通过检测“0 变 1”和“1 变 0”来切分连续线段。示例里对每条对角线调用 diag可读性好但性能一般。Np 在 1500 以下时够用超过这个规模就要考虑换更高效的扫描策略避坑章节会专门展开。有了长度分布RQA 指标就是简单的加权求和。下面是完整的指标计算函数function out rqa_metrics(rmat, lmin, vmin) % 计算RQA六项核心指标 % 输入: % rmat - Np×Np 逻辑递归矩阵 % lmin - 对角线最小统计长度, 默认2 % vmin - 垂直线最小统计长度, 默认2 % 输出: % out - 结构体, 包含 RR, DET, L, ENTR, LAM, TT if nargin 2, lmin 2; end if nargin 3, vmin 2; end N size(rmat, 1); RR sum(rmat(:)) / N^2; P line_hist(rmat, lmin); l (1:N).; rec_diag sum(l .* P); % 对角线上的总递归点数 det_num sum(l(lmin:end) .* P(lmin:end)); DET det_num / max(rec_diag, eps); L det_num / max(sum(P(lmin:end)), eps); pl P(lmin:end) / max(sum(P(lmin:end)), eps); ENTR -sum(pl .* log(pl eps)); % 取eps防止log(0) V line_hist(rmat., vmin); % 转置后统计垂直结构 v (1:N).; rec_vert sum(v .* V); lam_num sum(v(vmin:end) .* V(vmin:end)); LAM lam_num / max(rec_vert, eps); TT lam_num / max(sum(V(vmin:end)), eps); out struct(RR, RR, DET, DET, L, L, ... ENTR, ENTR, LAM, LAM, TT, TT); end这里有几个实现细节需要解释。第一对角线总递归点数 rec_diag 与 RR 在数学上相等因此 DET 本质上衡量的是“长度不小于 lmin 的递归点占比”。第二由于生成的递归矩阵是对称的统计垂直结构时直接转置后复用 line_hist无需单独写循环。第三eps在 MATLAB 里本身是系统常量这里用来防止分母或 log 参数为 0不改变统计结果。上述口径属于“含主对角线、全矩阵统计”如果要和某篇文献或另一个工具箱对齐必须先确认对方用的是 lmin1 还是 2、是否排除主对角线。3.4 三组示例序列周期、混沌、白噪声一眼看穿把上面三个函数串起来用三段性质完全不同的离散时间序列做验证正弦加噪序列代表周期动力学Logistic 映射代表混沌randn 代表白噪声。% main_demo.m % 用三段离散时间序列演示RQA全流程 N 600; t (0:N-1).; % 1) 周期信号加少量噪声 x1 sin(2*pi*0.05*t) 0.2*randn(N, 1); % 2) Logistic混沌序列 r 3.9; x2 zeros(N, 1); x2(1) 0.4; for i 2:N x2(i) r * x2(i-1) * (1 - x2(i-1)); end % 3) 白噪声 x3 randn(N, 1); m 3; tau 2; % 重构参数 for k 1:3 xk [x1, x2, x3]; % 按列取对应序列 xk xk(:, k); Y rp_reconstruct(xk, m, tau); D pdist2(Y, Y, euclidean); eps_th quantile(D(D 0), 0.3); % 用距离分位数定阈值 rmat rp_matrix(Y, eps_th); out{k} rqa_metrics(rmat, 2, 2); subplot(1, 3, k); imagesc(rmat); colormap(gray); axis square; title(sprintf(RR%.3f DET%.3f, out{k}.RR, out{k}.DET)); end运行之后会得到三种明显不同的递归图周期序列的图上有规则的平行对角线纹理DET 通常较高混沌序列的整体递归率偏低对角线短而碎ENTR 明显高于周期序列白噪声的图几乎找不到成形的结构RR 和 DET 都处于很低水平。数值会随 m、τ、ε 和噪声强度变化但三类序列的相对关系非常稳定这也是用 RQA 做模式判别的前提。4. 参数怎么选嵌入维m、延迟τ和阈值ε的实操口径4.1 时间延迟τ自相关与互信息两种取法τ 的物理含义是两次采样之间的“记忆跨度”。τ 太小状态向量的相邻分量高度相关递归图会被拉成一条粗宽的带子很多结构细节被抹掉τ 太大各分量之间失去关联吸引子被过度展开真实结构反而被噪声覆盖。常见做法有两种自相关函数降到 1/e 处的滞后或者平均互信息第一次出现局部极小值处的滞后。自相关法实现最简单几行代码[ac, lags] xcorr(xk - mean(xk), coeff); zero_lag N; % ac(zero_lag) 对应 lag0 ac_pos abs(ac(zero_lag:end)); % 取正半轴自相关 tau_ac find(ac_pos 1/exp(1), 1) - 1;这段代码先做去均值自相关再找到正半轴上第一个低于 1/e 的位置。1/e 约 0.3679是一个经验阈值低于这个值说明原序列与滞后 N 步的序列相关性已经弱到可以认为“失忆”。互信息法的结果一般比自相关法更可靠因为它刻画的是非线性依赖代价是要先分箱估计联合概率。我自己的使用习惯是先跑自相关快速估一个范围重要数据再用互信息法确认。两者结果如果差别很大多半是数据里有异常段或强趋势需要先清洗。4.2 嵌入维mCao方法和FNN该怎么用嵌入维度太低吸引子会被折叠原本不相邻的轨道在投影中交叉维度太高噪声又被放大。理论上 Takens 定理给出 m ≥ 2d1d 是吸引子维数但实际问题中几乎不会用到这么大的 m。经验上周期信号取 m2 到 3混沌和复杂信号取 m3 到 5高噪声实测数据有时需要取到 5 到 8。工程上常用的定阶方法是虚假最近邻法FNN和 Cao 方法。FNN 的思路是把每个状态向量在 m 维和 m1 维下分别找最近邻如果两者距离之比突然增大说明该近邻在低维下是“虚假”的m 还不够。Cao 方法是对 FNN 的改进它计算 E1(m) 从 m 到 m1 的变化量E1 曲线趋于平缓、E2 接近 1 时就可以停止增加维度。这两个方法在 MATLAB 里都有开源实现但我不建议在初跑 demo 时过度纠结。先用固定 m3、τ2 跑通流程再根据 DET、Lmax 是否随参数平稳来反推合理区间。记住一个耦合关系真正决定重构质量的是嵌入窗长度 (m-1)τ而不是 m 或 τ 单独一个。嵌入窗大致覆盖系统主导周期的四分之一到三分之一是一个比较实用的工程标尺。4.3 阈值ε固定值还是固定递归率阈值 ε 对递归图的影响最直观ε 太小递归图几乎全黑只有零星白点RQA 指标失去统计意义ε 太大图全白所有状态都被当成递归指标完全退化。最常见的做法是按距离矩阵的百分比分位数选阈值比如取 30% 分位让约三成距离被标记为递归点。另一种更可控的做法是固定目标递归率 RR通过二分搜索反推 ε。这在需要和文献结果对齐时特别有用因为不同论文的 RR 取值一致指标才可比。rr_target 0.05; % 目标递归率, 如5% lo 0; hi max(D(:)); while hi - lo 1e-8 mid (lo hi) / 2; rr_now mean(D(:) mid); if rr_now rr_target lo mid; else hi mid; end end eps_rr hi;这个二分搜索的逻辑很朴素每次取中间值递归率高于目标就把上界压低低于目标就把下界抬高。注意每次迭代都在对整个距离矩阵做一次比较Np 超过 2000 时开销很大更高效的办法是随机抽取一部分行做估算。以下是不同场景的阈值选择建议表。数据类型推荐做法备注平稳周期/准周期信号距离矩阵 30% 分位数简单快速递归图纹理清晰混沌信号固定 RR5% 左右便于与文献RQA指标对比高噪声实测信号固定物理阈值或先标准化建议先 z-score 去量纲再定阈值长序列(5000点)随机抽样估算分位数避免全矩阵内存和计算开销另外强调一个习惯无论选哪种阈值计算前最好把序列做标准化减均值除标准差。递归图对量纲很敏感不同物理单位的信号如果用同一个绝对阈值结果完全没有可比性标准化之后至少在同一类信号内部是稳定的。5. RQA分析避坑清单五个最常见的翻车现场5.1 避坑1序列稍长就内存爆掉现象N 在 8000 到 10000 时程序刚跑到距离矩阵阶段MATLAB 直接报内存不足甚至整个编辑器变得卡顿。原因距离矩阵是 Np 乘 Np 的 double 结构Np10000 时单是矩阵就约 800MB还没算后续逻辑矩阵和临时变量。解决使用分块计算加 single 压缩。分块行数降到 200D single(pdist2(...))把距离精度从 double 降到 single递归判定对大多数工程场景足够内存直接减半。blk 200; % 分块调小, 控制内存峰值 rmat false(Np, Np); for i 1:blk:Np idx i:min(iblk-1, Np); D single(pdist2(Y(idx, :), Y, euclidean)); rmat(idx, :) (D eps_th); end另一个思路是在生成距离矩阵之前先对序列做降采样或截取一段关注区间。RQA 对序列长度的需求没有傅里叶分析那么苛刻很多结构信息在两三百个点内就能显现硬啃全序列反而容易得不偿失。5.2 避坑2递归图全黑或全白现象代码能跑完但递归图图像是一整块纯色黑的完全看不到结构白的也看不到结构。原因ε 与序列的实际距离尺度不匹配。原始信号幅值跨度大固定阈值如 0.1 可能远小于所有状态间的距离递归率几乎为零反过来 ε 取 10 又可能把所有距离全部吞掉。解决不要凭感觉给 ε先看距离分布。跑一句histogram(D(D0), 100)观察距离集中在什么量级再用分位数定 ε。对幅值范围很夸张的数据先标准化再选阈值能省掉大半调参时间。5.3 避坑3同一份数据两个工具算的DET对不上现象自己算出的 DET 是 0.62换一个开源工具箱跑同一段序列结果是 0.47两者对不上怀疑代码有 bug。原因RQA 统计口径不统一。最典型的分歧有三个lmin 设为 1 还是 2主对角线是否参与统计统计全矩阵还是只统计上三角。主对角线的长度是 N把它算进分母会显著抬高 DET尤其当 N 不大时影响更明显。解决先定口径再谈数值。我一般在代码里强制固定lmin2、包含主对角线、统计全矩阵并在给别人的结果备注里写明这三个前提。对比文献或工具时先把对方源码里的这些默认值查清楚。5.4 避坑4信号带趋势RQA指标严重失真现象一段缓慢上升的传感器信号递归图左上角和右下角出现大片黑色区域RR 和 DET 都被压得很低和肉眼观察到的明显规律完全不符。原因趋势项在相空间里表现为状态向量持续偏离初始状态早期状态与晚期状态的距离被趋势拉大真实的递归关系被淹没。解决RQA 分析前先做趋势消除。MATLAB 里一行x_detrend detrend(xk, linear)就能去掉线性趋势非线性趋势可以考虑差分处理或分段标准化。处理之后重画递归图结构通常立刻变得清楚。5.5 避坑5对角线统计速度慢到无法接受现象N3000 时 line_hist 已经明显发慢N5000 要跑好几分钟每次调参都像在等进度条。原因line_hist 对每一条偏移对角线都调用 diag 和 find总操作量接近 O(N²)纯 MATLAB 循环的负担很重。解决优先用滑动窗口把长序列切成段每段几百点指标按段计算既能提速又能观察时序变化其次可以考虑把对角线统计改成只处理上三角区域或对列做并行 parfor再不行就寻找带 C-MEX 加速的实现。判断标准很简单N 在 1500 以上优先考虑切段不要硬顶。6. 进阶玩法把RQA指标变成分类特征的三个习惯当你已经能把递归图和 RQA 指标稳定算出来下一步通常是把它们用进识别或预测任务。这里分享三个我长期沿用的习惯。第一个习惯是滑动窗口出“RQA 时间曲线”。固定 m、τ 和阈值把长序列按窗长 W、步长 S 切成若干段每段算一组指标最后得到指标随时间变化的曲线。状态切换、间歇性恶化这类事件在单次全局指标里往往被平均掉但在滑动窗口下会表现为 RR 或 LAM 的显著跳变。注意窗口不宜过小否则每段的递归图太稀疏指标方差会大到失真。W 500; S 50; % 窗长与步长 m 3; tau 2; % 固定重构参数 for k 1:floor((length(xk)-W)/S) 1 seg xk((k-1)*S1 : (k-1)*SW); Y rp_reconstruct(seg, m, tau); rmat rp_matrix(Y, eps_th); feat(k, :) struct2array(rqa_metrics(rmat, 2, 2)); end第二个习惯是做特征选择时先算指标间的相关性。RR 和 DET 天然正相关L 和 Lmax 几乎同涨同跌LAM 和 TT 也高度绑定。如果把全部指标一股脑丢进分类器等于重复放了同一组信息。我用下来比较稳的组合是DET 或 ENTR 反映动力学复杂度LAM 或 TT 反映间歇状态再加一个 RR 反映整体活跃度。三个维度各取其一冗余小、解释性强。第三个习惯是先模拟后实测。每拿到一段真实信号我会先用已知性质的模拟数据——比如正弦加噪、Logistic 混沌、白噪声——在同样的参数口径下跑一遍确认算法没有系统性偏差再处理实测数据。我第一次拿真实信号跑 RQA第一版递归图几乎全黑折腾半天发现是忘了标准化数据幅值跨了两个数量级。现在每次都会先跑一遍最小 demo 把参数口径定下来再去碰实测数据。希望帮到你。本文还有配套的精品资源点击获取