分形参数计算程序分享

1、分形参数计算程序分享计算 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 的默

2、认值 .if nargin1error( 时间序列无效 .); endx=x(:);% N 是时间序列的长度N=length(x);endif nargin1error(n 必须是一个变化的标量或矢量 .); end% n 必须是个整数if n-round(n)=0error(n 必须是个整数 .); end% n 必须是确定if n=0error(n 必须是确定 .);endendif nargin2error(q必须是标量 .);end% q 必须是整数if q-round(q)=0error(q必须是整数 .);end% q 必须是确定if q0sq=sq+(1-k/(q+1)*v(k+

3、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 k0 sq1=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

4、参数中 q 必须是 2.); end 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 k0 & k0 & kq/2sq2=sq2+(1-(k/qF3)*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 的对数

5、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);endif nargout2% 估算 VV(i,1)=mean(rs)/sqrt(n(i);endend

6、计算 lyapunovCODE:function lambda_1=lyapunov_wolf(data,N,m,tau,P) % 该函数用来计算时间序列的最大 % m: 嵌入维数% tau: 时间延迟% data: 时间序列% N: 时间序列长度% P: 时间序列的平均周期J|P 的相点中搜寻Lyapunov 指数 -Wolf 方法, 选择演化相点距当前点的位置差，即若当前相点为I ，则演化相点只能在 |I% lambda_1: 返回最大 lyapunov 指数值 min_point=1 ; %& 要求最少搜索到的点数 MAX_CISHU=5 ; %&最大增加搜索范围次数 %FLYINGHA

7、WK最大相点距离最小相点距离相空间重构 重构相空间中相点的个数% 求最大、最小和平均相点距离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;end平均相点距离endav

8、g_d = 2*avg_dd/(M*(M-1);%dlt_eps = (avg_d - min_d) * 0.02 ;%max_eps的放宽幅度min_eps = min_d + dlt_eps / 2 ;%max_eps = min_d + 2 * dlt_eps ;%&若在min_epsmax_eps中找不到演化相点时，对演化相点与当前相点距离的最小限演化相点与当前相点距离的最大限DK = 1.0e+100;Loc_DK = 2;%从P+1M-1个相点中找与第一个相点最近的相点位置（Loc_DK）及其最短距离DK第 i 个相点到其最近距离点的距离 第 i 个相点对应的最近距离点的下标for

9、 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); end d = 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个Iog2

10、 (DK1/DK的累计和用于求i点的lambda值summd = 0 ;%存放前i个Iog2 ( 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); oId_Loc_DK = Loc_DK ;%保存原最近位置相点oId_DK=DK;% 计算前i个log2 (DK1/DK)的累计和以及保存i点的李氏指数 if (DK1 = 0)&( DK = 0)sum_Imd = sum_Imd +

11、Iog(DK1/DK) /Iog(2);endImd(i-1) = sum_Imd/(i-1);% 以下寻找i点的最短距离：要求距离在指定距离范围内尽量短，与DK1的角度最小point_num = 0 ; % &在指定距离范围内找到的候选相点的个数cos_sita = 0 ; %&夹角余弦的比较初值 要求一定是锐角zjfwcs=0 ;%&增加范围次数whiIe (point_num = 0)% *搜索相点for j = 1 : (M-1)if abs(j-i) =(P-1)%&候选点距当前点太近，跳过！continue;end%*计算候选点与当前点的距离dnew = 0;for k = 1 :

12、 mdnew = dnew + (Y(k,i)-Y(k,j)*(Y(k,i)-Y(k,j);enddnew = sqrt(dnew);if (dnew 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

13、 CTH cos_sita %& cos_sita = CTH; Loc_DK = j;DK = dnew;end point_num = point_num +1;endif point_num MAX_CISHU %& 超过最大放宽次数，改找最近的点DK = 1.0e+100;for ii = 1 : (M-1)if abs(i-ii) = (P-1) %& 候选点距当前点太近，跳过！ continue;endd = 0;for k = 1 : md = d + (Y(k,i)-Y(k,ii)*(Y(k,i)-Y(k,ii);endd = sqrt(d);if (d min_eps) DK

14、 = d;Loc_DK = ii;endendbreak;扩大距离范围后重新搜索endpoint_num = 0 ; %& cos_sita = 0;endendend %取平均得到最大李雅普诺夫指数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算法计算关联维数 De% N 是时间序列的长度% tau 是固定时间间隔% min_ m最小的嵌入维数% max_n最大的嵌入维数% ssr 的步长for m

15、=min_mmax_mY=reeonstitution(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 的距离endend max_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)=eorrelation_integral(Y,M,r)

16、;%计算关联积分函数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);endplot(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);计算 kolmogorovCopy to clipboard CODE: clear all x=load(are.dat) n=length(x) sd=std(x) r=0.2*sd for ii