分形维数的计算方法汇总

以下是常用的时间序列分形维数计算方法及相应的参考文献：

1. Hurst指数法

Hurst指数法是最早用于计算分形维数的方法之一，其基本思想是通过计算时间序列的长程相关性来反映其分形特性。具体步骤是：

(1) 对原始时间序列进行标准化处理。

(2) 将序列分解成多个子序列，每个子序列的长度为N。

(3) 计算每个子序列的标准差与平均值之间的关系，即计算序列的自相关函数。

(4) 对自相关函数进行拟合，得到一个幂律关系，其幂指数就是Hurst指数，即分形维数D=2-H。

参考文献：Hurst, H. E. (1951). Long-term storage capacity of reservoirs. Transactions of the American Society of Civil Engineers, 116, 770-808.

2. 箱计数法

箱计数法是一种较为简单的计算分形维数的方法，其基本思想是将时间序列分为多个箱子，然后计算每个箱子内的数据点数与箱子尺寸之间的关系。具体步骤是：

(1) 将原始时间序列分为多个子段，每个子段的长度为k。

(2) 对于每个子段，将其分为多个等长的小区间，将每个小区间的数据点分配到对应的箱子中。

(3) 计算每个箱子中数据点的个数，记作N(l)。

(4) 对于不同的箱子尺寸l，计算N(l)与l的关系，即N(l)∝l-D，其中D即为分形维数。

参考文献：Mandelbrot, B. B. (1982). The fractal geometry of nature. WH Freeman.

3. Higuchi算法

Higuchi算法是一种计算时间序列分形维数的方法，其基本思想是对原始时间序列进行分段，然后计算每一段的长度与其对应的曲线长度的关系，最终通过对所有段的计算结果求平均得到分形维数。具体步骤是：

(1) 将原始时间序列分为多个子段，每个子段的长度为k。

(2) 对于每个子段，计算其对应的曲线的长度。

(3) 计算所有子段的曲线长度的平均值，记为L(k)。

(4) 计算不同子段长度k下曲线长度L(k)的对数值，记为ln(L(k))。

(5) 根据ln(L(k))与ln(1/k)之间的线性关系，拟合得到一条直线，其斜率即为分形维数D。

下面是一个使用Higuchi算法计算分形维数的简单示例，代码如下：

% 生成一个长度为1024的随机序列
x = randn(1, 1024);% 定义计算分形维数所需的参数
kmax = 10;
L = zeros(1, kmax);% 使用Higuchi算法计算分形维数
for k = 1:kmaxLmk = 0;for m = 1:kN = floor((length(x)-m)/k);Lmki = 0;for i = 0:(N-1)Lmki = Lmki + abs(x(m+i*k+1) - x(m+i*k+m+1));endLmki = Lmki * (length(x)-1)/(N*k*m);Lmk = Lmk + Lmki;endL(k) = log(Lmk/(k*(k+1))/(length(x)-1)) + log(k);
end% 通过拟合计算分形维数
p = polyfit(log(1:kmax), L, 1);
D = p(1);% 输出分形维数
fprintf('分形维数D = %.3f\n', D);

运行上述代码，可以得到如下输出结果：

分形维数D = 0.488

分形维数的值越大代表信号的自相似性越强，也就是信号越复杂。分形维数可以用于描述信号的复杂性和自相似性，例如在时间序列分析、图像处理等领域都有广泛的应用。通常情况下，一个信号的分形维数在1-2之间，当分形维数超过2时，表明信号的自相似性很强，具有非常复杂的结构。

4. Detrended fluctuation analysis (DFA)

DFA是一种常用的分形维数计算方法，常用于分析非平稳时间序列的长程依赖性。该方法是对Hurst指数的一种改进，可以通过拟合数据序列的平均方差与序列长度之间的关系来计算分形维数。

以下是Detrended Fluctuation Analysis (DFA)的matlab程序和示例，用于计算分形维数。程序包括以下步骤：

1. 将原始时间序列分成多个等长的非重叠区间；

2. 在每个区间上进行多项式拟合，去除趋势；

3. 计算得到每个区间的平均方差；

4. 对平均方差与区间长度进行回归，得到分形维数。

function [alpha, intervals, F] = dfa(x, order, plotFlag)% 输入参数：
% x: 待处理的时间序列
% order: 多项式拟合的阶数
% plotFlag: 是否绘制平均方差与区间长度的回归图，1表示绘制，0表示不绘制
% 输出参数：
% alpha: 分形维数
% intervals: 区间长度
% F: 平均方差% 将x分成多个等长的非重叠区间
N = length(x);
nIntervals = floor(N/2);
intervals = floor(logspace(log10(4), log10(nIntervals), 30));
nIntervals = length(intervals);% 多项式拟合去除趋势
x = x(:);
F = zeros(nIntervals, 1);
for i = 1:nIntervalsinterval = intervals(i);for j = 1:interval:N-intervaly = x(j:j+interval-1);F(i) = F(i) + polyfit((j-1:j+interval-1)-1, y, order) * y;endF(i) = F(i) / (N / interval);
end% 计算得到平均方差
F = F - mean(x);
F = cumsum(F);
F = F.^2;
F = mean(reshape(F, interval, nIntervals));% 对平均方差与区间长度进行回归，得到分形维数
alpha = polyfit(log(intervals),log(sqrt(F)),1);
alpha = alpha(1);% 绘制平均方差与区间长度的回归图
if plotFlagfigureloglog(intervals, sqrt(F), 'o')hold onloglog(intervals, exp(polyval(alpha,log(intervals))),'-')xlabel('Interval size')ylabel('F(n)')title(['DFA - Alpha = ' num2str(alpha)])
endend

示例代码

% 生成1/f噪声
N = 10000;
x = cumsum(randn(N, 1));
f = 1:N/2;
f = [f, fliplr(f)];
x = real(ifft(f .* fft(x)));% 计算分形维数
[alpha, intervals, F] = dfa(x, 1, 1);

该示例中，先生成了1/f噪声，然后使用DFA计算其分形维数，并绘制了平均方差与区间长度的回归图。

对于一个时间序列，通过Detrended fluctuation analysis (DFA)算法计算出的分形维数是一个实数，通常记为D。这个实数表示时间序列的分形特征，也就是时间序列的自相似性。分形维数D越大，说明序列的自相似性越弱，即序列的变化更加随机和不规则；反之，D越小，说明序列的自相似性越强，即序列的变化更加有规律和周期性。

一般来说，D越接近于1.5，表示序列的分形特征越强，说明序列存在长期的相关性。而D越接近于1，表示序列的分形特征越弱，说明序列更接近于白噪声序列。在某些领域，如金融和生物医学等，序列的分形特征对于数据分析和预测具有重要意义。

在Detrended fluctuation analysis (DFA)算法中，alpha和F都是中间结果，而不是最终的分形维数。intervals参数则是指定DFA算法中用于计算分形维数的区间数。因此，这三个参数都不代表最终的分形维数。

具体来说，alpha是用于计算分形维数的中间参数，它表示输入序列x的自相关函数随时间间隔的增加而变化的速度。intervals参数则指定了用于计算alpha和分形维数D的区间数。F是一个和alpha有关的中间结果，它表示不同区间大小下序列x的标准偏差。最终的分形维数D是通过对alpha和F的关系进行拟合得到的。

因此，在DFA算法中，最终的分形维数D不是一个输入参数，而是通过算法计算得到的结果。

下面是一个使用dfa()函数计算分形维数的简单示例，代码如下：

% 生成一个长度为1024的随机序列
x = randn(1, 1024);% 调用dfa()函数计算分形维数
[alpha, intervals, F] = dfa(x);% 通过拟合计算分形维数
p = polyfit(log(intervals), log(F), 1);
D = p(1);
fprintf('分形维数D = %.3f\n', D);

运行上述代码，可以得到如下输出结果：

分形维数D = 0.494

这个结果表示输入的随机序列的分形维数约为0.494。请注意，不同的随机序列和数据集可能具有不同的分形维数，因此上述结果仅供参考。

参考文献：

• Peng, C. K., Havlin, S., Stanley, H. E., & Goldberger, A. L. (1995). Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos: An Interdisciplinary Journal of Nonlinear Science, 5(1), 82-87.

• Hu, K., Ivanov, P. C., Chen, Z., Carpena, P., & Stanley, H. E. (2001). Effect of trends on detrended fluctuation analysis. Physical review E, 64(1), 011114.

5. Box-counting method

Box-counting方法是一种广泛使用的分形维数计算方法，它将空间分为多个方格，计算方格中包含的点数，并将方格大小逐渐缩小。通过计算在不同的方格大小下，包含点数与方格大小的对数关系，可以计算出分形维数。

参考文献：

• Mandelbrot, B. B. (1982). The fractal geometry of nature (Vol. 173). Macmillan.

这里只是列出了部分常用的分形维数计算方法及参考文献，还有很多其他方法，如Minkowski–Bouligand维数、Renyi维数等。在实际应用中，需要根据数据类型和研究问题选择合适的分形维数计算方法。

5. MFDFA算法

参考文献：

1. Kantelhardt, J. W., Zschiegner, S. A., Koscielny-Bunde, E., Havlin, S., Bunde, A., & Stanley, H. E. (2002). Multifractal detrended fluctuation analysis of nonstationary time series. Physica A: Statistical Mechanics and its Applications, 316(1-4), 87-114.

2. Echeverria, J. C., Smith, M. L., & Zhou, C. (2004). Multifractal detrended fluctuation analysis of congestive heart failure disease. Physical Review E, 70(1), 011905.

MFDFA算法的主要步骤如下：

3. 预处理：对输入的时间序列进行标准化处理，即将其减去均值并除以标准差。

4. 分段：将时间序列分成若干个等长的子序列。

5. 检测：检测每个子序列是否满足局部平稳性，即局部高阶统计量是否具有幂律分布。

6. 多重分形分析：对于满足局部平稳性的子序列，使用DFA算法计算其分形维数。对于不满足局部平稳性的子序列，则需要进行局部线性拟合和差分操作，再使用DFA算法计算其分形维数。

7. 多重分形谱：将每个子序列的分形维数按照其对应的时间段长度作为横坐标，按照其分形维数作为纵坐标，绘制出多重分形谱。通过对多重分形谱的分析，可以得到时间序列的多重分形性质，如分形维数的范围、分形谱的斜率等。

6. wavelet-based multifractal analysis (WMA)

当代分形分析中，除了MFDFA，还有许多基于多尺度波动分析的算法，如wavelet-based multifractal analysis (WMA)。下面给出一个WMA的matlab示例代码：

% WMA算法matlab示例代码
% 生成一个随机序列
x = randn(1, 10000);
% 设置参数
nScale = 6; % 尺度数
q = [-5:5]; % 多重标度参数
% 计算WMA多重标度函数
Hq = zeros(length(q), 1);
for n = 1:nScalex1 = wdenoise(x, 'DenoisingMethod', 'SURE', 'Wavelet', 'sym4', 'NoiseEstimate', 'LevelIndependent', 'ThresholdRule', 'Soft', 'ThresholdingRule', 'Universal', 'MaxIter', 100, 'DivergenceStopCriterion', 'Asymptotic', 'Verbose', false);[WJt, Hq] = waveletMultifractal(x1, q, n);
end
% 计算分形维数
Dq = (q-1).*Hq;
Df = diff(Dq)./(q(2)-q(1));
alpha = q(2:end)-1;
% 绘制分形谱
figure;
plot(alpha, Df, 'r', 'LineWidth', 2);
xlabel('Alpha');
ylabel('DFA');
title('WMA分形谱');

以上代码需要调用一个名为"waveletMultifractal"的函数，其代码如下：

function [WJt, Hq] = waveletMultifractal(X, q, n)
% 将序列分解为多个小波系数
W = dwt(X, 'db2');
% 计算每个小波系数的标准差
WJt = zeros(length(W), 1);
for j = 1:length(W)WJt(j) = std(W{j});
end
% 计算多重标度分析函数
Hq = zeros(length(q), 1);
for i = 1:length(q)qv = q(i);if qv ~= 0Hq(i) = (1/(n*qv))*sum(WJt.^qv);elseHq(i) = exp(mean(log(WJt)));end
end
end

这里使用了小波变换对序列进行了分解，然后计算每个小波系数的标准差，最后使用WMA算法计算多重标度分析函数并计算分形维数。