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