竹园 发表于 2012-12-20 11:50

求助:小数据量方法求解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]
查看完整版本: 求助:小数据量方法求解le指数