www.pudn.com > MFDFAmethod.rar > MFDFAmethod.m

function DFA=MFDFAmethod(X,S,q) 
% 此函数为用多重分形DFA方法处理时间序列X 
% 定义一个划分时间序列的间隔数组:S 
% 计算波动函数时所用到的q值数组 
% S=[15,30,90,180,360,730,1100];   
% S=(30:20:1250); 
% q=(-1.6:0.6:2); 
figure(1),plot(X),xlabel('天数'),ylabel('气温(*0.1)'),title('原时间序列的图像'); 
 
% step1:start 
N=length(X); 
averageX=sum(X)./N; 
Y=[]; 
for i=1:N 
    temp=0; 
    for k=1:i 
        temp=temp+(X(k)-averageX);    % 依据公式(1) 
    end 
    Y(i)=temp; 
end 
figure(2),plot(Y),xlabel('天数'),title('经过预处理的时间序列图像'); 
% step1:end 
 
generalHurstExponent=[]; 
for i=1:length(q) 
    generalHurstExponent(i)=GeneralHurstExponent(Y,S,q(i)); 
end 
figure(4),plot(q,generalHurstExponent,'^b'),xlabel('q'),ylabel('h(q)') 
% 以上做出了h(q)与q的关系图 
 
% 下面分析指数q与多重分形标度指数T(q)的关系,可以利用公式:T(q)=q*h(q)-1 
% 具体来说,当T(q)与q满足非线性关系时(在图像上一般是一个上凸函数),则说明对象具有多重分形性 
T=[]; 
for i=1:length(q) 
    T(i)=q(i).*generalHurstExponent(i)-1; 
end 
figure(5),plot(q,T,'^b'),xlabel('q'),ylabel('\tau(q)') 
 
% 下面计算所对应的多重分形谱,可以利用公式:alpha=h(q)+q*h'(q) 与 f(alpha)=q*[alpha-h(q)]+1 
%  
DHurst=diff(generalHurstExponent)./diff(q); 
alpha=[]; f=[]; 
for i=1:length(DHurst) 
    alpha(i)=generalHurstExponent(i)+q(i).*DHurst(i); 
    f(i)=(q(i).^2).*DHurst(i)+1; 
end 
Alpha=[];F=[];j=1; 
for i=1:length(f) 
    if f(i)>0                    %去掉数组f中的负值 
       F(j)=f(i); 
       Alpha(j)=alpha(i); 
       j=j+1; 
    end 
end 
plot(Alpha,F,'bo'),xlabel('\alpha'),ylabel('f(\alpha)') 
DFA=[];s=1; 
DFA(s)=max(Alpha)-min(Alpha);s=s+1; 
DFA(s)=max(F);