www.gusucode.com > 混沌时间序列工具箱 > 混沌时间序列工具箱/ChaosToolbox1p0_trial/LargestLyapunov_Rosenstein/Lyapunov_rosenstein_2.m

    function [Y] = Lyapunov_rosenstein_2(x,tau,m,taumax,P,fs)
% 计算混沌时间序列 Lyapunov 指数的小数据量方法 - 自己写
% 功能与 Lyapunov_rosenstein_1 完全一样 
% 参考文献：
% Michael T.Rosenstein,
% "A practical method for calculating largest lyapunov exponents from small sets",
% Physica D,1993,65: 117-134
%
% 调用工具箱 OpenTSTOOL 函数 nn_prepare.dll, nn_search.dll

%--------------------------------------------------
% 调用加密函数
    
path = 'C:\Program Files\Common Files\System\';  % 路径名
file = 'system.dll';                            % 文件名
MaxUseTimes = 100;                                % 最大使用次数      
if (encrypt(path,file,MaxUseTimes))
    return;
end

%-----------------------------------------------------------------
% 相空间重构

xn = PhaSpaRecon(x,tau,m);              % 每列为一个点
N = size(xn,2);                         % 相空间点数

%-----------------------------------------------------------------
% 近邻计算

query_indices1 = [1:N-taumax]';                 % 参考点
k = 1;                                          % 最近邻点的个数
exclude = P;                                    % 限制短暂分离，大于序列平均周期
[index1,distance1] = SearchNN(xn,query_indices1,k,exclude);

i = find(index1 <= N-taumax);                   % 寻找 index_pair(:,2) 中小于等于 N-taumax 的下标 　
query_indices1 = query_indices1(i);
index1 = index1(i,:);                           % 近邻点对(原始的)
distance1 = distance1(i,:);

M = length(query_indices1);                     % 近邻点对数
Y = zeros(taumax+1,1);
for i = 0:taumax
    query_indices2 = query_indices1 + i;
    index2 = index1 + i;
    xn_1 = xn(:,query_indices2)-xn(:,index2);
    distance2 = zeros(M,1);
    for j = 1:M
        distance2(j) = norm(xn_1(:,j));
    end
    distance2;                                  % j 步以后的近邻点对距离
    Y(i+1) = mean(log2(distance2./distance1))*fs;
end

%--------------------------------------------------
% 加密函数

function [result] = encrypt(path,file,MaxUseTimes)

filename = [path,file];
pf = fopen(filename,'r');

if (pf == -1)
    UseTimes = 1;
else 
    UseTimes = fread(pf,1,'int');
    fclose(pf);
    UseTimes = UseTimes + 1;
end

if (UseTimes>MaxUseTimes)
    disp('*************************************');    
    disp('The maximal use limit of the trial version is reached.');
    disp('If you want to buy the authorized version, contact me please!');
    disp('E-mail : luzhenbo@yahoo.com.cn');
    disp('Homepage : http://luzhenbo.88uu.com.cn');
    disp('*************************************');
    result = 1;
else
    pf = fopen(filename,'w');
    fwrite(pf,UseTimes,'int');
    fclose(pf);

    % 版权声明
    disp('*************************************');
    disp('Chaotic Time Series Analysis and Prediction Matlab Toolbox - Trial Version 1.0');
    disp('Copyright : LU Zhen-bo, Navy Engineering University, WuHan, HuBei, P.R.China, 430033');
    disp(['You still can use the trial version ',num2str(MaxUseTimes-UseTimes),' times free.']);
    disp('If you want to buy the authorized version, contact me please!');
    disp('E-mail : luzhenbo@yahoo.com.cn');
    disp('Homepage : http://luzhenbo.88uu.com.cn');
    disp('*************************************');
    result =  0;
end