求助：小数据量方法求解le指数
本帖最后由 竹园 于 2012-12-20 11:56 编辑小数据量方法计算le指数，数据采用lorenz序列10001-13000样本区间，cc法算出m=21，tau=10，le指数算出来是0.0122，但是在使用这个程序算其他的数据却得不到结果，望路过的大侠指教，我是刚刚动手学编程，程序参考了yuling的贴，自己作了修改，data0是我的数据，就是算不出le指数，在运行的时候总是出现没定义
主程序
clear all
t0=0;
tf=130;
=ode45(@Lorenz,,[-1,0,1]);
Lorenz_data=x(10002:end,1);
m=21;
tau=10;
tic
lambda_1=largest_lyapunov_exponent(Lorenz_data,m,tau,t);
toc
function dy=Lorenz(t,y)
%方程的定义
a=16;
b=4.0;
c=45.92;
dy=zeros(3,1);
dy(1)=-a*(y(1)-y(2));
dy(2)=-y(1)*y(3)+c*y(1)-y(2);
dy(3)=y(1)*y(2)-b*y(3);
function T_mean=period_mean_fft(data)
%该函数使用快速傅里叶变换FFT计算序列平均周期
%data：时间序列
%T_mean：返回快速傅里叶变换FFT计算出的序列平均周期
Y = fft(data); %快速FFT变换
N = length(data); %FFT变换后数据长度
Y(1) = []; %去掉Y的第一个数据，它是data所有数据的和
power = abs(Y(1:N/2)).^2;%求功率谱
nyquist = 1/2;
freq = (1:N/2)/(N/2)*nyquist; %求频率
subplot(121)
plot(freq,power); grid on ; %绘制功率谱图
xlabel('频率');
ylabel('功率');
title('功率谱图');
period = 1./freq; %计算周期
subplot(122)
plot(period,power); grid on%绘制周期－功率谱曲线
ylabel('功率');
xlabel('周期');
title('周期—功率谱图');
= max(power); %求最高谱线所对应的下标
T_mean=period(index) %由下标求出平均周期
function lambda_1=largest_lyapunov_exponent(data,m,tau,t)
%本函数使用小数据量方法计算最大lyapunov指数
%data:时间序列
%m:嵌入维
%tau:时间延迟
%P:使用 FFT计算出的时间序列平均周期
%t:抽样时间间隔
%lambda_1:函数返回的最大lyapunov指数值
N=length(data);%序列长度
P=period_mean_fft(data);
P=fix(P)%时间序列的平均周期
delt_t=1;
Y=reconstitution(data,m,tau );%重构相空间
M=N-(m-1)*tau;%M是m维重构相空间中总点数
d=zeros(M-1,M);
for j=1:M
d_min=1e+100;
for jj=1:M %寻找相空间中每个点的最近距离点，并记下该点下标
if abs(j-jj)>P %限制短暂分离
d_s=norm(Y(:,j)-Y(:,jj));%计算分离后两点的距离
if d_s<d_min
d_min=d_s;
idx_j=jj;
end
end
end
max_i=min((M-j),(M-idx_j)); %计算点j的最大演化时间步长i
for k=1:max_i %计算点j与其最近邻点在i个离散步后的距离
d(k,j)=norm(Y(:,j+k)-Y(:,idx_j+k));
end
end
%对每个演化时间步长i，求所有的j的lnd(i,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
end
if q>0
y(i)=y_s/(q*delt_t);
end
end
I=1:length(y);
pp=polyfit(I,y,1);
lambda_1=pp(1);
yp=polyval(pp,I);
figure(1)
plot(I,y,'-');
xlabel('i');
ylabel('y(i)');
hold on;
for i=1:length(y)-1
x(i)=y(i+1)-y(i);
end
figure(2)
plot(I(1:(length(y)-1)),x);grid;
xlabel('i');
ylabel('y(i)-y(i-1)');
hold on;
linear=input('请输入线型部分的长度')
I1=1:linear;
PP=polyfit(I1,y(I1),1);
lambda_1=PP(1)
YP=polyval(PP,I1);
plot(I1,YP,'r--');
页:
[1]