% 產生 Henon 序列　
%x(n+1) = 1 - a * x(n)^2 + y(n);  
%y(n+1) = b * x(n)
clc
clear
close all

a = 1.4;
b = 0.3;

x0 = 0;
y0 = 0;

N1 = 2000;                   % 前面的迭代點數
N2 = 1000;                   % 后面的迭代點數

z = zeros(N1+N2,2);
for i = 1:N1+N2
    x = 1 - a * x0^2 + y0 ;
    y = b * x0;
    x0 = x;
    y0 = y;
    
    z(i,1) = x;
    z(i,2) = y;
end
z = z(N1+1:end,:);

xn = z(:,1);
Y = z(:,2);

% henon 入口參數%
tau = 1;                        % 時延
m = 3;                          % 嵌入維
P = 100;                        % 序列平均周期
fs = 1;                         % 采樣頻率
evolvement_length=20;           % 最大離散演化步數
norm=2;

figure;
radius=0.0005;
for m=3:5
    data = PhaSpaRecon(xn,tau,m);    % 每列為一個點
    data = data';
    Y = lyapunov_kantz_c(data,tau,m,evolvement_length,norm,P,fs,radius);
    plot(Y,'g.-');
    hold on;
end

radius=0.002;
for m=3:5
    data = PhaSpaRecon(xn,tau,m);    % 每列為一個點
    data = data';
    Y = lyapunov_kantz_c(data,tau,m,evolvement_length,norm,P,fs,radius);
    plot(Y,'r.-');
    hold on;
end

radius=0.008;
for m=3:5
    data = PhaSpaRecon(xn,tau,m);    % 每列為一個點
    data = data';
    Y = lyapunov_kantz_c(data,tau,m,evolvement_length,norm,P,fs,radius);
    plot(Y,'b.-');
    hold on;
end
xlabel('時間延遲tau');
ylabel('距離增長S(tau)');