function spectrum=timeseries_spectrum(X,l,Q,h1,h2)
% 此函數(shù)的功能為一個(gè)一維時(shí)間序列的多重分形特征
% 輸入的參數(shù)有時(shí)間序列X;劃分的時(shí)間間段(劃分尺度)數(shù)組l;
% 在計(jì)算配分函數(shù)時(shí)所用到的指數(shù)數(shù)組Q;
% h1表示在通過T(q)與q的關(guān)系圖求a的過程中，每一組q與T(q)值的個(gè)數(shù)，因?yàn)橐鼍€性擬合，h1顯然大于2,一般我們將其取為6到8
% h2表示了從a的最小值mina到最大值maxa的變化步長  一般取0.02到0.1之間的值
% l=[1200,1000,800,600,500,300,200,100,50,25]; %給出了10種尺度,即時(shí)間間隔,單位為天
% Q=(-50:50);                                  %統(tǒng)計(jì)矩函數(shù)中指數(shù)q的取值范圍
% h1=6;h2=0.02;
% 下面來求時(shí)間的冪譜頻率關(guān)系圖(雙對數(shù)關(guān)系圖)
E=[]; m=1;     
% 冪譜數(shù)組
w=[]; n=1;     
%頻率數(shù)組
for n=1:length(X)
    w(n)=1./n;
end
for n=1:length(w)
    e=0;        % 這是一個(gè)中間變量
    for j=1:length(X)
        e=e+X(j)*exp(-i*j*w(n));
    end
    E(m)=(abs(e).^2)./length(X);
    m=m+1;
end
u=log(w);       % 以自然對數(shù)為底
v=log(E);
figure(1),plot(u,v),xlabel('log(\omega)'),ylabel('log(E(\omega))') % 做出冪譜關(guān)于頻率的雙對數(shù)關(guān)系圖
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S=sum(X);       % 計(jì)算時(shí)間序列在總標(biāo)度范圍內(nèi)的和
M=[];           % 存儲(chǔ)同一q值的統(tǒng)計(jì)矩值的數(shù)組
H=[];           % 存儲(chǔ)統(tǒng)計(jì)矩關(guān)于每個(gè)尺度的矩陣
for k=1:length(Q)
    for t=1:length(l)
        n=floor(length(X)/l(t));
        s=0;  %定義一個(gè)臨時(shí)變量
        for r=1:n
            e=(sum(X(1+l(t)*(r-1):l(t)*r))./S).^Q(k);
            s=s+e;
        end
        if rem(length(X),l(t))~=0    % 當(dāng)length(X)不能除盡l(t)時(shí),還需要計(jì)算總區(qū)間剩余部分的統(tǒng)計(jì)矩 
           s=s+(sum(X(l(t)*n+1:length(X)))./S).^Q(k);
        end                          % 第42,43行是十分重要的一步
        M(t)=s;      % 對應(yīng)于某一特定指數(shù)值q的統(tǒng)計(jì)矩
    end
    H(k,:)=M;
end
l=l./length(X);
Logl=log(l);
LogH=log(H);   % 會(huì)出現(xiàn)H中的某一項(xiàng)為負(fù)數(shù)的情況，解決辦法可以將氣溫?cái)?shù)據(jù)修改為絕對溫度
k=floor(length(Q)/8);
figure(2)
for t=1:k:length(Q)
    s=s+1;
    plot(Logl,LogH(t,:)); hold on
    s=s+1;
end
xlabel('log(\lambda)'),ylabel('log(\chi_{q}(\lambda))')
%以上做出統(tǒng)計(jì)矩函數(shù)關(guān)于尺度的雙對數(shù)函數(shù)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 下面做出配分函數(shù)與指數(shù)q的關(guān)系圖,如果是多重分形，則應(yīng)該為一條上凸曲線
T=[];k=1;
for t=1:length(Q)
    Y=LogH(t,:);
    p=polyfit(Logl,Y,1);
    T(k)=p(1);
    k=k+1;
end
figure(3),plot(Q,T,'>'),xlabel('q'),ylabel('\tau(q)')         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 下面來計(jì)算多重分形譜的特征參數(shù)及做出譜曲線的圖像
spectrum=[];s=1;
a=[];           %用于存儲(chǔ)所有的alpha值
X=[];Y=[];
n=1;
maxi=fix(length(Q)/h1);
for t=1:maxi
    X=Q(h1*(t-1)+1:h1*t);
    Y=T(h1*(t-1)+1:h1*t);
    p=polyfit(X,Y,1);
    a(n)=p(1);
    n=n+1;
end
aa=[]; k=0;
for t=1:length(a)
   if a(t)==Inf
      t=t+1;
   else
      k=k+1;
      aa(k)=a(t);
   end
end
amin=min(aa);
amax=max(aa);
spectrum(s)=amax-amin; s=s+1;  % 由此可得一個(gè)特征參數(shù)W
F=[];
m=1;
for A=amin:h2:amax
    b=[];
    l=1;
    for t=1:length(Q)
        b(l)=A*Q(t)-T(t);
        l=l+1;
    end
    F(m)=min(b);
    m=m+1;
end
A=(amin:h2:amax);
f=[];
a=[];
k=1;
for t=1:length(F)
    if F(t)>=0 && F(t)<1
       f(k)=F(t);
       a(k)=A(t);
       k=k+1;
    end
end
spectrum(s)=max(f); s=s+1; % 從而求出第二個(gè)特征參數(shù),即多重分形譜的極大值fmax
figure(4),plot(a,f,'bo');
hold on
plot(a,f);
hold on
number=[];k=1;
for t=1:length(f)
    if f(t)==max(f)
       number(k)=t;
       k=k+1;
    end
end
t=1+fix(length(number)/2);
a_jizhi=a(number(t));
P=polyfit(a,f,2);
A=P(1);
B=P(2)+2*A*a_jizhi;
C=P(3)+B*a_jizhi-A*(a_jizhi)*(a_jizhi);
spectrum(s)=B;              %  從而求得第三個(gè)特征參數(shù)B
u=(amin:0.01:amax);
v=A*(u-a_jizhi).^2+B*(u-a_jizhi)+C;
plot(u,v,'r'),xlabel('\alpha'),ylabel('f(\alpha)')