timeseries_spectrum.m

此函数的功能为一个一维时间序列的多重分形特征

function spectrum=timeseries_spectrum(X,l,Q,h1,h2)
% 此函数的功能为一个一维时间序列的多重分形特征
% 输入的参数有时间序列X;划分的时间间段(划分尺度)数组l;
% 在计算配分函数时所用到的指数数组Q;
% h1表示在通过T(q)与q的关系图求a的过程中，每一组q与T(q)值的个数，因为要做线性拟合，h1显然大于2,一般我们将其取为6到8
% h2表示了从a的最小值mina到最大值maxa的变化步长  一般取0.02到0.1之间的值
% l=[1200,1000,800,600,500,300,200,100,50,25]; %给出了10种尺度,即时间间隔,单位为天
% Q=(-50:50);                                  %统计矩函数中指数q的取值范围
% h1=6;h2=0.02;
% 下面来求时间的幂谱频率关系图(双对数关系图)
E=[]; m=1;     
% 幂谱数组
w=[]; n=1;     
%频率数组
for n=1:length(X)
    w(n)=1./n;
end
for n=1:length(w)
    e=0;        % 这是一个中间变量
    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);       % 以自然对数为底
v=log(E);
figure(1),plot(u,v),xlabel('log(\omega)'),ylabel('log(E(\omega))') % 做出幂谱关于频率的双对数关系图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S=sum(X);       % 计算时间序列在总标度范围内的和
M=[];           % 存储同一q值的统计矩值的数组
H=[];           % 存储统计矩关于每个尺度的矩阵
for k=1:length(Q)
    for t=1:length(l)
        n=floor(length(X)/l(t));
        s=0;  %定义一个临时变量
        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    % 当length(X)不能除尽l(t)时,还需要计算总区间剩余部分的统计矩 
           s=s+(sum(X(l(t)*n+1:length(X)))./S).^Q(k);
        end                          % 第42,43行是十分重要的一步
        M(t)=s;      % 对应于某一特定指数值q的统计矩
    end
    H(k,:)=M;
end
l=l./length(X);
Logl=log(l);
LogH=log(H);   % 会出现H中的某一项为负数的情况，解决办法可以将气温数据修改为绝对温度
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))')
%以上做出统计矩函数关于尺度的双对数函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 下面做出配分函数与指数q的关系图,如果是多重分形，则应该为一条上凸曲线
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)')         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 下面来计算多重分形谱的特征参数及做出谱曲线的图像
spectrum=[];s=1;
a=[];           %用于存储所有的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;  % 由此可得一个特征参数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; % 从而求出第二个特征参数,即多重分形谱的极大值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;              %  从而求得第三个特征参数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)')