基于MATLAB的去趋势互相关分析(DCCA)算法源码

资 源 简 介

本项目旨在利用MATLAB平台开发一套完整的去趋势互相关分析（Detrended Cross-Correlation Analysis, DCCA）算法工具包。其核心功能是识别和量化两个非平稳时间序列之间的长程幂律互相关性。算法的具体实现流程包括：首先对两个输入的时间序列（如金融时间序列、生理信号等）进行累积求和积分，将其转化为随机游走序列；接着将序列分割为了不同长度n的重叠或非重叠窗口；在每个窗口内，利用最小二乘法对局部趋势进行多项式拟合（如线性、二次或更高阶多项式）并将其从原始数据中去除，以消除非平稳趋势的影响；随后计算两个去趋势序列在当前窗口内的协方差，并对所有窗口求平均以获得q阶（本项目主要针对q=2的情况）波动函数F_DCCA(n)。项目将遍历一系列的时间尺度n，计算对应的波动值。此外，项目还将实现DCCA互相联系数（DCCA Cross-Correlation Coefficient, rho_DCCA）的计算，用于量化两个非平稳信号在不同时间尺度下的具体相关强度，其值介于-1到1之间。最后，程序将生成双对数坐标图（log-log plot），通过线性回归拟合波动函数与时间尺度的关系，提取DCCA标度指数（scaling exponent），从而判断两个系统之间是否存在长程持久性（persistence）或反持久性（anti-persistence）的互相关联。该代码具有良好的封装性，便于用户替换数据源进行金融、气象或生物医学领域的复杂系统分析。

详 情 说 明

项目：基于MATLAB的去趋势互相关分析 (DCCA) 算法实现

项目简介

本项目基于MATLAB平台开发了一套去趋势互相关分析（Detrended Cross-Correlation Analysis, DCCA）算法工具。该工具旨在识别和量化两个非平稳时间序列之间的长程幂律互相关性。通过结合去趋势波动分析（DFA）的原理，本算法能够有效剔除信号中的多项式趋势，揭示隐藏在非平稳数据背后的内在相关行为。

该代码不仅实现了核心的DCCA算法，还涵盖了数据模拟、DFA自相关分析、DCCA互相关系数（$rho_{DCCA}$）计算以及详细的可视化输出，是一个完整的分析流程封装。

功能特性

• 非平稳数据模拟：内置随机游走（Random Walk）序列生成器，可生成具有特定相关性的模拟双变量时间序列，便于算法测试与验证。
• 多尺度波动分析：支持在对数分布的多个时间尺度上计算波动函数，覆盖从微观到宏观的动力学特征。
• 混合分析模式：同时计算两个输入序列的自相关波动（DFA）和它们之间的互相关波动（DCCA）。
• 去趋势处理：采用多项式拟合（默认为线性拟合，即DCCA-1）去除局部趋势，适应非平稳信号环境。
• 互相关系数 $rho_{DCCA}$：计算不同时间尺度下的互相关系数，量化相关性的强弱（范围 [-1, 1]）。
• 标度指数提取：通过双对数坐标下的线性回归，自动提取DCCA标度指数（$alpha_{DCCA}$），用于判断长程持久性。
• 可视化集成：自动生成包含原始序列、双对数波动图、DFA/DCCA对比图及互相关系数变化图的综合面板。

系统要求

• MATLAB R2016a 及以上版本（代码主要利用基础函数，对工具箱无特殊强依赖）。
• 推荐安装 Statistics and Machine Learning Toolbox（用于更高效的统计计算，但基础版亦可运行）。

使用方法

1. 将 main.m 文件放置于MATLAB的当前工作目录下。
2. 直接运行 main 函数。
3. 程序将自动执行以下流程：

* 生成模拟数据。 * 执行DCCA计算。 * 在控制台输出计算得到的标度指数。 * 弹出一个包含4个子图的图形窗口展示分析结果。

注意：若需分析自己的数据，请替换 Step 1 中的 ts1 和 ts2 变量（需确保通过 cumsum 转换为随机游走序列，或直接输入已积分为随机游走的非平稳数据）。

算法实现与逻辑详解

代码 main.m 的执行逻辑严格遵循 DCCA 的标准算法流程，具体步骤如下：

1. 数据准备 (Data Preparation)

• 模拟：首先生成两个长度为 $N=5000$ 的高斯白噪声序列，并通过线性组合引入相关性（系数0.6）。
• 构建：对噪声序列进行累积求和（cumsum），将其转化为非平稳的随机游走时间序列。这是DFA/DCCA分析的标准输入形式。
• 可视化：绘制原始的两个时间序列 $X$ 和 $Y$。

2. 参数设置 (Parameter Setup)

• 多项式阶数：设置 poly_order = 1，表示在去趋势过程中使用线性拟合。
• 时间尺度：定义了最小尺度 (min_scale=10) 和最大尺度 (max_scale=N/4)。
• 对数切分：利用 logspace 生成50个对数均匀分布的尺度 $n$，确保在双对数坐标图上的点分布均匀。

3. DCCA 核心计算 (Core Calculation)

程序首先对输入序列进行预处理：去均值并再次积分，得到用于分析的 Profile 序列。 随后进入主循环，遍历每一个时间尺度 $n$，调用子函数 compute_fluctuation 分别计算：

• $F_{DFA1}$：序列 X 的自相关波动。
• $F_{DFA2}$：序列 Y 的自相关波动。
• $F_{DCCA}$：序列 X 和 Y 的互相关波动。

互相关系数 $rho_{DCCA}(n)$ 计算： 利用公式 $rho(n) = frac{F_{DCCA}^2(n)}{F_{DFA1}(n) times F_{DFA2}(n)}$ 计算互相关系数。代码中对分母为0的情况进行了NaN处理。

4. 标度指数拟合 (Scaling Exponent)

• 对标度 $n$ 和波动函数 $F_{DCCA}$ 取以10为底的对数。
• 利用 polyfit 进行线性回归（一阶拟合）。
• 回归直线的斜率即为 DCCA 标度指数（Scaling Exponent），代码将其输出至控制台。

5. 结果可视化 (Visualization)

代码生成一个 $2 times 2$ 的图窗：

• 子图1：原始输入时间序列展示。
• 子图2：$F_{DCCA}$ 的双对数坐标图（log-log plot）及线性拟合线，标题显示计算出的指数 $alpha$。
• 子图3：同时展示 DFA(X), DFA(Y) 和 DCCA(X,Y) 的波动函数曲线，用于对比不同序列的标度行为。
• 子图4：展示 DCCA 互相关系数 $rho$ 随时间尺度 $n$ 的变化情况，纵轴范围固定在 [-1, 1]。

关键函数分析：compute_fluctuation

这是 main.m 中的核心子函数，负责计算指定尺度 $n$ 下的去趋势波动值。

输入参数：

• seq1, seq2：经过积分处理的 Profile 序列。
• n：当前的窗口长度（时间尺度）。
• order：拟合多项式的阶数。

实现细节与算法逻辑：

1. 窗口分割：将总长度 $N$ 的序列分割为 $lfloor N/n rfloor$ 个互不重叠的窗口（Box）。
2. 局部趋势拟合：

* 在通过循环遍历每个窗口。 * 构建局部时间轴 $t_{local}$。 * 使用 polyfit 对当前窗口内的 seq1 和 seq2 数据分别进行多项式拟合。

1. 去趋势：计算原始数据与拟合多项式值之间的残差（Residuals），即 $detrended = signal - trend$。
2. 协方差计算：

* 如果是 DCCA（即 seq1 不同于 seq2），计算两个残差序列的协方差：$frac{1}{n-1} sum (res_1 cdot res_2)$。 * 如果是 DFA（即 seq1 等于 seq2），该步骤等同于计算方差。

1. 全局平均：累加所有窗口的协方差值，除以窗口数量，得到平均波动 $F^2$。
2. 输出结果：返回 $sqrt{|F^2|}$。

* *注*：代码中对 $F^2$ 取了绝对值 abs(F2) 后再开根号。这是因为在DCCA中，局部协方差的均值可能通过统计波动呈现负值（表示反相关主导），为了保证波动幅度 $F$ 为实数用于对数绘图，通常取其模值。方向性信息保留在 $rho_{DCCA}$ 中。