分形参数计算程序分享.doc

分形参数计算程序分享计算hurst指数CODE: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 nargin1 error(时间序列无效.); end x=x(:); % N 是时间序列的长度 N=length(x);endif nargin1 error(n 必须是一个变化的标量或矢量.); end % n 必须是个整数 if n-round(n)=0 error(n 必须是个整数.); end % n 必须是确定 if n=0 error(n 必须是确定.); endendif nargin2 error(q 必须是标量.); end % q 必须是整数 if q-round(q)=0 error(q 必须是整数.); end % q 必须是确定 if q0 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 k0 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 k0 & k0 & kq/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 nargout1 % 开始计算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 nargout2 % 估算 V V(i,1)=mean(rs)/sqrt(n(i); endend计算lyapunovCODE:function lambda_1=lyapunov_wolf(data,N,m,tau,P)%该函数用来计算时间序列的最大Lyapunov 指数-Wolf 方法%m: 嵌入维数%tau:时间延迟%data:时间序列%N:时间序列长度%P:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|IJ|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 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_epsmax_eps中找不到演化相点时，对max_eps的放宽幅度 min_eps = min_d + dlt_eps / 2 ; %演化相点与当前相点距离的最小限 max_eps = min_d + 2 * dlt_eps; %&演化相点与当前相点距离的最大限 % 从P+1M-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 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 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 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 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算法计算关联维数Copy to clipboardCODE: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;endfid=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);计算kolmogorovCopy to clipboardCODE:clear allx=load(are.dat)n=length(x)sd=