简介这份MATLAB代码专注于离散时间序列的递归图RP与递归量化分析RQA面向信号处理、复杂系统分析和数据建模领域的研究者与工程师。代码利用递归图可视化时间序列中的相似结构能够识别周期、混沌、分岔等动态特征并计算复发率、最大线段长度、确定性、分层性、熵等RQA指标为评估系统稳定性、复杂性和可预测性提供定量依据。压缩包共含3个文件包括.m格式的MATLAB主程序、一份详细介绍递归分析方法的PDF说明文档以及一个文本说明文件整体仅1.29MB轻量紧凑便于下载后直接开展数据读取、递归图构建、RQA统计量计算和可视化分析。这套代码同样适用于生物医学信号、经济时间序列、物理实验数据等典型离散序列已有415人学习适合需要快速实现完整RQA分析流程、降低算法开发成本的研究人员。1. 用 RQA 量化离散时间序列的递归结构MATLAB 实现从哪里下手拿到一段离散时间序列多数人第一反应是画波形、看频谱或者直接丢给深度网络。但如果你想回答“这段序列内部是否存在重复出现的动力学模式”“系统是否在某一时刻回到了近似相同的状态”频域和时序模型都绕了远路。递归图和递归量化分析Recurrence Quantification Analysis, RQA把高维相空间中的状态访问情况压缩成一张二维图再用一组数值指标描述这张图的结构既不依赖分布假设也不要求序列平稳。对这个标题而言核心任务是把离散时间序列映射为递归图并计算 RQA 指标最后落成一套可以直接运行的 MATLAB 代码。适合做生理信号分析、机械振动状态识别和金融序列模式探测的工程师以及需要快速验证递归特性是否存在的算法研究者。2. 从离散时间序列到递归图相空间重构与递归矩阵的 MATLAB 构建2.1 为什么离散时间序列不能直接画递归图递归图的原始定义建立在系统状态向量上。对于连续系统状态向量来自微分方程对于时间序列我们只能拿到按采样间隔记录的一维观测值。直接把这些观测值当作状态会丢失系统内在的动力学信息。比如一个二阶振荡系统的位置观测值只占状态空间的一条投影曲线真实轨迹是在位置-速度平面上运动的。常见的做法是延迟嵌入Takens 嵌入定理用延迟坐标重构相空间。离散时间序列和连续采样序列在重构上有本质差别如果序列本身来自采样率足够高的物理信号相邻样本高度相关延迟 τ 可以选择较小的值如果序列来自经济数据或事件驱动的符号序列相邻值几乎无记忆那么重构时 τ 的选择就决定了递归结构是否真实存在。先明确这一点后面的 RQA 指标才有意义。2.2 延迟嵌入的最小代码MATLAB 里重构相空间不需要调用任何工具箱核心逻辑是构建一个 Hankel 矩阵。function [XS, nvec] phaseSpaceReconstruct(x, tau, m) % phaseSpaceReconstruct - 延迟嵌入重构相空间 % 输入: % x : 离散时间序列, 列向量 Nx1 % tau: 延迟步数, 正整数 % m : 嵌入维数, 正整数 % 输出: % XS : 相空间状态矩阵, mxL, 每列是一个状态向量 % nvec: 对应的时间索引 N length(x); L N - (m - 1) * tau; % 状态向量个数 XS zeros(m, L); for i 1:m XS(i, :) x(1 (i-1)*tau : L (i-1)*tau); end nvec (1:L); end这段代码的要点是x(1 (i-1)*tau : L (i-1)*tau)把序列按延迟步数错位切片第 i 行保存的是整段序列向左偏移(i-1)*tau个采样点的子序列。矩阵的每一列就是一个 m 维状态向量代表系统在某时刻的重构状态。参数tau和m不能乱选。tau太小会让状态向量的相邻坐标几乎线性相关重构后的轨迹贴近对角线tau太大则让相邻坐标失去关联递归图中噪声主导。m太小无法展开吸引子太小则放大测量噪声。对离散时间序列常见的做法是先根据样本数设定tau的范围再结合互信息函数或自相关函数来筛选。2.3 距离矩阵与递归矩阵的计算有了相空间矩阵 XS任意两个状态向量之间的距离构成了距离矩阵 D。递归矩阵 R 是距离阈值化的结果R(i,j) 1 当 D(i,j) epsilon否则 0% 计算欧氏距离矩阵 D squareform(pdist(XS, euclidean)); % 设定阈值 epsilon maxD max(D(:)); epsilon 0.1 * maxD; % 常用: 最大距离的10% R zeros(size(D)); R(D epsilon) 1;这里用了pdist函数它在 MATLAB 统计工具箱中。如果不希望依赖工具箱可以用两层循环替代但数据量上万时循环速度很慢。中间章的取舍是小规模序列优先可读性大规模序列改用向量化。向量化版本可以用分块计算避免一次性申请大矩阵n size(XS, 2); R false(n, n); blockSize 500; for ib 1:blockSize:n idx ib:min(ibblockSize-1, n); Db sqrt((XS - reshape(XS(:, idx), m, 1, [])) .^ 2 .* ones(1, 1, length(idx))); % 实际为: 计算 XS 与第 idx 块的状态向量两两距离 end这个写法用 reshape 配合隐式扩展避免一次性构造 n×n 的距离矩阵。大量递归图分析时内存瓶颈往往出在这里。2.4 离散序列的符号化处理如果序列是离散事件序列比如点击流、DNA 编码、设备状态码直接用数值做欧氏距离没有物理意义。这时需要对数值距离替换为符号距离D zeros(n, n); for i 1:n for j 1:n if XS(:, i) XS(:, j) D(i, j) 0; else D(i, j) 1; end end end这种汉明距离版本的递归图适合符号动力学场景。阈值 epsilon 也相应变成 0 和 1 之间的整数或取 0。序列类型距离度量阈值设定连续数值观测欧氏距离 / 切比雪夫距离最大距离的 5%15%符号 / 类别序列汉明距离0 或 1电压 / ECoG 等强噪声信号绝对值距离自适应分位数递归矩阵 R 的可视化直接imagesc(R)即可但离散短序列建议把轴改为axis square并关闭坐标轴便于观察拓扑结构。3. RQA 四个核心指标的计算逻辑与 MATLAB 实现3.1 RQA 指标从递归图里读什么递归图是一个二值方阵对角线方向上的连续线段称为对角线结构反映轨迹在某段时间内沿类似路径演化垂直方向上的连续线段垂直线结构反映系统停留在一个状态附近。RQA 便从这两种几何结构中提取数值描述。四个最常用的指标是递归率RR、确定性DET、平均对角线长度L、层流度LAM。RR 是递归点的密度DET 是落在对角线线段上的递归点占全部递归点的比例平均对角线长度表示系统保持相似演化路径的平均时长LAM 是落在垂直线段上的点的比例反映系统“卡住”的程度。3.2 对角线长度统计的向量化实现统计对角线线段长度不能简单对 R 求和需要对每条对角线按连续性分段。朴素实现是遍历每条对角线function histDiag diagHistogram(R, minLength) % diagHistogram - 统计递归图中对角线结构的长度分布 % minLength: 有效线段的最小长度 N size(R, 1); histDiag zeros(N, 1); for diagIdx 1-(N-1):N-1 seg diag(R, diagIdx); if numel(seg) 1 continue; end seg diff([0; seg(:); 0]); starts find(seg 1); ends find(seg -1) - 1; for k 1:numel(starts) len ends(k) - starts(k) 1; if len minLength histDiag(len) histDiag(len) 1; end end end enddiag(R, diagIdx)提取 R 的第 diagIdx 条对角线0 为主对角线正数向右上偏移。diff([0; seg(:); 0])的技巧把连续 1 段的前后边界找出来starts和ends分别对应线段的起点和终点索引。对每条线段计数长度并累加到直方图。参数minLength通常取 2因为长度为 1 的对角线段在噪声信号中大量存在不具动力学含义。3.3 四个指标代码function rqa computeRQA(R, minLine) % computeRQA - 计算递归量化分析指标 % 输入: R - 递归矩阵(0/1二值方阵) % minLine - 对角线/垂直线的最小长度, 通常设为2 % 输出: rqa - 结构体, 字段为各指标 N size(R, 1); % RR 递归率 RR sum(R(:)) / N^2; % 对角线长度统计 histDiag diagHistogram(R, minLine); N_diag_points sum(histDiag .* (1:N)); DET N_diag_points / max(sum(R(:)), eps); L_mean sum(histDiag .* (1:N)) / max(sum(histDiag), eps); % 垂直线统计: 先转置让同一逻辑复用 R_T R; histVert diagHistogram(R_T, minLine); N_vert_points sum(histVert .* (1:N)); LAM N_vert_points / max(sum(R(:)), eps); rqa.RR RR; rqa.DET DET; rqa.meanDiagLen L_mean; rqa.LAM LAM; end垂直线统计利用转置矩阵的“对角线”对应原矩阵的垂直方向这一招让代码量减半。sum(histDiag .* (1:N))计算所有有效线段的长度总和除以线段条数得到平均长度。对离散时间序列一个常见坑是 R 矩阵中存在大量孤立的递归点。这些点会让 DET 虚高或偏低取决于阈值设定的方式。建议先计算sum(R(:))/N^2若 RR 大于 0.2说明阈值过大递归图接近全黑DET 和 LAM 失去区分度。3.4 RQA 参数敏感性与选择策略RQA 指标高度依赖tau、m、epsilon三个参数。对离散时间序列信息论方法选tau是普遍推荐方案。MATLAB 实现互信息函数并不复杂function mi mutualInformation(x, tauMax) % mutualInformation - 计算序列与延迟序列的互信息 % 输出 mi: 长度为 tauMax 的向量 mi zeros(tauMax, 1); edges linspace(min(x), max(x), 16); for tau 1:tauMax pxy hist3([x(1:end-tau), x(1tau:end)], Edges, {edges, edges}); pxy pxy eps; px sum(pxy, 2); py sum(pxy, 1); mi(tau) sum(pxy .* log(pxy ./ (px * py)), all); end endtau取mi曲线第一个局部最小值对应的延迟。嵌入维数用 Cao 方法但考虑到离散序列的样本量通常有限m2或m3在工程中足够。tau和m确定后epsilon的典型选择使 RR 落在 1%10% 之间对离散符号序列直接固定epsilon0即可。4. 离散时间序列 RQA 分析完整流程与噪声信号处理4.1 标准脚本结构数据导入、参数估计、指标输出给出一个可跑的完整脚本框架。注意离散时间序列的预处理和连续信号存在差异这里把预处理单独列出。%% RQA主流程: 从离散序列到指标输出 % 加载数据 x load(your_timeseries.txt); x x(:); % 预处理 x detrend(x, constant); % 移除直流分量 x (x - mean(x)) / std(x); % 标准化 % 选参数 tauMax min(50, floor(length(x)/4)); miVec mutualInformation(x, tauMax); [~, tau] findpeaks(-miVec, NPeaks, 1); % 第一个局部最小 tau tau(1); % 嵌入维数固定为3 (样本量不充裕时) m 3; % 重构相空间 [XS, ~] phaseSpaceReconstruct(x, tau, m); n size(XS, 2); % 距离矩阵与R矩阵 (数据量大时用分块距离) D squareform(pdist(XS, euclidean)); RR_target 0.05; % 目标递归率5% epsilon quantile(D(:), RR_target); % 按递归率分位数确定阈值 R D epsilon; % 计算RQA rqa computeRQA(R, 2); % 显示 fprintf(RR %.3f%%, DET %.3f%%, L_mean %.3f, LAM %.3f\n, ... rqa.RR*100, rqa.DET*100, rqa.meanDiagLen, rqa.LAM);这段代码里的分位数阈值法替代了固定百分比优点是自动适配序列的动态范围不再需要人工尝试多个epsilon。findpeaks用于寻找互信息曲线的局部极小值若找不到峰值则退化为[~, tau] min(miVec)。这里有一个容易踩的坑对纯随机序列互信息曲线没有明显的局部极小findpeaks会返回空数组。所以在实际工程里先验证序列是否具有可递归的动力学特征——可以计算与随机洗牌序列的 RQA 指标对比而不是无条件信任参数搜索结果。4.2 噪声条件下的递归图改善策略离散时间序列常混有测量噪声或量化噪声噪声会以孤立点的形式出现在递归图中。处理方式分两层。第一层对状态向量做局部平均。例如把每个状态向量的邻域均值当作新状态k 3; % 邻域半径 XS_smooth zeros(size(XS)); for i 1:n lo max(1, i-k); hi min(n, ik); XS_smooth(:, i) mean(XS(:, lo:hi), 2); end第二层在阈值选择上使用“容差阈值”即状态向量之间的距离不超过 epsilon 即视为递归。这一层已经在前面的流程里体现。噪声还包括尖峰脉冲脉冲会在递归图上形成十字形结构直接干扰 DET。处理办法是先用中值滤波滤除尖刺再进入 RQA 流程x_filt medfilt1(x, 5);滤波器的窗口大小和采样步长相关这里取 5 是经验值。过滤后丢失的动力学细节有限但 DET 指标稳定性提升明显。4.3 批处理多个序列与滑窗 RQA在实际工程里通常是几十上百段序列同时分析。用 MATLAB 的循环结合前面的函数封装成批处理。另一个高频需求是滑动窗口 RQA用于观测指标随时间的变化winLen 200; step 20; rqaTimeSeries zeros(floor((length(x)-winLen)/step)1, 4); for wi 1:size(rqaTimeSeries, 1) idx (wi-1)*step 1 : (wi-1)*step winLen; xw x(idx); xw detrend(xw, constant); [~, tauW] min(mutualInformation(xw, 30)); m 3; [XS, ~] phaseSpaceReconstruct(xw, tauW, m); Dw squareform(pdist(XS, euclidean)); Rw Dw quantile(Dw(:), 0.05); r computeRQA(Rw, 2); rqaTimeSeries(wi, :) [r.RR, r.DET, r.meanDiagLen, r.LAM]; end滑窗的窗口长度不能太短。离散序列若周期成分明显窗口至少包含两个完整周期否则 RQA 指标会剧烈波动。步长选择不影响结果定性只影响时间的采样密度。quantile按 5% 递归率确定阈值保证了不同窗口间的 RR 一致这对 DET 的纵向比较非常重要。4.4 递归图可视化与结构化解读除了imagesc(R)推荐在图像基础上叠加网格线以辅助目视判断周期长度figure; imagesc(1:n, 1:n, R); colormap([1 1 1; 0 0 0]); axis square; xlabel(State index); ylabel(State index); title(sprintf(Recurrence Plot: RR%.2f%%, rqa.RR*100));当递归图的对角线方向出现周期相等的平行条纹时说明系统存在周期性复现大块黑色方块说明系统在局部状态空间内长时间徘徊常见于混沌系统的层流态。这些纹理判断是 RQA 数字指标的必要补充。5. 递归图滑窗特性在模式切换检测里的验证技巧对离散时间序列做滑窗 RQA 时真正有价值的是指标曲线的转折点。检测模式切换例如脑电从清醒到睡眠、机械设备从正常到故障不一定要训练模型RQA 指标的突变往往已经足够锋利。验证技巧是先对整段序列做一次全局 RQA 作为基准再计算滑窗指标的均值与标准差。当某个滑窗的 DET 超出全局均值 3 倍标准差时标记为候选切换点。这个标准比固定阈值更稳尤其对非平稳的离散序列。一个更细的验证手段是打乱窗口内的序列顺序重新计算 RQA。如果乱序后的 RR 和 DET 明显下降说明原始窗口内存在真实的递归结构如果指标几乎不变说明该窗口实际是噪声段应排除候选点。这个替换检验大约 5 行代码xPerm xw(randperm(length(xw))); [~, tauPerm] min(mutualInformation(xPerm, 30)); [XSP, ~] phaseSpaceReconstruct(xPerm, tauPerm, 3); DP squareform(pdist(XSP, euclidean)); RP DP quantile(DP(:), 0.05); rP computeRQA(RP, 2); if abs(rP.DET - r.DET) 0.05 * r.DET % 差异太小, 标记为噪声窗口 end这里的0.05是经验容差可以根据序列特性收紧到 0.01。验证之后对确认的切换点做标注再配合原始波形观察能快速确认 RQA 的检测结果是否对应到直观可见的模式变化。最后给一个实战建议离散时间序列的采样长度对 RQA 稳定性影响极大。少于 500 个点的短序列RQA 方差大适合用在事件对比而非绝对数值分析超过 5000 个点时直接用完整矩阵计算距离会造成内存溢出优先采用第 2 章给出的分块距离算法。内存节省不是优化的终点换来的是能够处理完整工业数据的自由度。本文还有配套的精品资源点击获取 离散时间序列递归量化分析（RQA）的MATLAB实现指南
简介这份MATLAB代码专注于离散时间序列的递归图RP与递归量化分析RQA面向信号处理、复杂系统分析和数据建模领域的研究者与工程师。代码利用递归图可视化时间序列中的相似结构能够识别周期、混沌、分岔等动态特征并计算复发率、最大线段长度、确定性、分层性、熵等RQA指标为评估系统稳定性、复杂性和可预测性提供定量依据。压缩包共含3个文件包括.m格式的MATLAB主程序、一份详细介绍递归分析方法的PDF说明文档以及一个文本说明文件整体仅1.29MB轻量紧凑便于下载后直接开展数据读取、递归图构建、RQA统计量计算和可视化分析。这套代码同样适用于生物医学信号、经济时间序列、物理实验数据等典型离散序列已有415人学习适合需要快速实现完整RQA分析流程、降低算法开发成本的研究人员。1. 用 RQA 量化离散时间序列的递归结构MATLAB 实现从哪里下手拿到一段离散时间序列多数人第一反应是画波形、看频谱或者直接丢给深度网络。但如果你想回答“这段序列内部是否存在重复出现的动力学模式”“系统是否在某一时刻回到了近似相同的状态”频域和时序模型都绕了远路。递归图和递归量化分析Recurrence Quantification Analysis, RQA把高维相空间中的状态访问情况压缩成一张二维图再用一组数值指标描述这张图的结构既不依赖分布假设也不要求序列平稳。对这个标题而言核心任务是把离散时间序列映射为递归图并计算 RQA 指标最后落成一套可以直接运行的 MATLAB 代码。适合做生理信号分析、机械振动状态识别和金融序列模式探测的工程师以及需要快速验证递归特性是否存在的算法研究者。2. 从离散时间序列到递归图相空间重构与递归矩阵的 MATLAB 构建2.1 为什么离散时间序列不能直接画递归图递归图的原始定义建立在系统状态向量上。对于连续系统状态向量来自微分方程对于时间序列我们只能拿到按采样间隔记录的一维观测值。直接把这些观测值当作状态会丢失系统内在的动力学信息。比如一个二阶振荡系统的位置观测值只占状态空间的一条投影曲线真实轨迹是在位置-速度平面上运动的。常见的做法是延迟嵌入Takens 嵌入定理用延迟坐标重构相空间。离散时间序列和连续采样序列在重构上有本质差别如果序列本身来自采样率足够高的物理信号相邻样本高度相关延迟 τ 可以选择较小的值如果序列来自经济数据或事件驱动的符号序列相邻值几乎无记忆那么重构时 τ 的选择就决定了递归结构是否真实存在。先明确这一点后面的 RQA 指标才有意义。2.2 延迟嵌入的最小代码MATLAB 里重构相空间不需要调用任何工具箱核心逻辑是构建一个 Hankel 矩阵。function [XS, nvec] phaseSpaceReconstruct(x, tau, m) % phaseSpaceReconstruct - 延迟嵌入重构相空间 % 输入: % x : 离散时间序列, 列向量 Nx1 % tau: 延迟步数, 正整数 % m : 嵌入维数, 正整数 % 输出: % XS : 相空间状态矩阵, mxL, 每列是一个状态向量 % nvec: 对应的时间索引 N length(x); L N - (m - 1) * tau; % 状态向量个数 XS zeros(m, L); for i 1:m XS(i, :) x(1 (i-1)*tau : L (i-1)*tau); end nvec (1:L); end这段代码的要点是x(1 (i-1)*tau : L (i-1)*tau)把序列按延迟步数错位切片第 i 行保存的是整段序列向左偏移(i-1)*tau个采样点的子序列。矩阵的每一列就是一个 m 维状态向量代表系统在某时刻的重构状态。参数tau和m不能乱选。tau太小会让状态向量的相邻坐标几乎线性相关重构后的轨迹贴近对角线tau太大则让相邻坐标失去关联递归图中噪声主导。m太小无法展开吸引子太小则放大测量噪声。对离散时间序列常见的做法是先根据样本数设定tau的范围再结合互信息函数或自相关函数来筛选。2.3 距离矩阵与递归矩阵的计算有了相空间矩阵 XS任意两个状态向量之间的距离构成了距离矩阵 D。递归矩阵 R 是距离阈值化的结果R(i,j) 1 当 D(i,j) epsilon否则 0% 计算欧氏距离矩阵 D squareform(pdist(XS, euclidean)); % 设定阈值 epsilon maxD max(D(:)); epsilon 0.1 * maxD; % 常用: 最大距离的10% R zeros(size(D)); R(D epsilon) 1;这里用了pdist函数它在 MATLAB 统计工具箱中。如果不希望依赖工具箱可以用两层循环替代但数据量上万时循环速度很慢。中间章的取舍是小规模序列优先可读性大规模序列改用向量化。向量化版本可以用分块计算避免一次性申请大矩阵n size(XS, 2); R false(n, n); blockSize 500; for ib 1:blockSize:n idx ib:min(ibblockSize-1, n); Db sqrt((XS - reshape(XS(:, idx), m, 1, [])) .^ 2 .* ones(1, 1, length(idx))); % 实际为: 计算 XS 与第 idx 块的状态向量两两距离 end这个写法用 reshape 配合隐式扩展避免一次性构造 n×n 的距离矩阵。大量递归图分析时内存瓶颈往往出在这里。2.4 离散序列的符号化处理如果序列是离散事件序列比如点击流、DNA 编码、设备状态码直接用数值做欧氏距离没有物理意义。这时需要对数值距离替换为符号距离D zeros(n, n); for i 1:n for j 1:n if XS(:, i) XS(:, j) D(i, j) 0; else D(i, j) 1; end end end这种汉明距离版本的递归图适合符号动力学场景。阈值 epsilon 也相应变成 0 和 1 之间的整数或取 0。序列类型距离度量阈值设定连续数值观测欧氏距离 / 切比雪夫距离最大距离的 5%15%符号 / 类别序列汉明距离0 或 1电压 / ECoG 等强噪声信号绝对值距离自适应分位数递归矩阵 R 的可视化直接imagesc(R)即可但离散短序列建议把轴改为axis square并关闭坐标轴便于观察拓扑结构。3. RQA 四个核心指标的计算逻辑与 MATLAB 实现3.1 RQA 指标从递归图里读什么递归图是一个二值方阵对角线方向上的连续线段称为对角线结构反映轨迹在某段时间内沿类似路径演化垂直方向上的连续线段垂直线结构反映系统停留在一个状态附近。RQA 便从这两种几何结构中提取数值描述。四个最常用的指标是递归率RR、确定性DET、平均对角线长度L、层流度LAM。RR 是递归点的密度DET 是落在对角线线段上的递归点占全部递归点的比例平均对角线长度表示系统保持相似演化路径的平均时长LAM 是落在垂直线段上的点的比例反映系统“卡住”的程度。3.2 对角线长度统计的向量化实现统计对角线线段长度不能简单对 R 求和需要对每条对角线按连续性分段。朴素实现是遍历每条对角线function histDiag diagHistogram(R, minLength) % diagHistogram - 统计递归图中对角线结构的长度分布 % minLength: 有效线段的最小长度 N size(R, 1); histDiag zeros(N, 1); for diagIdx 1-(N-1):N-1 seg diag(R, diagIdx); if numel(seg) 1 continue; end seg diff([0; seg(:); 0]); starts find(seg 1); ends find(seg -1) - 1; for k 1:numel(starts) len ends(k) - starts(k) 1; if len minLength histDiag(len) histDiag(len) 1; end end end enddiag(R, diagIdx)提取 R 的第 diagIdx 条对角线0 为主对角线正数向右上偏移。diff([0; seg(:); 0])的技巧把连续 1 段的前后边界找出来starts和ends分别对应线段的起点和终点索引。对每条线段计数长度并累加到直方图。参数minLength通常取 2因为长度为 1 的对角线段在噪声信号中大量存在不具动力学含义。3.3 四个指标代码function rqa computeRQA(R, minLine) % computeRQA - 计算递归量化分析指标 % 输入: R - 递归矩阵(0/1二值方阵) % minLine - 对角线/垂直线的最小长度, 通常设为2 % 输出: rqa - 结构体, 字段为各指标 N size(R, 1); % RR 递归率 RR sum(R(:)) / N^2; % 对角线长度统计 histDiag diagHistogram(R, minLine); N_diag_points sum(histDiag .* (1:N)); DET N_diag_points / max(sum(R(:)), eps); L_mean sum(histDiag .* (1:N)) / max(sum(histDiag), eps); % 垂直线统计: 先转置让同一逻辑复用 R_T R; histVert diagHistogram(R_T, minLine); N_vert_points sum(histVert .* (1:N)); LAM N_vert_points / max(sum(R(:)), eps); rqa.RR RR; rqa.DET DET; rqa.meanDiagLen L_mean; rqa.LAM LAM; end垂直线统计利用转置矩阵的“对角线”对应原矩阵的垂直方向这一招让代码量减半。sum(histDiag .* (1:N))计算所有有效线段的长度总和除以线段条数得到平均长度。对离散时间序列一个常见坑是 R 矩阵中存在大量孤立的递归点。这些点会让 DET 虚高或偏低取决于阈值设定的方式。建议先计算sum(R(:))/N^2若 RR 大于 0.2说明阈值过大递归图接近全黑DET 和 LAM 失去区分度。3.4 RQA 参数敏感性与选择策略RQA 指标高度依赖tau、m、epsilon三个参数。对离散时间序列信息论方法选tau是普遍推荐方案。MATLAB 实现互信息函数并不复杂function mi mutualInformation(x, tauMax) % mutualInformation - 计算序列与延迟序列的互信息 % 输出 mi: 长度为 tauMax 的向量 mi zeros(tauMax, 1); edges linspace(min(x), max(x), 16); for tau 1:tauMax pxy hist3([x(1:end-tau), x(1tau:end)], Edges, {edges, edges}); pxy pxy eps; px sum(pxy, 2); py sum(pxy, 1); mi(tau) sum(pxy .* log(pxy ./ (px * py)), all); end endtau取mi曲线第一个局部最小值对应的延迟。嵌入维数用 Cao 方法但考虑到离散序列的样本量通常有限m2或m3在工程中足够。tau和m确定后epsilon的典型选择使 RR 落在 1%10% 之间对离散符号序列直接固定epsilon0即可。4. 离散时间序列 RQA 分析完整流程与噪声信号处理4.1 标准脚本结构数据导入、参数估计、指标输出给出一个可跑的完整脚本框架。注意离散时间序列的预处理和连续信号存在差异这里把预处理单独列出。%% RQA主流程: 从离散序列到指标输出 % 加载数据 x load(your_timeseries.txt); x x(:); % 预处理 x detrend(x, constant); % 移除直流分量 x (x - mean(x)) / std(x); % 标准化 % 选参数 tauMax min(50, floor(length(x)/4)); miVec mutualInformation(x, tauMax); [~, tau] findpeaks(-miVec, NPeaks, 1); % 第一个局部最小 tau tau(1); % 嵌入维数固定为3 (样本量不充裕时) m 3; % 重构相空间 [XS, ~] phaseSpaceReconstruct(x, tau, m); n size(XS, 2); % 距离矩阵与R矩阵 (数据量大时用分块距离) D squareform(pdist(XS, euclidean)); RR_target 0.05; % 目标递归率5% epsilon quantile(D(:), RR_target); % 按递归率分位数确定阈值 R D epsilon; % 计算RQA rqa computeRQA(R, 2); % 显示 fprintf(RR %.3f%%, DET %.3f%%, L_mean %.3f, LAM %.3f\n, ... rqa.RR*100, rqa.DET*100, rqa.meanDiagLen, rqa.LAM);这段代码里的分位数阈值法替代了固定百分比优点是自动适配序列的动态范围不再需要人工尝试多个epsilon。findpeaks用于寻找互信息曲线的局部极小值若找不到峰值则退化为[~, tau] min(miVec)。这里有一个容易踩的坑对纯随机序列互信息曲线没有明显的局部极小findpeaks会返回空数组。所以在实际工程里先验证序列是否具有可递归的动力学特征——可以计算与随机洗牌序列的 RQA 指标对比而不是无条件信任参数搜索结果。4.2 噪声条件下的递归图改善策略离散时间序列常混有测量噪声或量化噪声噪声会以孤立点的形式出现在递归图中。处理方式分两层。第一层对状态向量做局部平均。例如把每个状态向量的邻域均值当作新状态k 3; % 邻域半径 XS_smooth zeros(size(XS)); for i 1:n lo max(1, i-k); hi min(n, ik); XS_smooth(:, i) mean(XS(:, lo:hi), 2); end第二层在阈值选择上使用“容差阈值”即状态向量之间的距离不超过 epsilon 即视为递归。这一层已经在前面的流程里体现。噪声还包括尖峰脉冲脉冲会在递归图上形成十字形结构直接干扰 DET。处理办法是先用中值滤波滤除尖刺再进入 RQA 流程x_filt medfilt1(x, 5);滤波器的窗口大小和采样步长相关这里取 5 是经验值。过滤后丢失的动力学细节有限但 DET 指标稳定性提升明显。4.3 批处理多个序列与滑窗 RQA在实际工程里通常是几十上百段序列同时分析。用 MATLAB 的循环结合前面的函数封装成批处理。另一个高频需求是滑动窗口 RQA用于观测指标随时间的变化winLen 200; step 20; rqaTimeSeries zeros(floor((length(x)-winLen)/step)1, 4); for wi 1:size(rqaTimeSeries, 1) idx (wi-1)*step 1 : (wi-1)*step winLen; xw x(idx); xw detrend(xw, constant); [~, tauW] min(mutualInformation(xw, 30)); m 3; [XS, ~] phaseSpaceReconstruct(xw, tauW, m); Dw squareform(pdist(XS, euclidean)); Rw Dw quantile(Dw(:), 0.05); r computeRQA(Rw, 2); rqaTimeSeries(wi, :) [r.RR, r.DET, r.meanDiagLen, r.LAM]; end滑窗的窗口长度不能太短。离散序列若周期成分明显窗口至少包含两个完整周期否则 RQA 指标会剧烈波动。步长选择不影响结果定性只影响时间的采样密度。quantile按 5% 递归率确定阈值保证了不同窗口间的 RR 一致这对 DET 的纵向比较非常重要。4.4 递归图可视化与结构化解读除了imagesc(R)推荐在图像基础上叠加网格线以辅助目视判断周期长度figure; imagesc(1:n, 1:n, R); colormap([1 1 1; 0 0 0]); axis square; xlabel(State index); ylabel(State index); title(sprintf(Recurrence Plot: RR%.2f%%, rqa.RR*100));当递归图的对角线方向出现周期相等的平行条纹时说明系统存在周期性复现大块黑色方块说明系统在局部状态空间内长时间徘徊常见于混沌系统的层流态。这些纹理判断是 RQA 数字指标的必要补充。5. 递归图滑窗特性在模式切换检测里的验证技巧对离散时间序列做滑窗 RQA 时真正有价值的是指标曲线的转折点。检测模式切换例如脑电从清醒到睡眠、机械设备从正常到故障不一定要训练模型RQA 指标的突变往往已经足够锋利。验证技巧是先对整段序列做一次全局 RQA 作为基准再计算滑窗指标的均值与标准差。当某个滑窗的 DET 超出全局均值 3 倍标准差时标记为候选切换点。这个标准比固定阈值更稳尤其对非平稳的离散序列。一个更细的验证手段是打乱窗口内的序列顺序重新计算 RQA。如果乱序后的 RR 和 DET 明显下降说明原始窗口内存在真实的递归结构如果指标几乎不变说明该窗口实际是噪声段应排除候选点。这个替换检验大约 5 行代码xPerm xw(randperm(length(xw))); [~, tauPerm] min(mutualInformation(xPerm, 30)); [XSP, ~] phaseSpaceReconstruct(xPerm, tauPerm, 3); DP squareform(pdist(XSP, euclidean)); RP DP quantile(DP(:), 0.05); rP computeRQA(RP, 2); if abs(rP.DET - r.DET) 0.05 * r.DET % 差异太小, 标记为噪声窗口 end这里的0.05是经验容差可以根据序列特性收紧到 0.01。验证之后对确认的切换点做标注再配合原始波形观察能快速确认 RQA 的检测结果是否对应到直观可见的模式变化。最后给一个实战建议离散时间序列的采样长度对 RQA 稳定性影响极大。少于 500 个点的短序列RQA 方差大适合用在事件对比而非绝对数值分析超过 5000 个点时直接用完整矩阵计算距离会造成内存溢出优先采用第 2 章给出的分块距离算法。内存节省不是优化的终点换来的是能够处理完整工业数据的自由度。本文还有配套的精品资源点击获取