largest_lyapunov_exponent.m

检测时间序列的平稳性

function lambda_1=largest_lyapunov_exponent(data,N,m,tau,P,delt_t)
d_length=[];
d_content=[];
Y=reconstitution(data,N,m,tau);
M=N-(m-1)*tau;
idx_j=0;
for j=1:M
    d_min=10000; 
    for jj=1:M                                              
        d_s=0;      %寻找相空间中每个点的最近距离点，并记下该点下标
        if abs(j-jj)>P                                      %限制短暂分离
            for i=1:m
                d_s=d_s+abs(Y(j,i)-Y(jj,i));
            end
            if d_s<d_min
               d_min=d_s;
               idx_j=jj;
            end
        end
    end
    if ((M-j)>(M-idx_j));%计算点j的最大演化时间步长i
        max_i=M-idx_j;
    else
        max_i=M-j;
    end
    d_length=[d_length,max_i];
    for k=1:max_i              %计算点j与其最近邻点在i个离散步后的距离
        d_j_i=0;
        for kk=1:m
            d_j_i=d_j_i+(Y(j+k,kk)-Y(idx_j+k,kk))^2;
        end
        d_content=[d_content,sqrt(d_j_i)];
    end
end

%对每个演化时间步长i，求所有的j的lnd(i,j)平均
y=[];
for i=1:max(d_length)
   S_j_i=0;
   sum_former=0;
   Count=0;
   for j=1:M
       if j==1
          former=0;
       else
          former=d_length(j-1);
       end    
       sum_former=sum_former+former;
       if i<=(d_length(j))
           if d_content(sum_former+i)~=0
              S_j_i=S_j_i+log(d_content(sum_former+i));
              Count=Count+1;
           end
       end
   end
   y=[y,S_j_i/(Count*delt_t)]; %对每个演化时间步长i，求所有的j平均
end
x=1:length(y);
%x=1:10;
pp=polyfit(x,y(x),1);
lambda_1=pp(1);
yp=polyval(pp,x);
plot(x,y(x),'-o',x,yp,'--');