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