tmulti.m

对一维时间序列采用盒计数法对其进行多重分形谱的计算

clear all;
x=rand(1,1024);                   %随机序列（可以换成你自己的序列）
stem(x);                           %显示序列
xlabel('随机序列')                 %显示序列后的标题
n=size(x,2);                       %求序列长度     2，求长度
s=sum(x);                          %求序列的和
Tt=log2(n);                        
e=2.^(1:Tt);                       %两部家一起就是区间长度    
N=n./e;                            %求区间个数
for k=1:Tt
    d=N(k);
    for m=1:e(k)
        A=x(d*(m-1)+[1:d]);
        s0=sum(A);
        pe(m,k)=s0/s;
    end
end                                %求概率

q=-10:1:10;    %q取-10~10 配分函数Xq和概率之间的关系
for q1=1:21
for c=1:Tt
Xq(q1,c)=0;
for i=1:e(c)
Xq(q1,c)=Xq(q1,c)+pe(i,c)^q(q1);    %根据公式 Xq=p的q次密求和
end;
end;
end;
                                             %求t套(q)
lnXq=log(Xq);
lne=log(1./e);
figure,plot(lne,lnXq,'r-*')                              %求lnXp和lne的图像
xlabel('lnxq-lne曲线')

for c=1:Tt
for q1=1:21
tq(q1,c)=lnXq(q1,c)/lne(c); 
end;
end; 
                             
for q1=1:21                    %此两段加在一起是求分形普fa
xp(q1,Tt)=0;
for i=1:e(c)
xp(q1)=xp(q1)+log(pe(i,Tt))*pe(i,Tt)^q(q1);
end;
end;
for q1=1:21
aai(q1)=0;
aai(q1)=(xp(q1)/Xq(q1,Tt))/log(4);
fa(q1)=0;
fa(q1)=-aai(q1)*q(q1)+log(Xq(q1,Tt))/log(4); 
end;

figure,plot(aai,fa,'b-*')
xlabel('f(a)-a曲线')