大神帮忙分析一下多重分形为程序呗

喜阳阳

有没有用过多重分形为呢，下载了一个程序，看不懂，不会用。大神帮我分析一下呗
clear all
close all;
p1=xlsread('E:matlabdatap1');
p=p1';
C=sum(p);
nn=size(p);
a=length(p);
s=0;
%求kk=[1/a,2/a...1/2,1]
kk=[];
kkk=1/a;
while kkk<=1
     kk=[kk,kkk];
     kkk=kkk*2;
end
for q=-20:4:20
       logyit=[];
       logT=[];
       logXq=[];
       aqfz=[];
       aqfm=[];
%定义了下面两个空集
       fqfz=[];
       fqfm=[];
       for n=1:1:length(kk)
       yit=kk(n);
%计算子集长度T，把3800改为现有数据长度
       T=a*yit;
       P=[];
       X=[];
       for i=1:1:1/yit
          c=p((1+(i-1)*T):(i*T));
          P(i)=sum(c)/C;
          X(i)=[P(i)]^q;
        end
       Xq=sum(X);
       logXq=[logXq,log(Xq)];
       logyit=[logyit,log(yit)];
       logT=[logT,log(T)];
       aqfz=[aqfz sum(X.*log(P))];
       aqfm=[aqfm sum(X)*log(yit)];
       fqfz=[fqfz (q*(sum(X.*log(P)))-sum(X)*log(sum(X)))];
       fqfm=[fqfm sum(X)*log(yit)];
end
s=s+1;
Tq(s)=mean(diff(logXq)./diff(logyit));
aq(s)=mean(diff(aqfz)./diff(aqfm));
fq(s)=mean(diff(fqfz)./diff(fqfm));
plot(logT,logXq,'v-');
hold on
end
figure
plot(aq,fq,'o-');
figure
q=-20:4:20;
%去掉了等号
plot(q,Tq,'x-');

Vickyvictoria

说明一下什么地方看不懂吧

Vickyvictoria

我这儿也有一个相关的程序，算多重分形谱的，你可以看看

1. 
   
2. %%
   
3. % Author: Tegy J. Vadakkan
   
4. % Date: 09/08/2009
   
5. % input a binary image
   
6. % the multifractal spectra is calculated based on the ideas in the paper by
   
7. % Posadas et al., Soil Sci. Soc. Am. J. 67:1361-1369, 2003
   
8. 
   
9. clc;
   
10. clear all;
    
11. close all;
    
12. %%
    
13. syms s;
    
14. indata=inputdlg({'input photo'});
    
15. a = imread(indata{1});
    
16. [rows, cols] = size(a);
    
17. figure;imshow(a);
    
18. npix = sum(sum(a));
    
19. %% calculates niL which is the number of pixels in the ith box of size L
    
20. % ideas from boxcount.m by F. Moisy have been borrowed here
    
21. width = rows;
    
22. p = log(width)/log(2);
    
23. max_boxes = power(rows,2)/power(2,2);
    
24. nL = double(zeros(max_boxes,p));
    
25. for g=(p-1):-1:0
    
26.             siz = 2^(p-g);
    
27.             sizm1 = siz - 1;
    
28.             index = log2(siz);
    
29.             count = 0;
    
30.             for i=1:siz:(width-siz+1)
    
31.                 for j=1:siz:(width-siz+1)
    
32.                     count = count + 1;
    
33.                     sums = sum(sum(a(i:i+sizm1,j:j+sizm1)));
    
34.                     nL(count,index) = sums;
    
35.                 end
    
36.             end
    
37. end
    
38. %%
    
39. qran = 1;
    
40. logl = zeros(p,1);
    
41. 
    
42. for l=1:p
    
43.     logl(l) = log(power(2,l));
    
44. end
    
45. %% normalized masses
    
46. pL = double(zeros(max_boxes,p));
    
47. for l=1:p
    
48.     nboxes = power(rows,2)/power(power(2,l),2);
    
49.     norm = sum(nL(1:nboxes,l));
    
50.     if(norm ~= npix)
    
51.         display('error');
    
52.     end
    
53.     for i=1:nboxes
    
54.         pL(i,l) = nL(i,l)/norm;
    
55.     end
    
56. end
    
57. %%
    
58. %falpha, alpha
    
59. for l=1:p
    
60.    
    
61.     count = 0;
    
62.     nboxes = power(rows,2)/power(power(2,l),2);
    
63.    
    
64.     for q = -qran:+0.1:qran      
    
65.         
    
66.         %denominator of muiql
    
67.         qsum = 0.0;
    
68.         for i=1:nboxes
    
69.             if(pL(i,l) ~= 0)
    
70.                 qsum = qsum + power(pL(i,l),q);
    
71.             end
    
72.         end
    
73. 
    
74.         fqnum = 0.0;
    
75.         aqnum = 0.0;
    
76.         smuiqL = 0.0;
    
77.         for i=1:nboxes
    
78.             if(pL(i,l) ~= 0)
    
79.                   muiqL = power(pL(i,l),q)/qsum;
    
80.                   fqnum = fqnum + (muiqL * log(muiqL));
    
81.                   aqnum = aqnum + (muiqL * log(pL(i,l)));
    
82.                   smuiqL = smuiqL + muiqL;
    
83.             end
    
84.         end
    
85.         if(uint8(smuiqL)~=1)
    
86.             display('error');
    
87.         end
    
88.         
    
89.         count = count + 1;
    
90.         fql(l,count) = fqnum;
    
91.         aql(l,count) = aqnum;
    
92.         qval(count) = q;
    
93.     end
    
94. 
    
95. end
    
96. %%
    
97. % alpha_q
    
98. for i=1:count
    
99.      line = polyfit(logl,aql(:,i),1);
    
100.      aq(i) = line(1);
     
101.      yfit = polyval(line,logl);
     
102.      sse = sum(power(aql(:,i)-yfit,2));
     
103.      sst = sum(power(aql(:,i)-mean(aql(:,i)),2));
     
104.      ar2(i) = 1-(sse/sst);
     
105. end
     
106. % f_q
     
107. for i=1:count
     
108.      line = polyfit(logl,fql(:,i),1);
     
109.      fq(i) = line(1);
     
110.      yfit = polyval(line,logl);
     
111.      sse = sum(power(fql(:,i)-yfit,2));
     
112.      sst = sum(power(fql(:,i)-mean(fql(:,i)),2));
     
113.      fr2(i) = 1-(sse/sst);
     
114. end
     
115. figure;plot(qval,aq,'r:o',qval,fq,'g:o');grid on;
     
116. h = legend('alpha(q)','f(q)',1);
     
117. xlabel('q','FontSize',14);
     
118. %gird on;
     
119. 
     
120. figure;plot(aq,fq,'r:o');grid on;
     
121. xlabel('alpha(q)','FontSize',14);
     
122. ylabel('f(q)','FontSize',14);
     
123. %%T=fmaxbnd(aq,2.00,4.00);
     
124. %xt=solve(aq);
     
125. %grid on;
     
126. dy=diff(aq);
     
127. xz=solve(dy);
     
128. 
     
129. 
     
130. 
     
131. line=polyfit(aq,fq,2);
     
132. pfit = polyval(line,aq);
     
133. figure;plot(aq,fq,'r:o',aq,pfit,'g:o');grid on;
     
134. h = legend('f(q)','Parabolic fit to f(q)',3);
     
135. xlabel('alpha(q)','FontSize',14);
     
136. %grid on;

复制代码

喜阳阳

Vickyvictoria 发表于 2015-12-16 09:42
我这儿也有一个相关的程序，算多重分形谱的，你可以看看

谢谢啦