%根據計盒維數原理編寫了求一維曲線分形維數的matlab程序 
 
function D=FractalDim(y,cellmax)
%求輸入一維信號的計盒分形維數
%y是一維信號
%cellmax:方格子的最大邊長,可以取2的偶數次冪次(1,2,4,8...),取大于數據長度的偶數
%D是y的計盒維數（一般情況下D>=1）,D=lim(log(N(e))/log(k/e)),
clear;
y=load('1d200.dat');
%y=[1,1,5,1,1,5];
cellmax=20480;
plot(y);
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=[0:L-1]./(L-1);
xx_ord=[0:cellmax]./(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=[begin:tail];%段坐標
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);