计算分形维数的matlab实现程序出现的问题

mumianhua

备注： 我是对May模型进行分析的，运行过程中，总是提醒有警告，结果中说“函数拟合polyfit”有问题，我计算的分维数，和书上计算的不太一致，我想请各位大师给予我指导和帮助！  下面我贴上我的程序，希望大家给我指导哦  非常感谢，也欢迎大家交流

mumianhua

1. clear
   
2. mu=0:0.1:4;
   
3. y=0.001*ones(2001,1);
   
4. yn=[];
   
5. for i=1:length(mu)
   
6.     for n=1:2000
   
7.         y(n+1)=mu(i)*y(n)*(1-y(n));
   
8.     end
   
9. yn=[yn,y(1500:2000)];   
   
10. i   
    
11. end
    
12. plot(mu,yn,'r*')
    
13. hold on
    
14. 
    
15. cellmax=2^9;
    
16. for i=1:length(mu)
    
17.     D(i)=FractalDim(yn(:,i),cellmax);
    
18. end
    
19. plot(mu,D)
    
20. %%%%%%%%%%%%%%%%%%%%%%%
    
21. %%%%%%%%%%%%%%%%%%%%%%%
    
22. function D=FractalDim(y,cellmax)
    
23. % y=load('C:\Users\Administrator\Desktop\分形维数-及其案例程序（关联维数）\[w=1.38].dat');%%%%%%%%%导入的数据
    
24. %求输入一维信号的计盒分形维数
    
25. %y是一维信号
    
26. %cellmax:方格子的最大边长,可以取2的偶数次幂次(1,2,4,8...),取大于数据长度的偶数
    
27. %D是y的计盒维数（一般情况下D>=1）,D=lim(log(N(e))/log(k/e)),
    
28. 
    
29. 
    
30. if cellmax<length(y)
    
31. error('cellmax must be larger than input signal!')
    
32. end
    
33. L=length(y);%输入样点的个数
    
34. y_min=min(y);
    
35. 
    
36. %移位操作，将y_min移到坐标0点
    
37. y_shift=y-y_min;
    
38. %重采样，使总点数等于cellmax+1
    
39. x_ord=[0:L-1]./(L-1);
    
40. xx_ord=[0:cellmax]./(cellmax);
    
41. y_interp=interp1(x_ord,y_shift,xx_ord);
    
42. %按比例缩放y，使最大值为2^^c
    
43. ys_max=max(y_interp);
    
44. factory=cellmax/ys_max;
    
45. yy=abs(y_interp*factory);
    
46. 
    
47. t=log2(cellmax)+1;%叠代次数
    
48. for e=1:t
    
49. Ne=0;%累积覆盖信号的格子的总数
    
50. cellsize=2^(e-1);%每次的格子大小
    
51. NumSeg(e)=cellmax/cellsize;%横轴划分成的段数
    
52. 
    
53. for j=1:NumSeg(e) %由横轴第一个段起通过计算纵轴跨越的格子数累积N(e)
    
54. begin=cellsize*(j-1)+1;%每一段的起始
    
55. tail=cellsize*j+1;
    
56. seg=[begin:tail];%段坐标
    
57. yy_max=max(yy(seg));
    
58. yy_min=min(yy(seg));
    
59. up=ceil(yy_max/cellsize);
    
60. down=floor(yy_min/cellsize);
    
61. Ns=up-down;% 本段曲线占有的格子数
    
62. Ne=Ne+Ns;%累加每一段覆盖曲线的格子数
    
63. 
    
64. end
    
65. 
    
66. N(e)=Ne;%记录每e下的N(e)
    
67. end
    
68. 
    
69. %对log(N(e))和log(k/e)进行最小二乘的一次曲线拟合,斜率就是D
    
70. 
    
71. r=-diff(log2(N));%去掉r超过2和小于1的野点数据
    
72. id=find(r<=2&r>=1);%保留的数据点
    
73. Ne=N(id);
    
74. e=NumSeg(id);
    
75. P=polyfit(log2(e),log2(Ne),1);%一次曲线拟合返回斜率和截距
    
76. D=P(1);
    
77. %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
78. %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
79. function [prr,pcr,p]=glws(x,m,t)
    
80. %函数名为关联维数的首字母,用于单串序列，多串到glsw；
    
81. %x为要分析的数据；
    
82. %x=xlsread('d:\matworks\dbin.xls');
    
83. [m1,n1]=size(x);
    
84. n=m1;
    
85. [mm1,mm]=size(m);
    
86. p=zeros(mm,2); %存放拟合系数的矩阵；
    
87. rr=zeros(20,mm);%rr是相当于筛子的那个距离，存放的是对数；
    
88. cr=zeros(20,mm);%cr是小于筛子距离的距离个数，存放的是对数；
    
89. %prr=zeros(20,mm);%rr是相当于筛子的那个距离，存放的是对数；
    
90. %pcr=zeros(20,mm);%cr是小于筛子距离的距离个数，存放的是对数；
    
91. scope=zeros(19,1);
    
92. msr=zeros(19,1);
    
93.    for k=1:mm
    
94.       tt=0;
    
95.       nm=n-(m(k)-1)*t;%Nm为列数；
    
96.       nr=(nm-1)*nm/2;%Nr为距离的总个数；
    
97.       juli=zeros(nr,1);%全部距离搞成一列的长矩阵；
    
98.       r=zeros(nm,nm);%各列之间距离矩阵；
    
99.       y=zeros(m(k),nm);%重构相矩阵的值yij；
    
100.      for j=1:nm
     
101.          for i=1:m(k)
     
102.             y(i,j)=x(j+(i-1)*t);
     
103.          end
     
104.      end
     
105.      for i=1:nm-1
     
106.          for j=i+1:nm
     
107.              for kk=1:m(k)
     
108.                r(i,j)=r(i,j)+(y(kk,j)-y(kk,i))^2;
     
109.              end
     
110.              r(i,j)=sqrt(r(i,j));
     
111.              tt=tt+1;
     
112.              juli(tt)=r(i,j);
     
113.          end
     
114.      end
     
115.   %进行r和cr个数的计算；
     
116.     rmin=min(juli);
     
117.     rmax=max(juli);
     
118.     for i=1:20  %每次把距离间隔分20分来慢慢加；
     
119.        rr(i,k)=(rmax-rmin)*(i+1)/21; %距离取法值得研究一下；
     
120.        for j=1:nr
     
121.            if juli(j)<=rr(i,k)
     
122.                 cr(i,k)=cr(i,k)+1;
     
123.            end
     
124.        end
     
125.        rr(i,k)=log(rr(i,k));
     
126.        cr(i,k)=log(cr(i,k)/nr);
     
127.      end
     
128.     %rr=rr';
     
129.     tt=0;
     
130.     for i=1:19
     
131.         scope(i)=(cr(i+1,k)-cr(i,k))/(rr(i+1,k)-rr(i,k));%每点的斜率；
     
132.         tt=tt+scope(i);
     
133.         plot(i,scope(i),'-bd'),hold on;
     
134.     end
     
135.     tt=tt/19;%各相邻点间斜率平均值；
     
136.     tshold=(max(scope)-min(scope))/2;%threshold，阈值;
     
137.     for i=1:19
     
138.           msr(i)=abs(scope(i)-tt);  %各斜率与平均值的均方根，mean square root;
     
139.     end
     
140.     tt=0;
     
141.     for i=2:18
     
142.           if (msr(i-1)>tshold && msr(i+1)>tshold)||(msr(i-1)<0.001 && msr(i+1)<0.001)
     
143.               continue
     
144.           else
     
145.               tt=tt+1;
     
146.               prr(tt)=rr(i,k);%符合条件的；
     
147.               pcr(tt)=cr(i,k);
     
148.           end
     
149.     end
     
150.      p(k,1:2)=polyfit(prr,pcr,1);%线性拟合，p为两个数，p1为斜率，p2为截距；
     
151.   end

复制代码

mumianhua

下面是我运行的结果和书上计算的结果对比如下图

mumianhua

关键问题是拟合的曲线为什么那么别扭呢，有其他方法吗？   

Derecksky

路过         看看