www.pudn.com > C-C.rar > lyapunov_rosenstein.m, change:2004-03-23,size:2407b

```function LCE1=lyapunov_rosenstein(x,m,t,P,p,ts)
% Syntax:LCE1=lyapunov_rosenstein(x,m,t,P,p)
% Calculate maximum lyapunov exponent from small sets
% "A practical method for calculating largest lyapunov
%  exponents from small sets"--Michael T.Rosenstein
%Input:
%     x----time series coloum format
%     m----embedding dimension
%     t----time delay
%     P----mean perieod time
%     p----p norm
%     ts---Sampling frequence
%Output:
%     LCE1-Dominant lyapunov exponents
%
%Usage:
%    x=[];!!!!Remember: in coloum format
%    m=30;
%    t=10;
%    P=20;%gained by FFT
%    p=inf;
%    ts=0.001;
%    LCE1=lyapunov_rosenstein(x,m,t,P,p,ts);
%Related routine:
%        phasespace.m
%Author:skyhawk
% Created on 2004.3.8
% Modified on 2004.3.23
% Air Force Engineering University
% Air Force Engineering Institute
% Dept.1, Shaan Xi, Xi'an 710038, PR China.
% Email:xunkai_wei@163.com
%
[Y,M]=phasespace(x,m,t);
%calculate two points' distance in phase space and determine the nearest neighbor
%count computation times
for i=1:M
%Calculate original distance
count_refer=1;
Refer=Y(i,:);
for j=1:M
if abs(j-i)<P     %constrain temporal seperation
continue;
end
dist(count_refer)=norm(Refer-Y(j,:),p);% record its value
index_count_refer(count_refer)=j;% remember its index
count_refer=count_refer+1;%modify counter
end
%sort distance
[dist,h]=sort(dist);
%get refer's nearest neighbor
Dist0(i)=dist(1);
Index_Neighbor(i)=index_count_refer(h(1));
end
%对于每个相点及其领域点离散k步后，计算当前的相空间距离
%
for j=1:M
orignal_point=Y(j,:)
neighbor_point=Y(Index_Neighbor(j),:);
for i=1:min(M-j,M-Index_Neighbor(j))%i discrete times
orignal_point_i=Y(j+i,:)
neighbor_point_i=Y(Index_Neighbor(j)+i,:);
dist1_k(j,i)=norm(orignal_point_i-neighbor_point_i,p);% distance accordingly
end
end
%work out <ln(dj(i))>
[row,col]=size(dist1_k);
for i=1:col
cout_j_for_i=1;
sum_j_for_i=0;
for j=1:row
if dist1_k(j,i)~=0
sum_j_for_i=sum_j_for_i+log(dist1_k(j,i));
cout_j_for_i=cout_j_for_i+1;
end
end
aver_j_for_i(i)=sum_j_for_i/cout_j_for_i;
end
X=1:col;
Y=aver_j_for_i/ts;
FIT=polyfit(X,Y,1);
%plot(X,Y)%,X,polyval(FIT,X));
LCE1=FIT(1);```