MF-DFA 代码求助

kingjames · 发表于 2018-9-13 10:07:08

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

我在论坛上复制了一个代码，是使用MF-DFA方法计算广义Hurst指数的，但不知道需填入的各个变量的名称，我是零基础，但任务时间很紧，来不及现学，所以想先用来算个数，请各位高手指教，不胜感激！
请问括号中的signal,m,scmin,scmax,ressc,qmin,qmax,qres各自是什么意思，都应该填写什么？最好能举个例子，实在是不胜感激！！！
具体代码如下：

function [s,q,Hq,h,Dh,logFq]=MFDFA(signal,m,scmin,scmax,ressc,qmin,qmax,qres)
%
% Multifractal detrended fluctuation analysis (MFDFA)
%
% [s,q,Hq,h,Dh,logFq]=MFDFA(signal,m,scmin,scmax,ressc,qmin,qmax,qres);
%
% INPUT PARAMETERS-----------------------------------------------------
%
% signal: input signal
% m: polynomial order for the detrending
% scmin: lower bound of the window size s
% scmax: upper bound of the window size s
% ressc: number of elements in s
% qmin lower bound of q
% qmax upper bound of q
% qres number of elements in q
%
% OUTPUT VARIABLES-----------------------------------------------------
%
% s: scale
% q: q-order
% Hq: q-generalized Hurst exponent
% logFq: q-generalized scaling function F(s,q) in log coordinates
% lin_fit linear least square fit to logFq
% h: H鰈der exponent
% Dh: Multifractal spectrum
%
% EXAMPLE-----------------------------------------------------
%
% T=40.96 ; rmin=0.02 ; Dt=rmin/2 ; Law=0 ; ParamLaw1=0.2; ParamLaw2=0; C=1 ; q1=(-10:10)/2;
% [Qr1,time,TtheoQ,q] = IDCnoise_epl(T,rmin,Dt,Law,ParamLaw1,C,q1);
% [Qr2,time,TtheoQ,q] = IDCnoise_epl(T,rmin,Dt,Law,ParamLaw2,C,q1);
% H=0.75 ; oversamprate = 16 ; printout = 0;
% [Ar,VH1,BH] = synthArVH(H,Qr1,Dt,oversamprate,printout) ;
% [Ar,VH2,BH] = synthArVH(H,Qr2,Dt,oversamprate,printout) ;
% [s,q,Hq1,h1,Dh1,logFq1]=MFDFA(diff(VH1),1,10,200,40,-3,3,61);
% [s,q,Hq2,h2,Dh2,logFq2]=MFDFA(diff(VH2),1,10,200,40,-3,3,61);
% figure;
% subplot(311)
% plot(1:4096,diff(VH1),'b-',1:4096,diff(VH2),'r-');xlabel('time');
% ylabel('amplitude'); title('multiplicative cascading noise \DeltaB_H(A(t))and fractional Gaussian noise \DeltaB_H(t)');
% legend('\DeltaB_H(A(t))','\DeltaB_H(t)');
% subplot(323);
% plot(log2(s),logFq1(:,1:5:end),'b-');hold all
% plot(log2(s),logFq2(:,1:5:end),'r-');
% xlabel('log_2(\Deltat)');ylabel('log_2(F_\Delta_t(q))');title('log-scaling function')
% subplot(324)
% plot(q,Hq1,'bo-',q,Hq2,'ro-');xlabel('q'); ylabel('H(q)'); title('q-generalized Hurst exponent');
% legend('\DeltaB_H(A(t))','\DeltaB_H(t)');
% subplot(325)
% plot(h1, Dh1,'bo-',h2, Dh2,'ro-');xlabel('h(q)');ylabel('D(h)');title('multifractal half-spectrum');
% legend('\DeltaB_H(A(t))','\DeltaB_H(t)');
%
% ---------------------------------------------------------------
% Written by Espen A. F. Ihlen (espen.ihlen@ntnu.no), 2009
warning off;
Fluct=cumsum(signal-mean(signal)./std(signal));
FluctRev=fliplr(Fluct);
N=length(Fluct);
ScaleNumb=linspace(log2(scmin),log2(scmax),ressc);
s=round(2.^ScaleNumb);
q=linspace(qmin,qmax,qres);
znumb=find(q==0);
Fq=zeros(length(s),length(q));
for ns=1:length(s),%disp(strcat('computing scale number_',num2str(ns)));
Ns=floor(s(ns)\N);
Var=zeros(Ns,length(q));
Varr=zeros(Ns,length(q));
for v=1:Ns,
SegNumb=((((v-1)*s(ns))+1):(v*s(ns)))';
Seg=Fluct(SegNumb);
SegRev=FluctRev(SegNumb);
poly=polyfit(SegNumb,Seg,m);
polyr=polyfit(SegNumb,SegRev,m);
fit=polyval(poly,SegNumb);
fitr=polyval(polyr,SegNumb);
for nq=1:length(q);
Var(v,nq)=((sum((Seg-fit).^2))/s(ns))^(q(nq)/2);
Varr(v,nq)=((sum((SegRev-fitr).^2))/s(ns))^(q(nq)/2);
end;
clear SegNumb Seg SedRev poly polyr fit fitr
end
for nq=1:length(q)
Fq(ns,nq)=((sum(Var(:,nq))+sum(Varr(:,nq)))/(2*Ns))^(1/q(nq));
end
Fq(ns,znumb)=(Fq(ns,znumb-1)+Fq(ns,znumb+1))./2;
clear Var Varr
end
logFq=log(Fq);
Hq=zeros(1,length(q));
lin_fit=zeros(length(s),length(q));
for nq=1:length(q);
P=polyfit(log(s),logFq(:,nq),1);
lin_fit(1:length(s),nq)=polyval(P,log(s));
Hq(nq)=P(1);
end;
tau=(q.*Hq)-1;
hh=diff(tau)./(q(2)-q(1));
Dh=(q(1:(end-1)).*hh)-tau(1:(end-1));
h=hh-1;

kingjames · 发表于 2018-9-13 10:16:29

奉献点下载，未调试的DFA代码

何佳田力 · 发表于 2019-1-11 14:41:12

正在学习，希望能够搞清楚

黎明的夕阳 · 发表于 2020-7-28 11:47:37

这个程序注释已经很清楚了：（1）signal是输入的数据（如果是气象数据，比如日降雨，就是一列数，数据的总个数有N个）；（2）m是多项式拟合的阶数，就是你划分的小段用几阶的多项式去拟合，一般1-3阶就可以；（3）scmin,scmax这两个指的是你划分的最小最大窗口的值，原文献是m+2~N/4,也有些文献说是最小窗格为10，最大N/10,ressc根据你划分的窗格大小来计算，可以不输入；（4）qmin,qmax,qres对应的q的范围，一般根据数据量的大小确定，-5~5,-10~10......

		自动登录	找回密码
密码			立即注册

[源程序] MF-DFA 代码求助

登录后查看更多精彩内容~