分形参数计算程序分享文档格式.docx
《分形参数计算程序分享文档格式.docx》由会员分享，可在线阅读，更多相关《分形参数计算程序分享文档格式.docx（9页珍藏版）》请在冰豆网上搜索。
end
x=x（:
%N是时间序列的长度
N=length（x）;
2|isempty（n）==1
n=1;
%n必须是一个变化的标量或矢量
ifmin（size（n））>
n必须是一个变化的标量或矢量.'
%n必须是个整数
ifn-round（n）~=0
error（'
n必须是个整数.'
%n必须是确定
ifn<
=0
n必须是确定.'
4|isempty（q）==1
q=0;
ifq=='
t=autocorr（x,1）;
t=t
（2）;
q=（（3*N/2）^（1/3））*（2*t/（1-t^2））^（2/3）;
else
%q必须是标量
ifsum（size（q））>
2
q必须是标量.'
%q必须是整数
ifq-round（q）~=0
q必须是整数.'
%q必须是确定
ifq<
q必须是确定.'
end
fori=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）;
%估算这个标准偏差
switchmethod
case'
%Hurst-Mandelbrot参数
stdev=std（X）;
%Lo参数
forj=1:
a
sq=0;
fork=0:
q
v（k+1）=sum（X（k+1:
n（i）,j）'
*X（1:
n（i）-k,j））/（n（i）-1）;
ifk>
sq=sq+（1-k/（q+1））*v（k+1）;
stdev（j）=sqrt（v
（1）+2*sq）;
%Moody-Wu参数
sq1=0;
sq2=0;
sq1=sq1+（1-k/（q+1））*（n（i）-k）/n（i）/n（i）;
sq2=sq2+（1-k/（q+1））*v（k+1）;
stdev（j）=sqrt（（1+2*sq1）*v
（1）+2*sq2）;
%Parzen参数
ifmod（q,2）~=0
在"
Parzen"
参数中q必须是2.'
0&
k<
=q/2
sq1=sq1+（1-6*（k/q）^2+6*（k/q）^3）*v（k+1）;
elseifk>
k>
q/2
sq2=sq2+（1-（k/q）^3）*v（k+1）;
stdev（j）=sqrt（v
（1）+2*sq1+2*sq2）;
otherwise
你应该付给"
method"
另一个值.'
%估算（R（t）/s（t））
rs=（max（cumdev）-min（cumdev））./stdev;
clearstdev
%取这个平均值R/S的对数
logRS（i,1）=log10（mean（rs））;
ifnargout>
%开始计算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））;
%估算V
V（i,1）=mean（rs）/sqrt（n（i））;
计算lyapunov
functionlambda_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;
%重构相空间中相点的个数
fori=1:
（M-1）
forj=i+1:
M
d=0;
fork=1:
m
d=d+（Y（k,i）-Y（k,j））*（Y（k,i）-Y（k,j））;
d=sqrt（d）;
ifmax_d<
d
max_d=d;
ifmin_d>
min_d=d;
avg_dd=avg_dd+d;
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个相点对应的最近距离点的下标
fori=（P+1）:
（M-1）
%限制短暂分离，从点P+1开始搜索
d=d+（Y（k,i）-Y（k,1））*（Y（k,i）-Y（k,1））;
if（d<
DK）&
（d>
min_eps）
DK=d;
Loc_DK=i;
以下计算各相点对应的李氏数保存到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）的累计和
fori=2:
（M-1）
%计算演化距离
DK1=0;
DK1=DK1+（Y（k,i）-Y（k,Loc_DK+1））*（Y（k,i）-Y（k,Loc_DK+1））;
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）;
lmd（i-1）=sum_lmd/（i-1）;
以下寻找i点的最短距离：
要求距离在指定距离范围内尽量短，与DK1的角度最小
point_num=0
%&
在指定距离范围内找到的候选相点的个数
cos_sita=0
夹角余弦的比较初值——要求一定是锐角
zjfwcs=0
增加范围次数
while（point_num==0）
%*搜索相点
forj=1:
ifabs（j-i）<
=（P-1）
候选点距当前点太近，跳过！
continue;
%*计算候选点与当前点的距离
dnew=0;
fork=1:
dnew=dnew+（Y（k,i）-Y（k,j））*（Y（k,i）-Y（k,j））;
dnew=sqrt（dnew）;
if（dnew<
min_eps）|（dnew>
max_eps）
不在距离范围，跳过！
continue;
%*计算夹角余弦及比较
DOT=0;
DOT=DOT+（Y（k,i）-Y（k,j））*（Y（k,i）-Y（k,old_Loc_DK+1））;
CTH=DOT/（dnew*DK1）;
ifacos（CTH）>
（3.14151926/4）
不是小于45度的角，跳过！
ifCTH>
cos_sita
新夹角小于过去已找到的相点的夹角，保留
cos_sita=CTH;
Loc_DK=j;
DK=dnew;
point_num=point_num+1;
end
ifpoint_num<
=min_point
max_eps=max_eps+dlt_eps;
zjfwcs=zjfwcs+1;
ifzjfwcs>
MAX_CISHU
超过最大放宽次数，改找最近的点
forii=1:
ifabs（i-ii）<
=（P-1）
d=0;
d=d+（Y（k,i）-Y（k,ii））*（Y（k,i）-Y（k,ii））;
d=sqrt（d）;
if（d<
DK=d;
Loc_DK=ii;
break;
;
扩大距离范围后重新搜索
cos_sita=0;
%取平均得到最大李雅普诺夫指数
lambda_1=sum（lmd）/length（lmd）;
G-P算法计算关联维数
[Copytoclipboard]
function[ln_r,ln_C]=G_P（data,N,tau,min_m,max_m,ss）
%此程序是用G-P算法计算关联维数Dc
%tau是固定时间间隔
%min_m最小的嵌入维数
%max_m最大的嵌入维数
%ssr的步长
form=min_mmax_m
%重建矢量空间
M=N-（m-1）tau;
%矢量空间的点数
fori=1M-1
forj=i+1M
d（i,j）=max（abs（Y（,i）-Y（,j）））;
%计算其余点到点Xi的距离
max_d=max（max（d））;
%是所有点中距离最远的点
d（1,1）=max_d;
min_d=min（min（d））;
%是所有点中距离最近的点
delt=（max_d-min_d）ss;
%是r的步长
fork=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）;
plot（ln_r（m,）,ln_C（m,））;
holdon;
fid=fopen（'
lnr.txt'
'
w'
fprintf（fid,'
%6.2f%6.2fn'
ln_r）;
fclose（fid）;
fid=fopen（'
lnC.txt'
ln_C）;
计算kolmogorov
clearall
x=load（'
are.dat'
）
n=length（x）
sd=std（x）
r=0.2*sd
forii=1:
m=ii+1;
num=zeros（n-m+1,1）;
n-m+1
ifj~=i
fork=1:
m
d（k）=abs（x（i+k-1）-x（j+k-1））;
d1=max（d）;
ifd1<
r
num（i,1）=num（i,1）+1;
end
c0（i）=num（i,1）/（n-m）;
c1（i）=log（c0（i））;
sc=sum（c1）;
fi（ii）=sc/（n-m+1）;
app=fi
（1）-fi
（2）;
升级会员 