%這里設采樣信號為leleccum
load leleccum;
%取其中的1024個點
signal=leleccum(1:1024);
s=sum(signal);
%y為101*11的二維數組，每一行對應不同q值的配分函數
y=zeros(101,11);
i=1;
%q取值范圍為-50~50，步進為1
for q=-50:50
    j=1;
    %以下循環用于求q相同而測量單元尺寸不同的配分函數，存于y的第i行中
    for n=0:10
        p=zeros(1,1024/2^n);
        m=1;
        %以下的循環用于求測量單元尺寸為2^n時，各個盒子的平均高度概率分布，存于數組p中
        for k=1:2^n:1024
            Ins=signal(k:2^n+k-1);
            p(m)=sum(Ins)/s;
            m=m+1;
        end;
        y(i,j)=sum(p.^q);
        j=j+1;
    end;
    %x為測量單元尺寸
    x=zeros(1,11);
    for n=1:11
        x(n)=2^n/1024;
    end
    n=1:11;
    xlabel('測量單元尺寸的自然對數');ylabel('配分函數的自然對數');
    plot(log(x),log(y(i,n)));
    hold on;
    i=i+1;
end;
Dy=diff(log(y),1,2);
De=diff(log(x),1,2);
%Tq為質量指數，是q的函數
Tq=zeros(101,1);
for n=1:101
    Tq(n)=mean(Dy(n,:)./De);
end;
Tq=Tq';
%廣義分形維數Dq
q=-50:1:50;
Dq=zeros(101,1);
for i=1:1:101
    Dq(i)=Tq(i)/(q(i)-1)
end;
figure;
xlabel('質量指數');ylabel('廣義分形維數');
plot(q,Dq);
q=-49:49;
z=diff(Tq,1,2);
%Alf為奇異指數
Alf=1:99;
for n=1:99
    Alf(n)=(z(n)+z(n+1))/2;
end;
Tq=Tq(2:100);
%fa為多重分形譜函數
fa=Alf.*q-Tq;
n=1:99;
figure;
xlabel('a');ylabel('f(a)');
plot(Alf,fa);