|
|
补充
我用的小数量算法程序如下：
function lambda_1=largest_lyapunov_p3(data,N,m,tau,P,delt_t)
m=8;tau=9;
delt_t=1;
%**********************************************************
x=data;
xPower=abs(fft(x)).^2;
NN=length(xPower);
xPower(1)=[];%去除直流分量
NN=floor(NN/2);
xPower=xPower(1:NN);
freq=(1:NN)/NN*0.5;
[mP,index]=max(xPower);
P=index;
%*************************************************************
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 %限制短暂分离
d_s=sum(abs(Y(j)-Y(jj)));
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;
d_j_i=norm(Y(j+k)-Y(idx_j+k));
d_content=[d_content,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
XX=1:length(y);
figure;
plot(XX,y(XX),'-','markersize',1);hold on;%——————————【上面的图是从这儿画出来的】
linear=input('请输入线形部分的长度');
XX1=1:linear;
pp=polyfit(XX1,y(XX1),1);
lambda_1=pp(1)
yp=polyval(pp,XX1);
figure;
plot(XX1,yp,'r--');
function X=reconstitution(data,N,m,tau)
%该函数用来重构相空间
% m为嵌入空间维数
% tau为时间延迟
% data为输入时间序列
% N为时间序列长度
% X为输出,是m*n维矩阵
M=N-(m-1)*tau;%相空间中点的个数
for j=1:M %相空间重构
for i=1:m
X(i,j)=data((i-1)*tau+j);
end
end |
|
发表于 2008-11-7 16:10
显身卡
