混沌时间序列参数计算工具箱

混沌时间序列参数计算工具箱，包含时间延迟、嵌入维数、关联维等关键参数的计算程序。

混沌时间序列参数计算工具箱

1. 主函数 - 混沌参数综合分析

function chaos_analysis_results = chaotic_time_series_analysis(time_series, varargin)
    % 混沌时间序列参数综合分析
    % 输入:
    %   time_series - 时间序列数据
    %   varargin - 可选参数
    % 输出:
    %   chaos_analysis_results - 包含所有分析结果的结构体
p = inputParser;

addRequired(p, 'time_series', @isnumeric);

addParameter(p, 'max_delay', 50, @isnumeric);

addParameter(p, 'max_dim', 10, @isnumeric);

addParameter(p, 'radius_range', logspace(-3, 0, 50), @isnumeric);

addParameter(p, 'show_plots', true, @islogical);
parse(p, time_series, varargin{:});
fprintf('开始混沌时间序列分析...\n');

fprintf('数据长度: %d\n', length(time_series));
% 初始化结果结构

results = struct();
%% 1. 时间延迟计算

fprintf('计算时间延迟...\n');

[results.tau_autocorr, results.autocorr_func] = calculate_time_delay_autocorr(time_series, p.Results.max_delay);

[results.tau_mutual, results.mutual_info] = calculate_time_delay_mutual(time_series, p.Results.max_delay);
%% 2. 嵌入维数计算

fprintf('计算嵌入维数...\n');

[results.embedding_dim, results.fnn_ratio] = calculate_embedding_dimension(time_series, results.tau_mutual, p.Results.max_dim);
%% 3. 关联维计算

fprintf('计算关联维...\n');

[results.correlation_dim, results.correlation_integral] = calculate_correlation_dimension(time_series, results.tau_mutual, results.embedding_dim, p.Results.radius_range);
%% 4. Lyapunov指数计算

fprintf('计算最大Lyapunov指数...\n');

[results.lyapunov_exponent, results.lyapunov_curve] = calculate_lyapunov_exponent(time_series, results.tau_mutual, results.embedding_dim);
%% 5. 其他混沌指标

fprintf('计算其他混沌指标...\n');

results.hurst_exponent = calculate_hurst_exponent(time_series);

[results.approximate_entropy, results.sample_entropy] = calculate_entropy_measures(time_series);
%% 6. 可视化结果

if p.Results.show_plots

plot_chaos_analysis_results(results, time_series);

end
%% 7. 生成报告

generate_analysis_report(results);
chaos_analysis_results = results;

fprintf('混沌分析完成!\n');

end

2. 时间延迟计算

2.1 自相关函数法

function [optimal_tau, autocorr_values] = calculate_time_delay_autocorr(time_series, max_delay)
    % 使用自相关函数法计算时间延迟
    % 选择自相关函数第一次过零或下降到1-1/e时的延迟
n = length(time_series);

autocorr_values = zeros(1, max_delay);
% 计算自相关函数

for tau = 1:max_delay

autocorr_values(tau) = abs(corr(time_series(1:end-tau)', time_series(1+tau:end)'));

end
% 寻找最优延迟

threshold = 1 - 1/exp(1); % 约0.6321
% 方法1: 第一次低于阈值

below_threshold = find(autocorr_values < threshold, 1);

if ~isempty(below_threshold)

optimal_tau = below_threshold;

else

% 方法2: 选择最小值

[~, optimal_tau] = min(autocorr_values);

end
fprintf('自相关函数法 - 最优时间延迟: %d\n', optimal_tau);

end

2.2 互信息法

function [optimal_tau, mutual_info_values] = calculate_time_delay_mutual(time_series, max_delay)
    % 使用互信息法计算时间延迟
    % 选择互信息第一次达到局部最小值时的延迟
n = length(time_series);

mutual_info_values = zeros(1, max_delay);
% 将时间序列离散化

num_bins = 10; % 直方图箱数

[~, edges] = histcounts(time_series, num_bins);
% 计算原始序列的概率分布

px = histcounts(time_series, edges, 'Normalization', 'probability');
for tau = 1:max_delay

% 延迟序列

x_t = time_series(1:end-tau);

x_t_tau = time_series(1+tau:end);
% 计算联合概率分布

joint_hist = histcounts2(x_t, x_t_tau, edges, edges, 'Normalization', 'probability');

py = histcounts(x_t_tau, edges, 'Normalization', 'probability');
% 计算互信息

mi = 0;

for i = 1:num_bins

for j = 1:num_bins

if joint_hist(i,j) > 0 && px(i) > 0 && py(j) > 0

mi = mi + joint_hist(i,j) * log(joint_hist(i,j) / (px(i) * py(j)));

end

end

end
mutual_info_values(tau) = mi;

end
% 寻找第一个局部最小值

optimal_tau = find_first_local_minimum(mutual_info_values);
if isempty(optimal_tau)

optimal_tau = 1; % 默认值

end
fprintf('互信息法 - 最优时间延迟: %d\n', optimal_tau);

end
function min_idx = find_first_local_minimum(data)

% 寻找第一个局部最小值

for i = 2:length(data)-1

if data(i) < data(i-1) && data(i) < data(i+1)

min_idx = i;

return;

end

end

min_idx = [];

end

3. 嵌入维数计算

function [optimal_dim, fnn_ratio] = calculate_embedding_dimension(time_series, tau, max_dim)
    % 使用伪最近邻法计算嵌入维数
    % 当伪最近邻比例低于阈值时的最小维数即为最优嵌入维数
n = length(time_series);

fnn_ratio = zeros(1, max_dim);
% 计算每个嵌入维数的伪最近邻比例

for m = 1:max_dim

% 重构相空间

embedded_data = phase_space_reconstruction(time_series, m, tau);
% 寻找最近邻

[fnn_ratio(m), ~] = calculate_fnn_ratio(embedded_data, m);
fprintf('嵌入维数 %d: 伪最近邻比例 = %.4f\n', m, fnn_ratio(m));

end
% 选择伪最近邻比例低于阈值的最小维数

threshold = 0.1; % 常用阈值

below_threshold = find(fnn_ratio < threshold, 1);
if ~isempty(below_threshold)

optimal_dim = below_threshold;

else

% 选择变化平缓的点

diff_ratio = diff(fnn_ratio);

if any(diff_ratio < 0.01)

optimal_dim = find(diff_ratio < 0.01, 1) + 1;

else

optimal_dim = max_dim;

end

end
fprintf('最优嵌入维数: %d\n', optimal_dim);

end
function embedded_data = phase_space_reconstruction(time_series, m, tau)

% 相空间重构

n = length(time_series);

embedded_length = n - (m-1)*tau;
embedded_data = zeros(embedded_length, m);
for i = 1:embedded_length

for j = 1:m

embedded_data(i, j) = time_series(i + (j-1)*tau);

end

end

end
function [fnn_ratio, neighbors] = calculate_fnn_ratio(embedded_data, m)

% 计算伪最近邻比例
[n_points, ~] = size(embedded_data);
% 寻找每个点的最近邻

distances = zeros(n_points, 1);

neighbor_indices = zeros(n_points, 1);
for i = 1:n_points

% 计算到其他点的距离（排除自身）

dist = zeros(n_points, 1);

for j = 1:n_points

if i ~= j

dist(j) = norm(embedded_data(i, 😃 - embedded_data(j, 😃);

else

dist(j) = inf;

end

end
[min_dist, min_idx] = min(dist);

distances(i) = min_dist;

neighbor_indices(i) = min_idx;

end
% 计算伪最近邻比例

fnn_count = 0;

threshold = 10; % 常用阈值
for i = 1:n_points

current_point = embedded_data(i, 😃;

neighbor_point = embedded_data(neighbor_indices(i), 😃;
% 计算距离比

if m > 1

% 在m-1维空间中的距离

dist_m_minus_1 = norm(current_point(1:end-1) - neighbor_point(1:end-1));

dist_m = distances(i);
if dist_m_minus_1 > 0

ratio = abs(dist_m - dist_m_minus_1) / dist_m_minus_1;

if ratio > threshold

fnn_count = fnn_count + 1;

end

end

end

end
fnn_ratio = fnn_count / n_points;

neighbors = neighbor_indices;

end

4. 关联维计算

function [correlation_dim, corr_integral] = calculate_correlation_dimension(time_series, tau, m, radius_range)
    % 使用G-P算法计算关联维
% 重构相空间

embedded_data = phase_space_reconstruction(time_series, m, tau);

[n_points, ~] = size(embedded_data);
corr_integral = zeros(length(radius_range), 1);
% 计算关联积分

for r_idx = 1:length(radius_range)

r = radius_range(r_idx);

count = 0;
% 使用随机采样来加速计算

sample_size = min(1000, n_points);

sample_indices = randperm(n_points, sample_size);
for i = 1:sample_size

idx_i = sample_indices(i);

for j = i+1:sample_size

idx_j = sample_indices(j);

distance_ij = norm(embedded_data(idx_i, 😃 - embedded_data(idx_j, 😃);

if distance_ij < r

count = count + 1;

end

end

end
% 计算关联积分

corr_integral(r_idx) = 2 * count / (sample_size * (sample_size - 1));

end
% 寻找线性区域并计算关联维

[correlation_dim, scaling_region] = estimate_correlation_dimension(radius_range, corr_integral);
fprintf('关联维数: %.4f\n', correlation_dim);

end
function [correlation_dim, scaling_region] = estimate_correlation_dimension(radius, corr_integral)

% 估计关联维数 - 通过线性区域的斜率
% 对数变换

log_r = log10(radius);

log_c = log10(corr_integral);
% 移除无效点（零或负值）

valid_idx = isfinite(log_r) & isfinite(log_c) & log_c > -inf;

log_r_valid = log_r(valid_idx);

log_c_valid = log_c(valid_idx);
if length(log_r_valid) < 10

correlation_dim = 0;

scaling_region = [];

return;

end
% 寻找线性区域（通过滑动窗口）

window_size = 10;

max_slope = -inf;

best_start = 1;
for i = 1:length(log_r_valid)-window_size+1

window_r = log_r_valid(i:i+window_size-1);

window_c = log_c_valid(i:i+window_size-1);
% 线性拟合

p = polyfit(window_r, window_c, 1);

slope = p(1);
% 检查拟合质量

residuals = window_c - polyval(p, window_r);

r_squared = 1 - sum(residuals.^2) / sum((window_c - mean(window_c)).^2);
if slope > max_slope && r_squared > 0.9

max_slope = slope;

best_start = i;

end

end
% 使用最佳线性区域计算关联维

scaling_region = best_start:best_start+window_size-1;

p_final = polyfit(log_r_valid(scaling_region), log_c_valid(scaling_region), 1);

correlation_dim = p_final(1);

end

5. Lyapunov指数计算

function [lyapunov_exp, divergence_curve] = calculate_lyapunov_exponent(time_series, tau, m)
    % 计算最大Lyapunov指数 - 使用小数据量法
% 重构相空间

embedded_data = phase_space_reconstruction(time_series, m, tau);

[n_points, ~] = size(embedded_data);
% 参数设置

min_neighbors = 5;

evolution_steps = min(50, floor(n_points/10));
% 为每个点寻找最近邻

divergence = zeros(evolution_steps, 1);

count = zeros(evolution_steps, 1);
for i = 1:n_points-evolution_steps

% 当前点

current_point = embedded_data(i, 😃;
% 寻找最近邻（排除太近的点）

distances = zeros(n_points, 1);

for j = 1:n_points

if abs(i - j) > evolution_steps % 避免时间相关

distances(j) = norm(current_point - embedded_data(j, 😃);

else

distances(j) = inf;

end

end
[~, neighbor_indices] = mink(distances, min_neighbors);
% 追踪每个最近邻的演化

for k = 1:min_neighbors

neighbor_idx = neighbor_indices(k);

initial_distance = norm(embedded_data(i, 😃 - embedded_data(neighbor_idx, 😃);
if initial_distance > 0

for step = 1:min(evolution_steps, n_points-max(i, neighbor_idx))

d_i = i + step;

d_j = neighbor_idx + step;
if d_j <= n_points

current_distance = norm(embedded_data(d_i, 😃 - embedded_data(d_j, 😃);

if current_distance > 0

divergence(step) = divergence(step) + log(current_distance / initial_distance);

count(step) = count(step) + 1;

end

end

end

end

end

end
% 计算平均发散率

valid_steps = count > 0;

divergence_curve = zeros(size(divergence));

divergence_curve(valid_steps) = divergence(valid_steps) ./ count(valid_steps);
% 线性拟合得到Lyapunov指数

time_steps = (1:evolution_steps)';

valid_range = find(valid_steps & time_steps > 0 & time_steps < evolution_steps/2);
if length(valid_range) > 5

p = polyfit(time_steps(valid_range), divergence_curve(valid_range), 1);

lyapunov_exp = p(1);

else

lyapunov_exp = 0;

end
fprintf('最大Lyapunov指数: %.6f\n', lyapunov_exp);

end

6. 其他混沌指标

function hurst_exp = calculate_hurst_exponent(time_series)
    % 计算Hurst指数 - 重标极差法(R/S方法)
n = length(time_series);

max_window = floor(n/10);

window_sizes = unique(round(logspace(log10(10), log10(max_window), 20)));
rs_values = zeros(length(window_sizes), 1);
for i = 1:length(window_sizes)

window_size = window_sizes(i);

num_windows = floor(n / window_size);
window_rs = zeros(num_windows, 1);
for j = 1:num_windows

start_idx = (j-1)window_size + 1;

end_idx = jwindow_size;

window_data = time_series(start_idx:end_idx);
% 计算累积偏差

mean_val = mean(window_data);

deviations = window_data - mean_val;

cumulative_deviations = cumsum(deviations);
% 范围

R = max(cumulative_deviations) - min(cumulative_deviations);

S = std(window_data);
if S > 0

window_rs(j) = R / S;

else

window_rs(j) = 0;

end

end
rs_values(i) = mean(window_rs(window_rs > 0));

end
% 线性拟合得到Hurst指数

valid_idx = rs_values > 0 & window_sizes' > 0;

if sum(valid_idx) > 3

p = polyfit(log(window_sizes(valid_idx)), log(rs_values(valid_idx)), 1);

hurst_exp = p(1);

else

hurst_exp = 0.5;

end
fprintf('Hurst指数: %.4f\n', hurst_exp);

end
function [approx_entropy, sample_entropy] = calculate_entropy_measures(time_series)

% 计算近似熵和样本熵
m = 2; % 嵌入维数

r = 0.2 * std(time_series); % 容忍度
% 近似熵

approx_entropy = approximate_entropy(time_series, m, r);
% 样本熵

sample_entropy = sample_entropy_calc(time_series, m, r);
fprintf('近似熵: %.4f, 样本熵: %.4f\n', approx_entropy, sample_entropy);

end
function apen = approximate_entropy(data, m, r)

% 计算近似熵

n = length(data);
% 重构相空间

embedded_m = zeros(n-m+1, m);

for i = 1:n-m+1

embedded_m(i, 😃 = data(i:i+m-1);

end
% 计算距离矩阵

phi_m = calculate_phi(embedded_m, n-m+1, r);
% 对于m+1

embedded_m1 = zeros(n-m, m+1);

for i = 1:n-m

embedded_m1(i, 😃 = data(i:i+m);

end
phi_m1 = calculate_phi(embedded_m1, n-m, r);
% 计算近似熵

apen = log(phi_m / phi_m1);

end
function phi = calculate_phi(embedded_data, n, r)

% 计算Phi函数

count = zeros(n, 1);
for i = 1:n

for j = 1:n

if i ~= j

if max(abs(embedded_data(i, 😃 - embedded_data(j, 😃)) <= r

count(i) = count(i) + 1;

end

end

end

end
C = count / (n-1);

phi = mean(C);

end
function sampen = sample_entropy_calc(data, m, r)

% 计算样本熵

n = length(data);
% 重构相空间 - m维

embedded_m = zeros(n-m, m);

for i = 1:n-m

embedded_m(i, 😃 = data(i:i+m-1);

end
% 计算匹配数 - m维

A = 0;

for i = 1:n-m

for j = i+1:n-m

if max(abs(embedded_m(i, 😃 - embedded_m(j, 😃)) <= r

A = A + 1;

end

end

end
% 重构相空间 - m+1维

embedded_m1 = zeros(n-m-1, m+1);

for i = 1:n-m-1

embedded_m1(i, 😃 = data(i:i+m);

end
% 计算匹配数 - m+1维

B = 0;

for i = 1:n-m-1

for j = i+1:n-m-1

if max(abs(embedded_m1(i, 😃 - embedded_m1(j, 😃)) <= r

B = B + 1;

end

end

end
% 计算样本熵

if B > 0 && A > 0

sampen = -log(B / A);

else

sampen = 0;

end

end

7. 可视化函数

function plot_chaos_analysis_results(results, time_series)
    % 绘制混沌分析结果
figure('Position', [100, 100, 1400, 1000]);
% 原始时间序列

subplot(3, 4, 1);

plot(time_series, 'b-', 'LineWidth', 1);

title('原始时间序列');

xlabel('时间');

ylabel('幅值');

grid on;
% 自相关函数

subplot(3, 4, 2);

plot(1:length(results.autocorr_func), results.autocorr_func, 'ro-', 'LineWidth', 1.5);

hold on;

plot([results.tau_autocorr, results.tau_autocorr], [0, 1], 'r--', 'LineWidth', 2);

title('自相关函数法');

xlabel('延迟');

ylabel('自相关');

legend('自相关', '最优延迟', 'Location', 'best');

grid on;
% 互信息

subplot(3, 4, 3);

plot(1:length(results.mutual_info), results.mutual_info, 'go-', 'LineWidth', 1.5);

hold on;

plot([results.tau_mutual, results.tau_mutual], [0, max(results.mutual_info)], 'g--', 'LineWidth', 2);

title('互信息法');

xlabel('延迟');

ylabel('互信息');

legend('互信息', '最优延迟', 'Location', 'best');

grid on;
% 伪最近邻

subplot(3, 4, 4);

plot(1:length(results.fnn_ratio), results.fnn_ratio, 'bo-', 'LineWidth', 1.5);

hold on;

plot([results.embedding_dim, results.embedding_dim], [0, max(results.fnn_ratio)], 'b--', 'LineWidth', 2);

title('伪最近邻法');

xlabel('嵌入维数');

ylabel('伪最近邻比例');

legend('FNN比例', '最优维数', 'Location', 'best');

grid on;
% 关联维

subplot(3, 4, 5);

loglog(results.correlation_integral.radius, results.correlation_integral.C, 'mo-', 'LineWidth', 1.5);

title('关联积分');

xlabel('距离 r');

ylabel('关联积分 C(r)');

grid on;
% 双对数图

subplot(3, 4, 6);

plot(log10(results.correlation_integral.radius), log10(results.correlation_integral.C), 'mo-', 'LineWidth', 1.5);

title('关联维计算');

xlabel('log(r)');

ylabel('log(C(r))');

grid on;
% Lyapunov指数

subplot(3, 4, 7);

plot(1:length(results.lyapunov_curve), results.lyapunov_curve, 'co-', 'LineWidth', 1.5);

title('Lyapunov指数计算');

xlabel('演化步数');

ylabel('平均发散率');

grid on;
% 相空间重构 (3D)

if results.embedding_dim >= 3

subplot(3, 4, 8);

embedded_3d = phase_space_reconstruction(time_series, 3, results.tau_mutual);

plot3(embedded_3d(:,1), embedded_3d(:,2), embedded_3d(:,3), 'b-', 'LineWidth', 0.5);

title('相空间重构 (3D)');

grid on;

end
% 结果汇总

subplot(3, 4, [9, 10, 11, 12]);

text(0.1, 0.9, sprintf('混沌参数分析结果汇总:'), 'FontSize', 12, 'FontWeight', 'bold');

text(0.1, 0.7, sprintf('时间延迟 (自相关法): %d', results.tau_autocorr), 'FontSize', 10);

text(0.1, 0.6, sprintf('时间延迟 (互信息法): %d', results.tau_mutual), 'FontSize', 10);

text(0.1, 0.5, sprintf('嵌入维数: %d', results.embedding_dim), 'FontSize', 10);

text(0.1, 0.4, sprintf('关联维数: %.4f', results.correlation_dim), 'FontSize', 10);

text(0.1, 0.3, sprintf('最大Lyapunov指数: %.6f', results.lyapunov_exponent), 'FontSize', 10);

text(0.1, 0.2, sprintf('Hurst指数: %.4f', results.hurst_exponent), 'FontSize', 10);

text(0.1, 0.1, sprintf('近似熵: %.4f, 样本熵: %.4f', results.approximate_entropy, results.sample_entropy), 'FontSize', 10);
% 混沌特性判断

if results.lyapunov_exponent > 0

chaos_status = '具有混沌特性';

else

chaos_status = '无显著混沌特性';

end

text(0.5, 0.9, sprintf('系统状态: %s', chaos_status), 'FontSize', 12, 'Color', 'red', 'FontWeight', 'bold');
axis off;

end
function generate_analysis_report(results)

% 生成分析报告
fprintf('\n=== 混沌时间序列分析报告 ===\n');

fprintf('1. 时间延迟参数:\n');

fprintf('   - 自相关函数法: τ = %d\n', results.tau_autocorr);

fprintf('   - 互信息法: τ = %d\n', results.tau_mutual);
fprintf('2. 嵌入参数:\n');

fprintf('   - 最优嵌入维数: m = %d\n', results.embedding_dim);
fprintf('3. 分形特性:\n');

fprintf('   - 关联维数: D₂ = %.4f\n', results.correlation_dim);

fprintf('   - Hurst指数: H = %.4f\n', results.hurst_exponent);
fprintf('4. 动力学特性:\n');

fprintf('   - 最大Lyapunov指数: λ₁ = %.6f\n', results.lyapunov_exponent);

fprintf('   - 近似熵: %.4f\n', results.approximate_entropy);

fprintf('   - 样本熵: %.4f\n', results.sample_entropy);
% 混沌判断

fprintf('5. 混沌特性判断:\n');

if results.lyapunov_exponent > 0

fprintf('   - 系统具有混沌特性 (λ₁ > 0)\n');

if results.correlation_dim > 0 && ~isinf(results.correlation_dim)

fprintf('   - 系统具有分形吸引子\n');

end

else

fprintf('   - 系统无显著混沌特性\n');

end
fprintf('================================\n\n');

end

8. 使用示例和测试

function chaos_analysis_demo()
    % 混沌时间序列分析演示
fprintf('混沌时间序列参数计算工具箱演示\n');
% 生成测试数据

fprintf('生成测试时间序列...\n');
% 1. Lorenz系统

lorenz_data = generate_lorenz_system(1000);
% 2. 混沌逻辑映射

logistic_data = generate_logistic_map(0.4, 3.9, 1000);
% 3. 随机序列

random_data = randn(1000, 1);
% 分析Lorenz系统

fprintf('\n分析Lorenz系统:\n');

results_lorenz = chaotic_time_series_analysis(lorenz_data(:,1), 'max_delay', 30, 'max_dim', 8);
% 分析逻辑映射

fprintf('\n分析混沌逻辑映射:\n');

results_logistic = chaotic_time_series_analysis(logistic_data, 'max_delay', 20, 'max_dim', 6);
% 分析随机序列

fprintf('\n分析随机序列:\n');

results_random = chaotic_time_series_analysis(random_data, 'max_delay', 20, 'max_dim', 6);
% 比较结果

compare_results(results_lorenz, results_logistic, results_random);

end
function lorenz_data = generate_lorenz_system(n_points)

% 生成Lorenz系统数据

sigma = 10;

rho = 28;

beta = 8/3;
dt = 0.01;

x = zeros(n_points, 1);

y = zeros(n_points, 1);

z = zeros(n_points, 1);
% 初始条件

x(1) = 1; y(1) = 1; z(1) = 1;
% 数值积分

for i = 1:n_points-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
lorenz_data = [x, y, z];

end
function logistic_data = generate_logistic_map(x0, r, n_points)

% 生成逻辑映射数据

logistic_data = zeros(n_points, 1);

logistic_data(1) = x0;
for i = 2:n_points

logistic_data(i) = r * logistic_data(i-1) * (1 - logistic_data(i-1));

end

end
function compare_results(results1, results2, results3)

% 比较不同系统的分析结果
fprintf('\n=== 不同系统混沌特性比较 ===\n');

fprintf('系统\t\tLyapunov指数\t关联维数\tHurst指数\t混沌特性\n');

fprintf('----------------------------------------------------------------\n');
systems = {'Lorenz系统', '逻辑映射', '随机序列'};

results = {results1, results2, results3};
for i = 1:3

lyap = results{i}.lyapunov_exponent;

corr_dim = results{i}.correlation_dim;

hurst = results{i}.hurst_exponent;
if lyap > 0

chaos_char = '是';

else

chaos_char = '否';

end
fprintf('%-12s\t%.6f\t\t%.4f\t\t%.4f\t\t%s\n', ...

systems{i}, lyap, corr_dim, hurst, chaos_char);

end

end

参考工具箱 含有各种混沌时间序列参数（时间延迟，嵌入维数，关联维等）计算程序 www.3dddown.com/cna/82714.html

工具箱特性总结

功能模块	计算方法	输出参数
时间延迟	自相关函数法、互信息法	最优延迟τ
嵌入维数	伪最近邻法	最优嵌入维数m
关联维数	G-P算法	关联维D₂
Lyapunov指数	小数据量法	最大Lyapunov指数λ₁
其他指标	R/S法、熵计算	Hurst指数、近似熵、样本熵