%function LE1=lya_wolf(data,tau,min_m,max_m,P)
%Wolf在一維數據中提取Lyapunov指數的方法
%(基于Lyapunov指數的觀測數據短期預測.陳繼光.山東大學 土建及水利學院)
min_m=1;max_m=20;
tau=1;P=2;
k=1;%取時間步長等于tau

N=length(data);

for m=min_m:max_m
    m
    Y=reconstitution(data,N,m,tau);
    M=N-(m-1)*tau;
    %-----------求距每個向量最短距離相點和距離-------------
    for j=1:M
        h=1;
        for jj=1:M         %尋找相空間中每個點的最近距離點，并記下該點下標限制短暫分離
            if abs(j-jj)>P
                dd(h,1)=jj;
                dd(h,2)=norm(Y(:,j)-Y(:,jj));
                h=h+1;
            end
        end
        min_d=min(dd(:,2));
        for s=1:length(dd(:,2))
            if min_d==dd(s,2)
                idx_j=dd(s,1);
                break
            end
        end
        D(j,1)=idx_j;%第一列存放序號
        D(j,2)=min_d;%第二列存放最短距離
    end

    %---------計算lya指數-----------------------
    lmd=0;
    n=0;%迭代次數初值
    for i=1:M
        if i+k<=M
            lmd=lmd+log2(D(i+k,2)/D(i,2))/k;
            n=n+1;
        end
    end
    LE1(m-min_m+1)=lmd/n;%取各指數增長率的平均值為最大Lyapunov指數估計
    clear dd D
end

%----------畫圖--------------
plot(min_m:max_m,LE1)