clear all
close all
clc
a=load('F:\硕士毕业论文\各城市转加权矩阵\三个地区原数据均值.txt');
%y=a(:,1);
pvc1=a(:,1);
%pvc1=dlmread('cu-01.txt','\t',[5000,1,8000,1]); %100(male 69) 正常
X = [pvc1];
x = 1:3001
plot(x,X);
q = [-30:1:-8,-7.2:0.9:-5.4,-5.3:0.4:-2.9,-2.76:0.08:-0.2,-0.1:0.2:1,1.6:0.8:7.2,8:1:30];
% q = [-10:2:10];
% averageX=mean(X);
% N=length(X);
S = 4:10:500;
% for i=1:N %%%%%%%%%%%%%%%%MF-DFA 算法的第一步,构造去均值的和序列
%for i2 = 1:size(X,2); %%size(X,2)返回X第二维的长度
N=length(X);
averageX=sum(X)./N;
Y = [];
for a=1:N
temp=0;
for k=1:a
temp=temp+(X(k)-averageX);
end
Y(a)=temp;
end
%%%%%%%%%%%%%%%%%第二次求离差
averageY = mean(Y);
Y_1 = [];
for a_1 = 1:N
temp_1 = 0;
for k_1 = 1:a_1
temp_1 = temp_1+(Y(k_1)-averageY);
end
Y_1(a_1)=temp_1;
end
M=length(Y_1);
for q_i = 1:length(q)
for i = 1:length(S)
N1=floor(M/S(i));
y1 = reshape(Y_1(1:S(i)*N1),S(i),N1); %以下实现从前到后,从后到前划分序列,得S行,2*N1列
y2 = reshape(Y_1(M:-1:M-S(i)*N1+1),S(i),N1);
y = [y1,y2];
N_totle = 2*N1;
coefPoly = zeros(2,N_totle); %设置空数组为拟合做准备
valueFit = zeros(size(y));
for segment_ii = 1:N_totle %对划分后的每一个区间进行拟合
coefPoly(:,segment_ii) = polyfit([1:S(i)]',y(:,segment_ii),1); %分别去每一段进行拟合
valueFit(:,segment_ii) = polyval(coefPoly(:,segment_ii),1:S(i));
end
diffMatrix = valueFit - y; %都是矩阵形式,得到多个数据
diffMMSE = mean(diffMatrix.^2);
FuncQs(i) = (mean(diffMMSE.^(q(q_i)/2)))^(1/q(q_i)); %获得q阶波动函数Fq(s)
end %figure(1),plot(log10(S),log10(FuncQs),'-ok','MarkerSize',8), xlabel('lgS'),ylabel('lgFq(S)');hold on
coef = polyfit(log10(S),log10(FuncQs),1);
valuepoly = polyval(coef,log10(S));
%figure(2),plot(log10(S),valuepoly,'-ok','MarkerSize',4), xlabel('lgS'),ylabel('valuepoly');hold on
hq(q_i) = coef(1);%figure(1),plot(q(q_i),hq(q_i),'-ok','MarkerSize',4), xlabel('q'),ylabel('hq');hold on
end
figure(1),plot(q,hq,'-ok','MarkerSize',4), xlabel('q'),ylabel('hq');hold on
%%%%%%%%%%%%%%%%%%%%求T(q)-q的关系图
T = [];
for q_m = 1:length(q)
T(q_m) = q(q_m).*hq(q_m)-1;
end
%figure(2),plot(q,T,'-ok','MarkerSize',4), xlabel('q'),ylabel('T(q)');hold on
figure(2),plot(q,T,'-ok','MarkerSize',4), xlabel('q'),ylabel('T(q)');hold on
%%%%%%%%%%%%%%%%%%%%%%%%%求多重分形谱
Dhq = diff(hq)./diff(q);
alpha = [];f_alpha = [];
for i = 1:length(Dhq)
alpha(i) = hq(i) + q(i).* Dhq(i);
f_alpha(i) = (q(i).^2).* Dhq(i) + 1;
% F_alpha = abs(f_alpha);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%去掉f(a)中的负值
Alpha=[];F_alpha=[];
j=1;
for i=1:length(f_alpha)
if f_alpha(i)>0
F_alpha(j) = f_alpha(i);
Alpha(j)=alpha(i);
j=j+1;
end
end
%figure(3),plot(Alpha,F_alpha,'-ok','MarkerSize',4), xlabel('Alpha'),ylabel('F_alpha');hold on
figure(3),plot(Alpha,F_alpha,'-ok','MarkerSize',4),xlabel('Alpha'),ylabel('F_alpha');hold on
% DFA=[];s=1;
% DFA(s)=max(Alpha)-min(Alpha);s=s+1;
% DFA(s)=max(F);
评论1