计算李雅普诺夫特性指数的matlab函数

function lambda_1=lyapunov（dataNmtauP）
%  该函数用来计算时间序列的最大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（dataNmtau）;                    %相空间重构
    M=N-（m-1）*tau;                                     %重构相空间中相点的个数
    for i = 1 : （M-1）
        for j = i+1 : M
            d = 0;
            for k = 1 : m
                d = d + （Y（ki）-Y（kj））*（Y（ki）-Y（kj））;
            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（ki）-Y（k1））*（Y（ki）-Y（k1））;
        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（ki）-Y（kLoc_DK+1））*（Y（ki）-Y（kLoc_DK+1））;
        end
        DK1 = sqrt（DK1）;
        old_Loc_DK = Loc_DK ;                  % 保存原最近位置相点
        Loc_DK=DK;

%     计算前i个log2（DK1/DK）的累计和以及保存i点的李氏指数
        if （DK1 ~= 0）&&（ DK ~= 0）
           sum_lmd = sum_lmd + log（DK1/DK） /log（2）;
        end
        lm