function lambda_1=lyapunov_wolf(data,N,m,tau,P)
%  該函數用來計算時間序列的最大Lyapunov 指數--Wolf 方法
%  m: 嵌入維數
%  tau:時間延遲
%  data:時間序列
%  N:時間序列長度
%  P:時間序列的平均周期,選擇演化相點距當前點的位置差，即若當前相點為I，則演化相點只能在|I－J|>P的相點中搜尋
%  lambda_1:返回最大lyapunov指數值
min_point=1  ; %&&要求最少搜索到的點數
MAX_CISHU=5 ;  %&&最大增加搜索范圍次數
%FLYINGHAWK
%   求最大、最小和平均相點距離
    max_d = 0;                                         %最大相點距離
    min_d = 1.0e+100;                                  %最小相點距離
    avg_dd = 0;
    Y=reconstitution(data,N,m,tau);                    %相空間重構
    M=N-(m-1)*tau;                                     %重構相空間中相點的個數
    for i = 1 : (M-1)
        for j = i+1 : M
            d = 0;
            for k = 1 : m
                d = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
            end
            d = sqrt(d);
            if max_d < d
               max_d = d;
            end
            if min_d > d
               min_d = d;
            end
            avg_dd = avg_dd + d;
        end
    end
    avg_d = 2*avg_dd/(M*(M-1));                %平均相點距離
    
    dlt_eps = (avg_d - min_d) * 0.02 ;         %若在min_eps～max_eps中找不到演化相點時，對max_eps的放寬幅度
    min_eps = min_d + dlt_eps / 2 ;            %演化相點與當前相點距離的最小限
    max_eps = min_d + 2 * dlt_eps  ;           %&&演化相點與當前相點距離的最大限
    
%     從P+1～M-1個相點中找與第一個相點最近的相點位置(Loc_DK)及其最短距離DK
    DK = 1.0e+100;                             %第i個相點到其最近距離點的距離
    Loc_DK = 2;                                %第i個相點對應的最近距離點的下標
    for i = (P+1):(M-1)                        %限制短暫分離，從點P+1開始搜索
        d = 0;
        for k = 1 : m
            d = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));
        end
        d = sqrt(d);
        if (d < DK) & (d > min_eps) 
           DK = d;
           Loc_DK = i;
        end
    end
%     以下計算各相點對應的李氏數保存到lmd()數組中
%     i 為相點序號，從1到(M-1)，也是i-1點的演化點；Loc_DK為相點i-1對應最短距離的相點位置，DK為其對應的最短距離
%     Loc_DK+1為Loc_DK的演化點，DK1為i點到Loc_DK+1點的距離，稱為演化距離
%     前i個log2（DK1/DK）的累計和用于求i點的lambda值
    sum_lmd = 0 ;                              % 存放前i個log2（DK1/DK）的累計和
    for i = 2 : (M-1)                          % 計算演化距離      
        DK1 = 0;
        for k = 1 : m
            DK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));
        end
        DK1 = sqrt(DK1);
        old_Loc_DK = Loc_DK ;                  % 保存原最近位置相點
        old_DK=DK;

%     計算前i個log2（DK1/DK）的累計和以及保存i點的李氏指數
        if (DK1 ~= 0)&( DK ~= 0)
           sum_lmd = sum_lmd + log(DK1/DK) /log(2);
        end
        lmd(i-1) = sum_lmd/(i-1);
%     以下尋找i點的最短距離：要求距離在指定距離范圍內盡量短，與DK1的角度最小
        point_num = 0  ; % &&在指定距離范圍內找到的候選相點的個數
        cos_sita = 0  ; %&&夾角余弦的比較初值 ——要求一定是銳角
        zjfwcs=0     ;%&&增加范圍次數
         while (point_num == 0)
           % * 搜索相點
            for j = 1 : (M-1)
                if abs(j-i) <=(P-1)      %&&候選點距當前點太近，跳過！
                   continue;     
                end
                
                %*計算候選點與當前點的距離
                dnew = 0;
                for k = 1 : m
                   dnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
                end
                dnew = sqrt(dnew);
                
                if (dnew < min_eps)|( dnew > max_eps )   %&&不在距離范圍，跳過！
                  continue;             
                end
                               
                %*計算夾角余弦及比較
                DOT = 0;
                for k = 1 : m
                    DOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));
                end
                CTH = DOT/(dnew*DK1);
                
                if acos(CTH) > (3.14151926/4)      %&&不是小于45度的角，跳過！
                  continue;
                end
                
                if CTH > cos_sita   %&&新夾角小于過去已找到的相點的夾角，保留
                    cos_sita = CTH;
                    Loc_DK = j;
                    DK = dnew;
                end

                point_num = point_num +1;
                
            end        
        
            if point_num <= min_point
               max_eps = max_eps + dlt_eps;
               zjfwcs =zjfwcs +1;
               if zjfwcs > MAX_CISHU    %&&超過最大放寬次數，改找最近的點
                   DK = 1.0e+100;
                   for ii = 1 : (M-1)
                      if abs(i-ii) <= (P-1)      %&&候選點距當前點太近，跳過！
                       continue;     
                      end
                      d = 0;
                      for k = 1 : m
                          d = d + (Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii));
                      end
                      d = sqrt(d);
        
                      if (d < DK) & (d > min_eps) 
                         DK = d;
                         Loc_DK = ii;
                      end
                   end
                   break; 
               end
               point_num = 0          ;     %&&擴大距離范圍后重新搜索
               cos_sita = 0;
            end
        end
   end

%取平均得到最大李雅普諾夫指數
lambda_1=sum(lmd)/length(lmd);