function lambda_1=largest_lyapunov_exponent(data,fs)
%the function is used to calcultate largest lyapunov exponent with the mended algorithm,which put forward by lv jing hu.
%data:the time series
%N:the length of data
%m:enbedding dimention
%tau:time delay
%P:the mean period of the time series,calculated with FFT
%lambda_1:return the largest lyapunov exponent
%skyhawk
    N = length(data); 
    p = priod(data,fs);%時間序列平均周期
    t = 1;%7;%tau(data'); % 時延
    m = 1;%12;  % 分析GEZ信號時m=12，正弦信號 m =2;            
    % m = cao(data,t); %  嵌入維    
    delt_t=1/fs;
    Y=reconstitution(data,N,m,t);%reconstitute state space
    M=N-(m-1)*t;%M is the number of embedded points in m-dimensional space
for j=1:M
    d_max=1e+100;  
    d_min=d_max;
    for jj=1:M                                              %尋找相空間中每個點的最近距離點，并記下該點下標
        d_s=0;                                              %
        if abs(j-jj) > p                                     %限制短暫分離,P為時間序列平均周期
            for i=1:m
                d_s=d_s+(Y(i,j)-Y(i,jj))*(Y(i,j)-Y(i,jj)); 
            end  %for i=1:m
                if d_s^0.5 < d_min
                   d_min=d_s^0.5;
                   idx_j=jj;
                end  %if d_s^0.5 < d_min
          
        end  %if abs(j-jj)>P
    end  %for jj=1:M   

    max_i=min((M-j),(M-idx_j));%計算點j的最大演化時間步長i
    for k=1:max_i              %計算點j與其最近鄰點在i個離散步長后的距離
        d_j_i=0;
        for kk=1:m
            d_j_i=d_j_i+(Y(kk,j+k)-Y(kk,idx_j+k))*(Y(kk,j+k)-Y(kk,idx_j+k)); 
        end  % for kk=1:m
        d(k,j)=d_j_i^0.5;  
    end  %  for k=1:max_i 
end  %for j=1:M

%對每個演化時間步長i，求所有的j的lnd(i,j)平均
[l_i,l_j]=size(d);
for i=1:l_i
    q=0;
    y_s=0;
    for j=1:l_j
        if d(i,j)~=0
            q=q+1;
            y_s=y_s+log(d(i,j));
        end  % if d(i,j)~=0
    end  % for j=1:l_j
    y(i)=y_s/(q*delt_t);
end  %for i=1:l_i
x=1:length(y);
pp=polyfit(x,y,1);
lambda_1=pp(1);
% yp=polyval(pp,x);
%plot(x,y,'-o',x,yp,'--')