|
给提一个多重分形谱算法的参考小程 。程序如下:- A=imread('cameraman.tif');
- L=length(A);
- i=1;
- modify=1;
- tmin=2; % 边框间距,“※”
- tmax=10;
- ttmin=-10;
- ttmax=10; % 自定义 q 的范围
- for r=tmin:1:tmax
- c(i,1)=mod(L,r);
- i=i+1;
- end
- c'; % 计算不能被边长r整除的余数
- a=L-c'; % 计算并剔除掉不能被边长r整除的原始数据
- n=length(a); % 求解格网化边长的个数,即为 n
- TT=[];
- j=1;
- r=tmin; % 自定义项,“※-2”
- for i=1:1:n % 即n=25-10+1,自定义的结果
- B=A(1:a(i),1);
- U=reshape(B,r,length(B)/r);
- T=mean(U);
- T=T'.*r^3;
- TT(1: length(T),i)=[T];
- modifying(modify,1)=length(T);
- r=r+1;
- if r>= tmin+n % 自定义项,“※-3”
- break; % 边长超过10+n,超过初始限制,则程序自动终止
- end
- modify=modify+1;
- end
- modifying;
- TT= nthroot(TT,1);
- % 或者不缩小
- % TT;
- TT=nonzeros(TT);
- for cugb=1:1:n
- modifying_modifying(cugb,1)=sum(modifying(1:cugb)); % 有影响的新加卷 % & * % ¥ # @ !) ……
- end
- modifying_modifying;
- j=1;
- % q 为任意数,这里取1到n,为n,与 k取值保持一致,q过大,计算机无法
- %识别,默认为无穷大,q过小,结果接近0,则意义不明确
- for q=ttmin:ttmax %这里取 q=-10:1:10
- for k=1:1:n
- X=TT(1:modifying_modifying(k,1),1).^q;
- if k>1
- X=TT(modifying_modifying(k-1,1)+1:modifying_modifying(k,1),1).^q;
- end
- t=sum(X);
- XX(k,j)=[t]; % 这里用到两个循环,即考虑到了幂函数,又需考虑求和
- end
- j=j+1;
- end
- XX; % 得到质量分配函数,Xq(ξ),
- X=log(tmin:1:tmax);
- % X=log(tmax:-1:tmin); % 此系以前的自定义输入结果,“※—4”
- Y=log(XX);
- % figure(1)
- % plot(X',Y,'o-k') % 至此,计算多重分形谱的第一步,分配函数构建完毕
- side_length= tmin:1:tmax; % 自定义网格边长,“※—5”
- side_length=side_length';
- q=ttmin:ttmax; % q=-5:1:n-10,q=-10:1:10
- m=1;
- [ha,hb]=size(XX);
- for i=1:1:hb
- % XX=XX’; % or not
- s=XX(:,i); % XX
- b=polyfit(log(side_length),log(s),1); % 在对数尺度下计算斜率
- slope(m,1)=b(1,1);
- m=m+1;
- end
- slope; % 这里的Slope即为质量指数,τ(q)
- N=polyfit(q', slope,1);
- plot(q', slope) % 此步是考察τ(q)-q 之间的关系,
- a=diff(slope)./diff(q'); % 第三步计算,diff函数求偏导确实少一列
- q=q';
- f_a=a.*q(1:end-1)-slope(1:end-1);
- % figure(2)
- % a=sort(a,'ascend');
- plot(a,f_a,'o-k')
- xlabel('α','FontSize',12);
- ylabel('f(α)','FontSize',12);
- % polyfit(a,f_a,3)
- a=sort(a,'ascend');
- % a+1
- % f_a+1
- a % 奇异性指数
复制代码 另外还有一个是盒子维数的计算的程序 (备注:运行的时候,把文件名改一下,不要出现中文)以及多重分形谱及其计算的小论文(见附件)。
|
|