分形参数计算程序分享.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 nargin1error('时间序列无效.');endx=x(:);% N 是时间序列的长度N=length(x);endif nargin1error('n 必须是一个变化的标量或矢量.');end% n 必须是个整数if n-round(n)~=0error('n 必须是个整数.');end% n 必须是确定if n2error('q 必须是标量.');end% q 必须是整数if q-round(q)~=0error('q 必须是整数.');end% q 必须是确定if q0sq=sq+(1-k/(q+1))*v(k+1);endendstdev(j)=sqrt(v(1)+2*sq);endcase 'MW'% Moody-Wu 参数for j=1:asq1=0;sq2=0;for k=0:qv(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);if k>0sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);sq2=sq2+(1-k/(q+1))*v(k+1);endendstdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);endcase 'Parzen'% Parzen 参数if mod(q,2)~=0error('在"Parzen" 参数中 q 必须是 2.');endfor j=1:asq1=0;sq2=0;for k=0:qv(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);if k>0 & k0 & k>q/2sq2=sq2+(1-(k/q)^3)*v(k+1);endendstdev(j)=sqrt(v(1)+2*sq1+2*sq2);endotherwiseerror('你应该付给 "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));endif nargout>2% 估算 VV(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，则演化相点只能在|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 : Md = 0;for k = 1 : md = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd = sqrt(d);if max_d dmin_d = d;endavg_dd = avg_dd + d;endendavg_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)及其最短距离 DKDK = 1.0e+100; %第 i个相点到其最近距离点的距离Loc_DK = 2; %第 i个相点对应的最近距离点的下标 for i = (P+1):(M-1) %限制短暂分离，从点 P+1开始搜索d = 0;for k = 1 : md = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));endd = sqrt(d);if (d min_eps) DK = d;Loc_DK = i;endend% 以下计算各相点对应的李氏数保存到 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 : mDK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));endDK1 = 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);endlmd(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) max_eps ) %&&不在距离范围，跳过！continue; end%*计算夹角余弦及比较DOT = 0;for k = 1 : mDOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));endCTH = DOT/(dnew*DK1);if acos(CTH) > (3.14151926/4) %&&不是小于 45度的角，跳过！continue;endif CTH > cos_sita %&&新夹角小于过去已找到的相点的夹角，保留cos_sita = CTH;Loc_DK = j;DK = dnew;endpoint_num = point_num +1;end if point_num MAX_CISHU %&&超过最大放宽次数，改找最近的点DK = 1.0e+100;for ii = 1 : (M-1)if abs(i-ii) min_eps) DK = d;Loc_DK = ii;endendbreak; endpoint_num = 0 ; %&&扩大距离范围后重新搜索cos_sita = 0;endendend%取平均得到最大李雅普诺夫指数lambda_1=sum(lmd)/length(lmd);G-P算法计算关联维数[Copy to clipboard]CODE: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_mY=reconstitution(data,N,m,tau);%重建矢量空间M=N-(m-1)tau;%矢量空间的点数for i=1M-1for j=i+1Md(i,j)=max(abs(Y(,i)-Y(,j)));%计算其余点到点 Xi的距离 end endmax_d=max(max(d));%是所有点中距离最远的点d(1,1)=max_d;min_d=min(min(d));%是所有点中距离最近的点delt=(max_d-min_d)ss;%是 r的步长for k=1ssr=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);%lnrfprintf('%d%d%d%dn',k,ss,m,max_m);endplot(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');fpr。