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曲線')