MFDFA1_时间序列_趋势_MF-DFAmatlab_

function [Hq,tq,hq,Dq,Fq] = MFDFA1(signal,scale,q,m,Fig) % Multifractal detrended fluctuation analysis (MFDFA) % % [Hq,tq,hq,Dq,Fq]=MFDFA(signal,scale,q,m,Fig); % % INPUT PARAMETERS--------------------------------------------------------- % % signal: input signal % scale: vector of scales % q: q-order that weights the local variations % m: polynomial order for the detrending % Fig: 1/0 flag for output plot of Fq, Hq, tq and multifractal % spectrum (i.e. Dq versus hq). % % OUTPUT VARIABLES--------------------------------------------------------- % % Hq: q-order Hurst exponent % tq: q-order mass exponent % hq: q-order singularity exponent % Dq: q-order dimension % Fq: q-order scaling function % % EXAMPLE------------------------------------------------------------------ % % load fractaldata % scmin=16; % scmax=1024; % scres=19; % exponents=linspace(log2(scmin),log2(scmax),scres); % scale=round(2.^exponents); % q=linspace(-5,5,101); % m=1; % signal1=multifractal; % signal2=monofractal; % signal3=whitenoise; % [Hq1,tq1,hq1,Dq1,Fq1]=MFDFA1(signal1,scale,q,m,1); % [Hq2,tq2,hq2,Dq2,Fq2]=MFDFA1(signal2,scale,q,m,1); % [Hq3,tq3,hq3,Dq3,Fq3]=MFDFA1(signal3,scale,q,m,1); %-------------------------------------------------------------------------- warning off X=cumsum(signal-mean(signal)); if min(size(X))~=1||min(size(scale))~=1||min(size(q))~=1; error('Input arguments signal, scale and q must be a vector'); end if size(X,2)==1; X=transpose(X); end if min(scale)<m+1 error('The minimum scale must be larger than trend order m+1') end for ns=1:length(scale), segments(ns)=floor(length(X)/scale(ns)); for v=1:segments(ns), Index=((((v-1)*scale(ns))+1):(v*scale(ns))); C=polyfit(Index,X(Index),m); fit=polyval(C,Index); RMS_scale{ns}(v)=sqrt(mean((X(Index)-fit).^2)); end for nq=1:length(q), qRMS{nq,ns}=RMS_scale{ns}.^q(nq); Fq(nq,ns)=mean(qRMS{nq,ns}).^(1/q(nq)); end Fq(q==0,ns)=exp(0.5*mean(log(RMS_scale{ns}.^2))); end for nq=1:length(q), C = polyfit(log2(scale),log2(Fq(nq,:)),1); Hq(nq) = C(1); qRegLine{nq} = polyval(C,log2(scale)); end tq = Hq.*q-1; hq = diff(tq)./(q(2)-q(1)); Dq = (q(1:end-1).*hq)-tq(1:end-1); %OUTPUT FIGURE------------------------------------------------------------- if Fig==1, qindex=[1,round(length(q)/2),length(q)]; qindex2=[1,round((length(q)-1)/2),length(q)-1]; %Variable settings--------- qstart1=q(qindex(1)); qmid1=q(qindex(2)); qstop1=q(qindex(3)); Hqstart1=Hq(qindex(1)); Hqmid1=Hq(qindex(2)); Hqstop1=Hq(qindex(3)); tqstart1=tq(qindex(1)); tqmid1=tq(qindex(2)); tqstop1=tq(qindex(3)); hqstart2=hq(qindex2(1)); hqmid2=hq(qindex2(2)); hqstop2=hq(qindex2(3)); Dqstart1=Dq(qindex2(1)); Dqmid1=Dq(qindex2(2)); Dqstop1=Dq(qindex2(3)); for nq=1:length(qindex), qRegFit(nq,:)=qRegLine{qindex(nq)}; end X1=log2(scale); YMatrix1=[log2(Fq(qindex,:));qRegFit]; X2=q; Y1=Hq; Y3=tq; X4=hq; Y5=Dq; %--------------------------- q_end=num2str(max(q)); if length(find(q==0))==1 q_middle=0; else qm0=min(q)+((max(q)-min(q))/2); qm1=[q(find(q<qm0, 1, 'last' )),q(find(q>qm0, 1 ))]; qm2=qm1-qm0; midind= qm2==min(qm2); q_middle=qm2(midind(1)); end q_mid=num2str(q_middle); q_start=num2str(min(q)); % Create figure figure1 = figure('PaperSize',[20.98 29.68],'Color',[1 1 1]); scaleInd=floor(min(log2(scale))):ceil(max(log2(scale))); for ns=1:length(scaleInd), scaletick{ns}=num2str(2^scaleInd(ns)); end Fqind=floor(log2(min(min(Fq)))):ceil(log2(max(max(Fq)))); for nf=1:length(Fqind), Fqtick{nf}=num2str(Fqind(nf)); end % Create subplot subplot1 = subplot(2,2,1,'Parent',figure1,... 'YTickLabel',Fqtick,... 'XTickLabel',scaletick,... 'XTick',scaleInd,... 'LineWidth',2,... 'FontSize',14); hold(subplot1,'all'); % Create multiple lines using matrix input to plot plot1 = plot(X1,YMatrix1,'Parent',subplot1); set(plot1(1),'MarkerFaceColor',[1 0 0],'MarkerEdgeColor',[0 0 0],... 'Marker','o',... 'LineStyle','none',... 'Color',[1 0 0],... 'DisplayName',strcat('q = ',num2str(min(q)))); set(plot1(2),'MarkerFaceColor',[0 0 1],'MarkerEdgeColor',[0 0 0],... 'Marker','o',... 'LineStyle','none',... 'Color',[0 0 1],... 'DisplayName',strcat('q = ',num2str(q_middle))); set(plot1(3),'MarkerFaceColor',[0 0.498 0],'MarkerEdgeColor',[0 0 0],... 'Marker','o',... 'LineStyle','none',... 'Color',[0 0.498 0],... 'DisplayName',strcat('q = ',num2str(max(q)))); set(plot1(4),'LineWidth',2,'Color',[1 0 0],'DisplayName',strcat('q = ',num2str(min(q)))); set(plot1(5),'MarkerFaceColor',[1 0 0],'MarkerEdgeColor',[1 1 1],... 'LineWidth',2,... 'Color',[0 0 1],... 'DisplayName',strcat('q = ',num2str(q_middle))); set(plot1(6),'LineWidth',2,'Color',[0 0.498 0],'DisplayName',strcat('q = ',num2str(max(q)))); % Create xlabel xlabel('Scale (segment sample size)','FontSize',14); % Create ylabel ylabel('log_2(Fq)','FontSize',14); % Create title title('Scaling function Fq (q-order RMS)','FontSize',14); qInd=floor(min(q)):ceil(max(q)); for nq=1:length(qInd), qtick{nq}=num2str(qInd(nq)); end % Create subplot subplot2 = subplot(2,2,2,'Parent',figure1,... 'XTickLabel',qtick,... 'XTick',qInd,... 'LineWidth',2,... 'FontSize',14); hold(subplot2,'all'); % Create plot plot(X2,Y1,'Parent',subplot2,'Marker','o','LineStyle','none',... 'Color',[0 0 1]); % Create xlabel xlabel('q','FontSize',16); % Create ylabel ylabel('Hq','FontSize',16); % Create title title('q-order Hurst exponent','FontSize',14); % Create plot plot(qstart1,Hqstart1,'Parent',subplot2,'MarkerFaceColor',[1 0 0],... 'MarkerEdgeColor',[0 0 0],... 'MarkerSize',8,... 'Marker','o',... 'LineWidth',2,... 'LineStyle','none',... 'Color',[1 0 0],... 'DisplayName',strcat('Hq(',num2str(min(q)),') = ',num2str(Hq(q==min(q))))); % Create plot plot(qmid1,Hqmid1,'Parent',subplot2,'MarkerFaceColor',[0 0 1],... 'MarkerEdgeColor',[0 0 0],... 'MarkerSize',8,... 'Marker','o',... 'LineWidth',2,... 'LineStyle','none',... 'Color',[0 0 1],... 'DisplayName',strcat('Hq(',num2str(q_middle),') = ',num2str(Hq(q==q_middle)))); % Create plot plot(qstop1,Hqstop1,'Parent',subplot2,'MarkerFaceColor',[0 0.498 0],... 'MarkerEdgeColor',[0 0 0],... 'MarkerSize',8,... 'Marker','o',... 'LineWidth',2,... 'LineStyle','none',... 'Color',[0 0.498 0],... 'DisplayName',strcat('Hq(',num2str(max(q)),') = ',num2str(Hq(q==max(q))))); % Create subplot subplot3 = subplot(2,2,3,'Parent',figure1,... 'XTickLabel',qtick,... 'XTick',qInd,... 'LineWidth',2,... 'FontSize',14); hold(subplot3,'all'); % Create plot plot(X2,Y3,'Parent',subplot3,'Marker','o','LineStyle','none',... 'Color',[0 0 1]); % Create xlabel xlabel('q','FontSize',16); % Create ylabel ylabel('tq','FontSize',16); % Create title title('q-order Mass exponent','FontSize',14); % Create plot plot(qstop1,tqstop1,'Parent',subplot3,'MarkerFaceColor',[0 0.498 0],... 'MarkerEdgeColor',[0 0 0],...