DFA（去趋势波动分析）法计算程序

hillside

         DFA（去趋势波动分析）法已经在气候、水文、地质等领域有了一些应用，在其他多个学科也都有广泛的应用价值。
         先从科学网博客转一个贴，页面下方的“源代码下载”是科学网博主于同奎提供的程序，然后合并2个其他来源的程序，制作成附件。这样就有3个matlab 程序进行对比、分析了。   
        本帖以下内容完全复制自于同奎的网页，特此说明。
                                                             DFA(Detrended fluctuation analysis) matlab 程序
                                                                   http://blog.sciencenet.cn/blog-4716-46086.html

The method of detrended fluctuation analysis has proven useful in revealing the extent of long-range correlations in time series.

DFA 是检测时间序列长程相关的范围的有效方法。

网上有很多源代码。但Matlab的代码不多。我从mathworks下载了一个，并改造了一番，分享来用。

附件中三个文件，一个是DFA.m，是用来计算对应给定时间间隔n的F(n)值；另一个是runDFA.m，用来执行计算的，其中调用了前面的DFA函数，这个要根据你的需要自己稍作修改，其中我在这里将时间取值设置为log坐标下等距，让画出来的图形更好看；另一个是data.mat，是演示数据。

进入包含着三个文件的目录，在Matlab命令窗口，直接输入runDFA，就可以得到结果。

------------------------------

runDFA.m

------------------------------

function F_n=DFA(DATA,win_length,order)
      
       N=length(DATA);   
       n=floor(N/win_length);
       N1=n*win_length;
       y=zeros(N1,1);
       Yn=zeros(N1,1);
      
       fitcoef=zeros(n,order+1);
            mean1=mean(DATA(1:N1));
       for i=1:N1
           y(i)=sum(DATA(1:i)-mean1);
       end
       y=y';
       for j=1:n
           fitcoef(j,:)=polyfit(1:win_length,y(((j-1)*win_length+1):j*win_length),order);
       end
       for j=1:n
           Yn(((j-1)*win_length+1):j*win_length)=polyval(fitcoef(j,:),1:win_length);
       end
       sum1=sum((y'-Yn).^2)/N1;
       sum1=sqrt(sum1);
       F_n=sum1;
         

----------------------------------------------------------------------------------

runDFA.m

------------------------------

%导入演示数据
load data

%设置时间
t=1:0.05:3;
n=zeros(1,length(t));
for i=1:length(t)
n(i)=10^t(i);
end
n=floor(n);
n=n';

%初始化
len=length(n);
F_n=zeros(len,1);

%对每个n值求F_n
for i=1:len
     F_n(i)=DFA(data,n(i),1);
end

%线性拟合
p=polyfit(log10(n),log10(F_n),1);

%画图
plot(n,F_n(:,1),'o');

%画拟合直线
x=n;
y=p(2)+p(1)*log10(x);
for i=1:len
   y(i)=10^y(i);
end
hold on
plot(x,y,'black');

%设置坐标
axis([8 1200 0.003 0.5])
set(gca,'XSCALE','log');
set(gca,'YSCALE','log');

xlabel('n','FontSize',16,'FontAngle','Italic')
ylabel('F(n)','FontSize',16,'FontAngle','Italic')

%添加文字
str=['

F(n)=′num2str(10p(2))′×n′num2str(p(1))′

']
text('Interpreter','latex','String',str,'Position',[50 0.01],'FontSize',16)

hold off

----------------------------------------------------------------------------------

                               
登录/注册后可看大图

源代码下载

kongfeng0824

感谢分享                          

rceclx

LZ真强，谢谢！

黎明的夕阳

正需要，谢谢楼主

tanker517

谢谢楼主分享，很有用

zhu

谢谢分享

scwjx520

好东西，以后可以研究{:5_213:}

黎明的夕阳

感谢分享  

lgsh1988

谢谢楼主分享！

lss_lucky

感谢分享，还是需要琢磨一下怎么用……