计算分形维数的matlab实现程序出现的问题
备注： 我是对May模型进行分析的，运行过程中，总是提醒有警告，结果中说“函数拟合polyfit”有问题，我计算的分维数，和书上计算的不太一致，我想请各位大师给予我指导和帮助！下面我贴上我的程序，希望大家给我指导哦非常感谢，也欢迎大家交流
本帖最后由 牛小贱 于 2015-3-15 15:21 编辑
clear
mu=0:0.1:4;
y=0.001*ones(2001,1);
yn=[];
for i=1:length(mu)
for n=1:2000
y(n+1)=mu(i)*y(n)*(1-y(n));
end
yn=;
i
end
plot(mu,yn,'r*')
hold on
cellmax=2^9;
for i=1:length(mu)
D(i)=FractalDim(yn(:,i),cellmax);
end
plot(mu,D)
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
function D=FractalDim(y,cellmax)
% y=load('C:\Users\Administrator\Desktop\分形维数-及其案例程序（关联维数）\.dat');%%%%%%%%%导入的数据
%求输入一维信号的计盒分形维数
%y是一维信号
%cellmax:方格子的最大边长,可以取2的偶数次幂次(1,2,4,8...),取大于数据长度的偶数
%D是y的计盒维数（一般情况下D>=1）,D=lim(log(N(e))/log(k/e)),
if cellmax<length(y)
error('cellmax must be larger than input signal!')
end
L=length(y);%输入样点的个数
y_min=min(y);
%移位操作，将y_min移到坐标0点
y_shift=y-y_min;
%重采样，使总点数等于cellmax+1
x_ord=./(L-1);
xx_ord=./(cellmax);
y_interp=interp1(x_ord,y_shift,xx_ord);
%按比例缩放y，使最大值为2^^c
ys_max=max(y_interp);
factory=cellmax/ys_max;
yy=abs(y_interp*factory);
t=log2(cellmax)+1;%叠代次数
for e=1:t
Ne=0;%累积覆盖信号的格子的总数
cellsize=2^(e-1);%每次的格子大小
NumSeg(e)=cellmax/cellsize;%横轴划分成的段数
for j=1:NumSeg(e) %由横轴第一个段起通过计算纵轴跨越的格子数累积N(e)
begin=cellsize*(j-1)+1;%每一段的起始
tail=cellsize*j+1;
seg=;%段坐标
yy_max=max(yy(seg));
yy_min=min(yy(seg));
up=ceil(yy_max/cellsize);
down=floor(yy_min/cellsize);
Ns=up-down;% 本段曲线占有的格子数
Ne=Ne+Ns;%累加每一段覆盖曲线的格子数
end
N(e)=Ne;%记录每e下的N(e)
end
%对log(N(e))和log(k/e)进行最小二乘的一次曲线拟合,斜率就是D
r=-diff(log2(N));%去掉r超过2和小于1的野点数据
id=find(r<=2&r>=1);%保留的数据点
Ne=N(id);
e=NumSeg(id);
P=polyfit(log2(e),log2(Ne),1);%一次曲线拟合返回斜率和截距
D=P(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function =glws(x,m,t)
%函数名为关联维数的首字母,用于单串序列，多串到glsw；
%x为要分析的数据；
%x=xlsread('d:\matworks\dbin.xls');
=size(x);
n=m1;
=size(m);
p=zeros(mm,2); %存放拟合系数的矩阵；
rr=zeros(20,mm);%rr是相当于筛子的那个距离，存放的是对数；
cr=zeros(20,mm);%cr是小于筛子距离的距离个数，存放的是对数；
%prr=zeros(20,mm);%rr是相当于筛子的那个距离，存放的是对数；
%pcr=zeros(20,mm);%cr是小于筛子距离的距离个数，存放的是对数；
scope=zeros(19,1);
msr=zeros(19,1);
for k=1:mm
tt=0;
nm=n-(m(k)-1)*t;%Nm为列数；
nr=(nm-1)*nm/2;%Nr为距离的总个数；
juli=zeros(nr,1);%全部距离搞成一列的长矩阵；
r=zeros(nm,nm);%各列之间距离矩阵；
y=zeros(m(k),nm);%重构相矩阵的值yij；
for j=1:nm
for i=1:m(k)
y(i,j)=x(j+(i-1)*t);
end
end
for i=1:nm-1
for j=i+1:nm
for kk=1:m(k)
r(i,j)=r(i,j)+(y(kk,j)-y(kk,i))^2;
end
r(i,j)=sqrt(r(i,j));
tt=tt+1;
juli(tt)=r(i,j);
end
end
%进行r和cr个数的计算；
rmin=min(juli);
rmax=max(juli);
for i=1:20%每次把距离间隔分20分来慢慢加；
rr(i,k)=(rmax-rmin)*(i+1)/21; %距离取法值得研究一下；
for j=1:nr
if juli(j)<=rr(i,k)
cr(i,k)=cr(i,k)+1;
end
end
rr(i,k)=log(rr(i,k));
cr(i,k)=log(cr(i,k)/nr);
end
%rr=rr';
tt=0;
for i=1:19
scope(i)=(cr(i+1,k)-cr(i,k))/(rr(i+1,k)-rr(i,k));%每点的斜率；
tt=tt+scope(i);
plot(i,scope(i),'-bd'),hold on;
end
tt=tt/19;%各相邻点间斜率平均值；
tshold=(max(scope)-min(scope))/2;%threshold，阈值;
for i=1:19
msr(i)=abs(scope(i)-tt);%各斜率与平均值的均方根，mean square root;
end
tt=0;
for i=2:18
if (msr(i-1)>tshold && msr(i+1)>tshold)||(msr(i-1)<0.001 && msr(i+1)<0.001)
continue
else
tt=tt+1;
prr(tt)=rr(i,k);%符合条件的；
pcr(tt)=cr(i,k);
end
end
p(k,1:2)=polyfit(prr,pcr,1);%线性拟合，p为两个数，p1为斜率，p2为截距；
end
下面是我运行的结果和书上计算的结果对比如下图
关键问题是拟合的曲线为什么那么别扭呢，有其他方法吗？ 路过 看看
页:
[1]