mfdfamethod.m

计算多重分形非趋势波动分析matlab程序

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);