分形参数计算程序分享&求帮助

zxp2008 · 发表于 2006-12-5 10:17

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

最近有想法在课题里使用下面这几个参数,收集了几个程序,但不知道咋使用!
求大家帮帮忙,我的QQ419398633油箱是:zxp200866@sina.com

计算hurst指数

function [logRS,logERS,V]=RSana(x,n,method,q)
% 用 R/S 方法分析间序列
% logRS 是 log(R/S).
% logERS 是 log(R/S)期望.
% V 是统计量.
% x 是时间序列.
% n 是这个数列的子集.
% method 可以取下列值
% 'Hurst' 为了Hurst-Mandelbrot变量
% 'Lo' 是Lo变量.
% 'MW' 是Moody-Wu变量.
% 'Parzen' 是Parzen变量.
% q 可以是任意值
% a 是非0整数.
% 'auto' 是 Lo的默认值.
if nargin<1 | isempty(x)==1
error('你应该给出一个时间序列.');
else
% x 必须是变量
if min(size(x))>1
error('时间序列无效.');
end
x=x(:);
% N 是时间序列的长度
N=length(x);
end
if nargin<2 | isempty(n)==1
n=1;
else
% n 必须是一个变化的标量或矢量
if min(size(n))>1
error('n 必须是一个变化的标量或矢量.');
end
% n 必须是个整数
if n-round(n)~=0
error('n 必须是个整数.');
end
% n 必须是确定
if n<=0
error('n 必须是确定.');
end
end
if nargin<4 | isempty(q)==1
q=0;
else
if q=='auto'
t=autocorr(x,1);
t=t(2);
q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);
else
% q 必须是标量
if sum(size(q))>2
error('q 必须是标量.');
end
% q 必须是整数
if q-round(q)~=0
error('q 必须是整数.');
end
% q 必须是确定
if q<0
error('q 必须是确定.');
end
end
end
for i=1:length(n)
% 计算这个子序列
a=floor(N/n(i));
% 仓健这个子序列的矩阵
X=reshape(x(1:a*n(i)),n(i),a);
% 估算这个子序列的平均值
ave=mean(X);
% 给这个序列的每一个值除以平均值
cumdev=X-ones(n(i),1)*ave;
% 估算累计离差
cumdev=cumsum(cumdev);
% 估算这个标准偏差
switch method
case 'Hurst'
% Hurst-Mandelbrot 参数
stdev=std(X);
case 'Lo'
% Lo 参数
for j=1:a
sq=0;
for k=0:q
v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
if k>0
sq=sq+(1-k/(q+1))*v(k+1);
end
end
stdev(j)=sqrt(v(1)+2*sq);
end
case 'MW'
% Moody-Wu 参数
for j=1:a
sq1=0;
sq2=0;
for k=0:q
v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
if k>0
sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);
sq2=sq2+(1-k/(q+1))*v(k+1);
end
end
stdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);
end
case 'Parzen'
% Parzen 参数
if mod(q,2)~=0
error('在"Parzen" 参数中q 必须是2.');
end
for j=1:a
sq1=0;
sq2=0;
for k=0:q
v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
if k>0 & k<=q/2
sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);
elseif k>0 & k>q/2
sq2=sq2+(1-(k/q)^3)*v(k+1);
end
end
stdev(j)=sqrt(v(1)+2*sq1+2*sq2);
end
otherwise
error('你应该付给 "method"另一个值.');
end
% 估算(R(t)/s(t))
rs=(max(cumdev)-min(cumdev))./stdev;
clear stdev
% 取这个平均值 R/S的对数
logRS(i,1)=log10(mean(rs));
if nargout>1
% 开始计算log(E(R/S))
j=1:n(i)-1;
s=sqrt((n(i)-j)./j);
s=sum(s);
% 估算log(E(R/S))
logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
%其它估算log(E(R/S))
%logERS(i,1)=log10((n(i)-0.5)/n(i)*s/sqrt(n(i)*pi/2));
%logERS(i,1)=log10(sqrt(n(i)*pi/2));
end
if nargout>2
% 估算 V
V(i,1)=mean(rs)/sqrt(n(i));
end
end

复制代码

计算lyapunov

function lambda_1=lyapunov_wolf(data,N,m,tau,P)
% 该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法
% m: 嵌入维数
% tau:时间延迟
% data:时间序列
% N:时间序列长度
% P:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|I－J|>P的相点中搜寻
% lambda_1:返回最大lyapunov指数值
min_point=1 ; %&&要求最少搜索到的点数
MAX_CISHU=5 ; %&&最大增加搜索范围次数
%FLYINGHAWK
% 求最大、最小和平均相点距离
max_d = 0; %最大相点距离
min_d = 1.0e+100; %最小相点距离
avg_dd = 0;
Y=reconstitution(data,N,m,tau); %相空间重构
M=N-(m-1)*tau; %重构相空间中相点的个数
for i = 1 : (M-1)
for j = i+1 : M
d = 0;
for k = 1 : m
d = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
end
d = sqrt(d);
if max_d < d
max_d = d;
end
if min_d > d
min_d = d;
end
avg_dd = avg_dd + d;
end
end
avg_d = 2*avg_dd/(M*(M-1)); %平均相点距离
dlt_eps = (avg_d - min_d) * 0.02 ; %若在min_eps～max_eps中找不到演化相点时，对max_eps的放宽幅度
min_eps = min_d + dlt_eps / 2 ; %演化相点与当前相点距离的最小限
max_eps = min_d + 2 * dlt_eps ; %&&演化相点与当前相点距离的最大限
% 从P+1～M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离DK
DK = 1.0e+100; %第i个相点到其最近距离点的距离
Loc_DK = 2; %第i个相点对应的最近距离点的下标
for i = (P+1):(M-1) %限制短暂分离，从点P+1开始搜索
d = 0;
for k = 1 : m
d = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));
end
d = sqrt(d);
if (d < DK) & (d > min_eps)
DK = d;
Loc_DK = i;
end
end
% 以下计算各相点对应的李氏数保存到lmd()数组中
% i 为相点序号，从1到(M-1)，也是i-1点的演化点；Loc_DK为相点i-1对应最短距离的相点位置，DK为其对应的最短距离
% Loc_DK+1为Loc_DK的演化点，DK1为i点到Loc_DK+1点的距离，称为演化距离
% 前i个log2（DK1/DK）的累计和用于求i点的lambda值
sum_lmd = 0 ; % 存放前i个log2（DK1/DK）的累计和
for i = 2 : (M-1) % 计算演化距离
DK1 = 0;
for k = 1 : m
DK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));
end
DK1 = sqrt(DK1);
old_Loc_DK = Loc_DK ; % 保存原最近位置相点
old_DK=DK;
% 计算前i个log2（DK1/DK）的累计和以及保存i点的李氏指数
if (DK1 ~= 0)&( DK ~= 0)
sum_lmd = sum_lmd + log(DK1/DK) /log(2);
end
lmd(i-1) = sum_lmd/(i-1);
% 以下寻找i点的最短距离：要求距离在指定距离范围内尽量短，与DK1的角度最小
point_num = 0 ; % &&在指定距离范围内找到的候选相点的个数
cos_sita = 0 ; %&&夹角余弦的比较初值 ——要求一定是锐角
zjfwcs=0 ;%&&增加范围次数
while (point_num == 0)
% * 搜索相点
for j = 1 : (M-1)
if abs(j-i) <=(P-1) %&&候选点距当前点太近，跳过！
continue;
end
%*计算候选点与当前点的距离
dnew = 0;
for k = 1 : m
dnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
end
dnew = sqrt(dnew);
if (dnew < min_eps)|( dnew > max_eps ) %&&不在距离范围，跳过！
continue;
end
%*计算夹角余弦及比较
DOT = 0;
for k = 1 : m
DOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));
end
CTH = DOT/(dnew*DK1);
if acos(CTH) > (3.14151926/4) %&&不是小于45度的角，跳过！
continue;
end
if CTH > cos_sita %&&新夹角小于过去已找到的相点的夹角，保留
cos_sita = CTH;
Loc_DK = j;
DK = dnew;
end
point_num = point_num +1;
end
if point_num <= min_point
max_eps = max_eps + dlt_eps;
zjfwcs =zjfwcs +1;
if zjfwcs > MAX_CISHU %&&超过最大放宽次数，改找最近的点
DK = 1.0e+100;
for ii = 1 : (M-1)
if abs(i-ii) <= (P-1) %&&候选点距当前点太近，跳过！
continue;
end
d = 0;
for k = 1 : m
d = d + (Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii));
end
d = sqrt(d);
if (d < DK) & (d > min_eps)
DK = d;
Loc_DK = ii;
end
end
break;
end
point_num = 0 ; %&&扩大距离范围后重新搜索
cos_sita = 0;
end
end
end
%取平均得到最大李雅普诺夫指数
lambda_1=sum(lmd)/length(lmd);

复制代码

G-P算法计算关联维数

function [ln_r,ln_C]=G_P(data,N,tau,min_m,max_m,ss)
% 此程序是用G-P算法计算关联维数Dc
% N 是时间序列的长度
% tau 是固定时间间隔
% min_m最小的嵌入维数
% max_m最大的嵌入维数
% ssr的步长
for m=min_mmax_m
Y=reconstitution(data,N,m,tau);%重建矢量空间
M=N-(m-1)tau;%矢量空间的点数
for i=1M-1
for j=i+1M
d(i,j)=max(abs(Y(,i)-Y(,j)));%计算其余点到点Xi的距离
end
end
max_d=max(max(d));%是所有点中距离最远的点
d(1,1)=max_d;
min_d=min(min(d));%是所有点中距离最近的点
delt=(max_d-min_d)ss;%是r的步长
for k=1ss
r=min_d+kdelt;
C(k)=correlation_integral(Y,M,r);%计算关联积分函数
ln_C(m,k)=log(C(k));%lnC(r)
ln_r(m,k)=log(r);%lnr
fprintf('%d%d%d%dn',k,ss,m,max_m);
end
plot(ln_r(m,),ln_C(m,));
hold on;
end
fid=fopen('lnr.txt','w');
fprintf(fid,'%6.2f %6.2fn',ln_r);
fclose(fid);
fid = fopen('lnC.txt','w');
fprintf(fid,'%6.2f %6.2fn',ln_C);
fclose(fid);

复制代码

计算kolmogorov

clear all
x=load('are.dat')
n=length(x)
sd=std(x)
r=0.2*sd
for ii=1:2
m=ii+1;
num=zeros(n-m+1,1);
for i=1:n-m+1
for j=1:n-m+1
if j~=i
for k=1:m
d(k)=abs(x(i+k-1)-x(j+k-1));
end
d1=max(d);
if d1<r
num(i,1)=num(i,1)+1;
end
end
end
c0(i)=num(i,1)/(n-m);
c1(i)=log(c0(i));
end
sc=sum(c1);
fi(ii)=sc/(n-m+1);
end
app=fi(1)-fi(2);
end

复制代码

gghhjj · 发表于 2006-12-16 08:43

这些都是matlab程序，相关的参数在注释中都有说明

建议先找一些相关算法的资料看看，把算法搞清楚了才能用

无水1324 · 发表于 2006-12-16 09:46

好东西
感谢分享！

一笑英雄 · 发表于 2006-12-16 15:54

请问，为什么我下载不了？

gghhjj · 发表于 2007-1-7 00:54

原帖由 一笑英雄 于 2006-12-16 15:54 发表
请问，为什么我下载不了？

都是代码，又没有附件，你要下载什么？

无水1324 · 发表于 2007-6-2 13:52

很好的东西！
提升上来，供大家参考学习！

gghhjj · 发表于 2007-6-3 09:01

原帖由 无水1324 于 2007-6-2 13:52 发表

 登录/注册后可看大图

很好的东西！
提升上来，供大家参考学习！

和http://www.chinavib.com/forum/viewthread.php?tid=32137是一样的
只是注释翻译成中文的了

[ 本帖最后由无水1324 于 2007-6-3 11:19 编辑 ]

无水1324 · 发表于 2007-6-3 11:19

原帖由 gghhjj 于 2007-6-3 09:01 发表

 登录/注册后可看大图

和http://www.chinavib.com/forum/viewthread.php?tid=32137是一样的
只是注释翻译成中文的了

哦
我没有仔细看别的版，只是觉得还不错，就提上来了

cr810317 · 发表于 2010-5-13 10:00

好强大，感谢分享:lol

11ella · 发表于 2011-4-16 16:59

不错，很好的帖子。

账号		自动登录	找回密码
密码			我要加入

[分形与混沌] 分形参数计算程序分享&求帮助

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐: