matlab 计算DFN（去趋势波动分析） regress回归分析

梦里梦外; 2022-12-26 11:23 73阅读 0赞

实现方法来自论文：

**Quantification of scaling exponents and crossover phenomena in nonstationary**

**heartbeat time series**

C.‐K. Peng ,  Shlomo Havlin ,  H. Eugene Stanley , and  Ary L. Goldberger

输入：一个数组B，长度为len

输出：斜率b（2）

第一步：计算y(k),其中Bave位B的均值。

![20201205142132790.png][]

第二步：计算yn(k)，n是能被len整除的数，比如我这是3500个数，那么n就可以取1，2，5，10，35等。

根据n的大小划分，比如我取n=500，那么y就被划分成了7块，对应了1-500；501-1000；1001-1500；1501-2000；2001-2500；2501-3000；3001-3500。

然后计算每一块的regress（即利用回归分析得到预测值）。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1RoZVdpc2VfbHp5_size_16_color_FFFFFF_t_70][]

第三步：根据如下公式计算F（n）

![20201205142146868.png][]

第四步：得到n和F（n），以log10为底，然后画图，再次用regress得到斜率b。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1RoZVdpc2VfbHp5_size_16_color_FFFFFF_t_70 1][]

function out=dfa(a)
        avg=mean(a);
        sum=0.0;
        len=length(a);
        for i=1:3500
            sum=sum+a(i)-avg;
            y(i)=sum;
        end
        %x=1:3500;
        %X=[ ones(length(y),1), x' ];
        %b=regress(y',X)
        index=1;
        for n=10:3500
            if(rem(3500,n)==0)    
                count=0;
                k=3500/n;
                for i=1:k
                    begin=(i-1)*n+1;
                    finish=i*n;
                    q=y(begin:finish);
                    p=begin:finish;
                    X=[ ones(length(q),1), p' ];
                    b=regress(q',X); %回归 y=b1+b2*x;
                    for j=begin:finish
                        cc=b(1)+b(2)*j;
                        op=y(j)-cc;
                        count=count+op*op;
                    end
                end
                f=(count/3500)^0.5;
                K(index)=log10(n);
                F(index)=log10(f);
                index=index+1;
            end
        end
        
        X=[ ones(length(F),1), K' ];
        [b,bint,r,rint,stats]=regress(F',X);
        out=b(2);
        z=b(1)+b(2)*K;
        figure,plot(K,F,'rp',K,z,'b');
        title('散点图及回归线');
    end

[20201205142132790.png]: /images/20221120/dcd2f910190c4f238faf3f86b9187614.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1RoZVdpc2VfbHp5_size_16_color_FFFFFF_t_70]: /images/20221120/a8019030408a4714bb0f54d158f65795.png
[20201205142146868.png]: https://img-blog.csdnimg.cn/20201205142146868.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1RoZVdpc2VfbHp5_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/202012051429553.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1RoZVdpc2VfbHp5,size_16,color_FFFFFF,t_70