递归图分析：从MATLAB工具包crptool.zip到非线性时间序列量化实践

简介：CRPTOOL是一个面向非线性动力学与复杂系统研究者的MATLAB交叉复发图（Cross Recurrence Plot, CRP）分析工具箱，专为时间序列同步性、混沌关联性及多变量动态关系建模而设计，适用于生物医学信号处理、气候序列比对、神经网络行为分析等科研场景。压缩包共76个文件，主体为67个MATLAB函数（.m），涵盖数据预处理（normalize.m、taucrp.m）、CRP构建（crp.m、jrp.m）、量化分析（crqad.m、rrspec.m、entropy.m）及GUI交互模块（mgui.m、.gpl.mgui），另含说明文档（crp_man.pdf）、示例数据（logo.mat）和配置文件（info.xml、mgui.rc），整体仅753KB，轻量易部署。已有325人学习下载，资源结构清晰，包含完整工具链：从相空间重构（phasespace.m）、阈值优化（crpclean.m）到可视化渲染（show_crp.m、trackplot.m）及统计指标计算（rpde.m、mi.m），并提供多篇方法注释与典型调用范例（如phasesynchro.m、twinsurr.m），可直接支撑科研复现与教学实践。

1. 项目概述：从crptool.zip到非线性时间序列分析

最近在整理一个老硬盘，翻出来一个名为crptool.zip的压缩包，里面是一套基于MATLAB的非线性时间序列分析工具集。这个工具包的核心围绕着“递归图”展开，包含了think4nn和uppju等模块。对于从事信号处理、复杂系统分析或者金融时间序列研究的朋友来说，这类工具包就像一把瑞士军刀，虽然小众，但在特定场景下能解决大问题。它本质上不是教你如何从零开始写代码，而是提供了一个经过封装、可以直接调用的分析框架，让你能快速地对一维时间序列数据（比如股票价格波动、心率信号、气象数据）进行非线性动力学特征的可视化与量化分析。

如果你手头有一串看起来杂乱无章、但可能蕴含内在规律的数据，想看看它背后是否存在确定性混沌、周期性或者突变行为，那么递归图及其相关的量化分析就是一个非常直观的切入点。这个crptool工具包，就是帮你省去从理论公式到代码实现那漫长且容易出错的过程，直接进入“观察-分析-解读”的环节。接下来，我会结合这个工具包，拆解递归图的核心原理、在MATLAB中的实操流程，以及如何利用think4nn这类模块进行更深层的分析，希望能为相关领域的研究者和工程师提供一份实用的参考指南。

2. 递归图核心原理与工具包设计思路

2.1 什么是递归图？

递归图是一种将时间序列的递归特性可视化为二维二元图像的方法。听起来有点绕，我们可以用一个简单的类比来理解：想象你正在记录一个人每天下午3点的心情得分（1到10分），记录了一个月。递归图要回答的问题是：在这么多天里，是否有某一天的心情状态与另一天“非常相似”？所谓“相似”，是指在多维的状态空间里，这两天的状态向量距离很近。

具体到技术实现，对于一个时间序列{x_i}, i=1...N，我们首先需要通过时间延迟嵌入法重构其相空间。这涉及到两个关键参数：嵌入维度m和时间延迟τ。重构后的相空间中的每个点都是一个向量：Y(i) = [x_i, x_{i+τ}, ..., x_{i+(m-1)τ}]。然后，我们计算所有相空间点对之间的欧氏距离，形成一个距离矩阵。最后，设定一个阈值ε，如果两点之间的距离小于ε，则在递归图矩阵的对应位置(i, j)标记一个点（通常为黑色）。这样生成的二维黑白图像，就是递归图。

注意：阈值ε的选择至关重要。太小，则图中几乎全是空白，无法捕捉递归特性；太大，则图中几乎全黑，失去了分辨能力。通常，ε会选择为相空间点距离分布的一个百分比（例如，使递归密度达到某个预定值），或者是时间序列标准差的某个倍数。

2.2 crptool.zip 工具包的架构解析

crptool.zip作为一个集成工具包，其设计必然遵循了非线性时间序列分析的标准流程。根据其包含的模块名推断，其核心功能模块可能包括：

1. 数据预处理模块：负责对原始时间序列进行去趋势、归一化、去噪等操作，为后续分析准备干净的数据。
2. 相空间重构模块：实现时间延迟τ（常用互信息法第一极小值确定）和嵌入维度m（常用虚假最近邻法确定）的自动或半自动计算。
3. 递归图生成模块：核心模块，根据上述参数和设定的阈值ε，计算并绘制递归图。
4. 递归量化分析模块：这就是think4nn和uppju等模块可能发挥作用的地方。它们不满足于只看图，还要从图中提取定量的特征指标。
   • think4nn：这个模块名强烈暗示了它与“最近邻”思考有关。在递归分析中，除了整体的递归点，对角线结构（表示系统在某个时刻的状态与之后某个时刻的状态相似，即可预测性）和垂直线/水平线结构（表示系统在某个状态停留，即层流或间歇性）是分析重点。think4nn可能专注于分析这些线状结构的长度、分布等，用于量化确定性、预测性等特性。
   • uppju：这个名称较为隐晦，可能是某个特定量化指标（如递归率、确定性、层流度、熵等）的缩写或组合，也可能是某种特定算法的实现。在递归量化分析中，常见的指标包括：
     • 递归率：递归图中黑色点占总点数的比例，反映系统状态总体上的递归概率。
     • 确定性：在对角线方向上形成的线段长度占总递归点数的比例，反映系统的确定性程度。
     • 层流度：垂直或水平方向上形成的线段平均长度，反映系统在某个状态停留的时长。
     • 递归熵：基于递归点分布的香农熵，反映系统动力学的复杂性。uppju很可能就是计算其中某一个或一组这样的指标。
5. 可视化与输出模块：将递归图、量化指标结果以图形和文本形式输出。

这种模块化设计的好处是流程清晰，用户既可以进行全流程分析，也可以单独调用某个模块进行特定计算。

3. 核心细节解析与实操要点

3.1 相空间重构：一切分析的基础

相空间重构的质量直接决定了后续递归图分析的有效性。crptool工具包内应该封装了相应的算法，但理解其原理对正确使用和解读结果至关重要。

• 时间延迟τ的选择：目标是使重构后的各个维度之间既不完全相关（τ太小），也不完全不相关（τ太大）。常用方法是计算时间序列的自互信息函数，并选择其第一个极小值对应的延迟。在MATLAB中，虽然没有内置的互信息函数，但可以基于直方图或核密度估计来实现。工具包可能提供了类似mutual_information(x, maxlag)的函数。
• 嵌入维度m的选择：目标是找到一个足够高的维度，使得重构的相空间能够完全展开动力系统的吸引子，避免不同轨迹的虚假交叉。虚假最近邻法是标准方法。其思想是：随着维度增加，由于投影而显得是“最近邻”的虚假邻居会逐渐消失。当虚假最近邻的比例下降到某个阈值（如5%）以下时，对应的维度就是合适的m。工具包中的think4nn模块很可能就包含了FNN算法的实现。

实操心得： 在自动计算τ和m时，对于噪声较大的数据，互信息函数的第一个极小值可能不明显，FNN曲线可能下降缓慢。此时，盲目相信自动结果可能导致重构不佳。我的经验是：

1. 先对数据进行适当的平滑滤波（如移动平均、小波去噪），但要注意不要过度平滑而抹掉非线性特征。
2. 将自动计算的结果作为参考，手动尝试其附近的一组参数(τ, m)，观察生成的递归图是否具有清晰的结构（如对角线、棋盘格状）。通常，一个“好”的重构参数下，递归图会呈现出相对清晰、有组织的纹理。

3.2 递归图阈值ε的选取策略

阈值ε是递归图的“分辨率旋钮”。crptool可能提供几种设定方式：

1. 固定值：直接给定一个数值。这要求你对数据的尺度有先验知识，不推荐。
2. 标准差倍数：ε = k * σ，其中σ是时间序列的标准差。k通常取0.1到1之间。这是一种常用且稳健的方法，crptool很可能默认采用这种方式。
3. 递归密度固定：指定一个期望的递归密度（如1%，5%，10%），然后反推出对应的ε。这是最科学的方法，因为它使得不同时间序列之间的递归图具有可比性。

实操要点： 在MATLAB中，如果工具包提供了按密度设定阈值的功能，其内部逻辑大致如下：

% 假设已重构相空间点集 Y (N x m 矩阵) distances = []; % 用于存储上三角距离 for i = 1:N-1 for j = i+1:N % 避免重复和零距离 d = norm(Y(i,:) - Y(j,:)); distances = [distances; d]; end end % 对距离排序，找到对应目标密度的分位数 target_density = 0.05; % 目标递归密度 5% num_pairs = N*(N-1)/2; target_num_points = target_density * num_pairs; sorted_dists = sort(distances); epsilon = sorted_dists(round(target_num_points));

使用固定密度法时，建议对同一类数据（如所有的心电图信号）使用相同的密度值，以保证结果的可比性。

4. 实操过程与核心环节实现

4.1 环境准备与数据加载

假设你已经将crptool.zip解压，并将其文件夹路径添加到MATLAB的搜索路径中。我们使用一个经典的混沌时间序列——洛伦兹系统的x分量——作为示例数据。

% 1. 添加工具包路径 addpath(genpath('/你的路径/crptool/')); % 2. 生成或加载示例数据（这里用洛伦兹系统仿真） sigma = 10; rho = 28; beta = 8/3; dt = 0.01; T = 100; % 总时长 steps = floor(T/dt); x = zeros(1, steps); y = zeros(1, steps); z = zeros(1, steps); x(1)=1; y(1)=1; z(1)=1; % 初始值 for i=1:steps-1 dx = sigma*(y(i)-x(i)); dy = x(i)*(rho-z(i))-y(i); dz = x(i)*y(i)-beta*z(i); x(i+1)=x(i)+dx*dt; y(i+1)=y(i)+dy*dt; z(i+1)=z(i)+dz*dt; end % 使用x分量作为分析的时间序列 data = x(1:10:end); % 降采样，使数据点约1000个 time_series = data - mean(data); % 去均值

4.2 全流程分析：从数据到量化指标

根据工具包的设计，可能会有一个主函数来协调整个流程。我们假设这个主函数叫做crp_analysis。

% 3. 调用主分析函数（函数名和参数为假设，需根据实际工具包调整） % 假设函数原型：[RP, metrics] = crp_analysis(data, 'method', 'fnn', 'density', 0.05, 'plot', true); [RP, metrics] = crp_analysis(time_series, ... 'delay_method', 'mutual_info', ... % 延迟选取方法 'dim_method', 'fnn', ... % 维度选取方法 'threshold_method', 'fix_density', ... % 阈值方法 'density', 0.05, ... % 递归密度5% 'normalize', true); % 归一化数据 % RP 是递归图矩阵（二值图像） % metrics 是一个结构体，包含了计算出的各种量化指标 % 4. 绘制递归图 figure; imagesc(1:size(RP,2), 1:size(RP,1), RP); colormap([1 1 1; 0 0 0]); % 黑白配色 axis square; xlabel('Time Index j'); ylabel('Time Index i'); title('Recurrence Plot of Lorenz System (x-component)');

执行上述代码后，你应该能看到一张黑白点阵图。对于混沌的洛伦兹系统，其递归图会呈现出不规则但具有明显纹理的图案，能看到短对角线（短程可预测性）和大量单点或小团块（混沌特性），同时由于系统的有界性，点阵整体分布相对均匀。

4.3 深入量化分析：使用 think4nn 与 uppju 模块

接下来，我们利用工具包中的专门模块进行深入分析。假设think4nn是一个用于分析递归图中对角线结构的函数。

% 5. 分析对角线结构（确定性分析） % 假设函数原型：diag_stats = think4nn(RP, 'min_diag_length', 2); diag_stats = think4nn(RP, 'min_diag_length', 2); disp('--- Diagonal Line Statistics (think4nn) ---'); disp(['Determinism (DET): ', num2str(diag_stats.DET)]); disp(['Average Diagonal Length (L): ', num2str(diag_stats.L)]); disp(['Maximum Diagonal Length (L_max): ', num2str(diag_stats.Lmax)]); disp(['Entropy of Diagonal Lengths (ENTR): ', num2str(diag_stats.ENTR)]);

• 确定性：DET值越高，说明系统的确定性越强，可预测性越好。纯随机噪声的DET会很低。
• 平均对角线长度：L反映了系统平均的可预测时间尺度。
• 对角线长度熵：ENTR反映了对角线长度分布的复杂性，熵值高意味着动力学行为更复杂。

对于uppju模块，我们假设它计算的是与垂直/水平线结构相关的指标，如层流度。

% 6. 分析垂直线结构（层流分析） % 假设函数原型：vert_stats = uppju(RP, 'min_vert_length', 2); vert_stats = uppju(RP, 'min_vert_length', 2); disp('--- Vertical Line Statistics (uppju) ---'); disp(['Laminarity (LAM): ', num2str(vert_stats.LAM)]); disp(['Average Vertical Length (TT): ', num2str(vert_stats.TT)]); disp(['Maximum Vertical Length (V_max): ', num2str(vert_stats.Vmax)]);

• 层流度：LAM表示系统在某个状态“停滞”或缓慢演化的倾向。高LAM可能意味着间歇性行为或状态切换。
• 平均垂直长度：TT即平均 trapping time，量化了系统被困在某个状态的平均时间。

实操现场记录： 在对一段股票收益率序列进行分析时，我发现当市场处于平稳震荡期时，递归图呈现出较密集的短对角线，DET值中等，LAM值也较高，说明存在一定的短期趋势和状态持续性。而在市场暴跌或暴涨的剧烈波动期，递归图变得稀疏且结构破碎，DET和LAM值显著下降，ENTR可能升高，这反映了动力学特性的突变和不可预测性的增加。这种从递归量化指标中捕捉“状态转变”的能力，正是该方法在金融等领域应用的价值所在。

5. 常见问题与排查技巧实录

在实际使用crptool或类似工具包进行递归分析时，你肯定会遇到各种问题。下面是我踩过的一些坑和对应的解决方案。

5.1 数据预处理不当导致分析失效

• 问题现象：递归图一片模糊，几乎全黑或全白，量化指标值异常（如DET接近1或0）。
• 排查思路：
  1. 检查数据平稳性：强烈的趋势会主导距离计算。先对数据进行去趋势处理。可以使用简单的线性拟合去趋势，或更高级的差分、经验模态分解等方法。

     % 示例：线性去趋势 t = 1:length(time_series); p = polyfit(t, time_series, 1); trend = polyval(p, t); detrended_data = time_series - trend;

  2. 检查数据尺度：如果数据绝对值非常大（如股价），计算距离时可能溢出或导致数值问题。务必进行归一化或标准化。

     % Z-score 标准化 normalized_data = (time_series - mean(time_series)) / std(time_series);

  3. 检查噪声水平：过高的噪声会淹没真实的动力学信号。考虑使用平滑滤波器，但需谨慎评估滤波对非线性结构的影响。

5.2 相空间重构参数选择困难

• 问题现象：自动计算的τ或m结果不合理（如τ=1,m非常大或非常小），导致递归图无法揭示任何有意义的结构。
• 排查技巧：
  1. 可视化辅助决策：不要完全依赖自动算法。手动绘制互信息函数和FNN比例随维度变化的曲线，直观判断拐点。

     % 假设工具包提供了 mutual_info_plot 和 fnn_plot 函数 mutual_info_plot(time_series, 50); % 查看前50个延迟的互信息 fnn_plot(time_series, 10); % 查看嵌入维度1到10的FNN比例

  2. 参数扫描：在一个合理的范围内（如τ从1到20，m从2到10），遍历多组参数，生成递归图并观察其纹理变化。选择那个能产生最清晰、最稳定结构的参数组。这虽然计算量大，但对于关键分析是值得的。
  3. 参考领域经验：对于特定类型的数据（如生理信号、气候数据），文献中常有推荐的参数范围起点，可以作为你的初始猜测。

5.3 量化指标解读歧义

• 问题现象：得到了DET,LAM等指标的值，但不知道是高是低，是否显著，如何与不同系统或不同状态进行比较。
• 解决方案：
  1. 建立参考基准：
     • 随机噪声：分析一段相同长度的高斯白噪声序列，得到其指标范围。你的数据指标如果明显超出这个范围，说明检测到了非随机结构。
     • 周期信号：分析一个正弦波，其DET会非常高，且对角线很长。
     • 混沌基准：分析像洛伦兹、罗斯勒这样的标准混沌系统，熟悉其典型指标值。
  2. 使用替代数据检验：这是非线性时间序列分析中的一种重要方法。生成多组与你原始数据具有相同线性特性（如均值、方差、自相关函数）但随机化了相位（从而破坏了非线性结构）的替代数据。分别计算原始数据和所有替代数据的递归量化指标。如果原始数据的指标值落在替代数据指标值分布范围之外（例如，使用95%置信区间），则可以认为检测到了显著的非线性动力学特征。
  3. 关注相对变化，而非绝对值：在许多应用中（如故障监测、状态识别），我们更关心指标随时间或条件的变化趋势。例如，在机械设备振动分析中，当DET下降而ENTR上升时，可能预示着系统从有序运行向混沌故障状态过渡。

5.4 工具包兼容性与报错处理

• 问题：crptool是较老的MATLAB工具包，可能在新版MATLAB（如R2020b以后）中遇到函数兼容性问题，例如某些旧的图形句柄操作或已废弃的函数。
• 排查与修复：
  1. 查看错误信息：MATLAB的命令行窗口会给出具体的错误文件和行号。这是第一线索。
  2. 常见替换：
     • 将findobj('Type', 'figure', ...)等旧的图形对象查找方式，检查是否需要更新。
     • 将str2num替换为更安全的str2double。
     • 检查是否使用了已移除的rand('state', ...)，应改为rng(...)。
  3. 逐函数调试：如果主函数报错，可以尝试将其内部调用的子函数单独拿出来，用示例数据测试，定位具体问题函数。
  4. 社区求助：像crptool这类学术工具包，有时会在MATLAB File Exchange或研究者的个人主页上发布更新版本。可以搜索一下是否有其他人维护的更新版。

最后，我想分享的一点个人体会是，递归图及其量化分析是一个强大的“可视化显微镜”，它能让你“看到”时间序列中隐藏的动力学模式。然而，它也是一把需要精心调校的仪器。参数的选择、数据的预处理、结果的解读，每一步都需要结合具体的物理背景或业务知识进行判断，没有放之四海而皆准的“最佳设置”。最好的学习方式，就是拿一个你熟悉其特性的系统（比如一个简单的正弦波加噪声，一个逻辑斯蒂映射）的数据，用这个工具包从头到尾跑一遍，观察参数变化如何影响最终的图和指标，建立起直观感受。当你对工具的行为有了预期，再用它去探索未知的数据时，才会更有把握，也更容易发现真正有意义的信息。

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