www.pudn.com > timeseries_spectrum.rar > 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)')