|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
小数据量方法计算le指数，数据采用lorenz序列10001-13000样本区间，cc法算出m=21，tau=10，le指数算出来是0.0122，但是在使用这个程序算其他的数据却得不到结果，望路过的大侠指教，我是刚刚动手学编程，程序参考了yuling的贴，自己作了修改，data0是我的数据，就是算不出le指数，在运行的时候总是出现没定义
主程序
clear all
t0=0;
tf=130;
[t,x]=ode45(@Lorenz,[t0:0.01:tf],[-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('周期—功率谱图');
[mp,index] = 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)平均
[l_i,l_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--');
data0.txt
(40.32 KB, 下载次数: 1)
Lozenz
|
|
发表于 2012-12-20 18:10
提升卡
变色卡
显身卡