for i=1:96
for j=1:96
xdata=[4 8 12 16 20 24 28 32];
if ratio1(i,j)*ratio2(i,j)*ratio3(i,j)*ratio4(i,j)*ratio5(i,j)*ratio6(i,j)*ratio7(i,j)*ratio8(i,j)==0;
DT(i,j)=0;
DMean(i,j)=0;
else
ydata(1)=ratio1(i,j);ydata(2)=ratio2(i,j); ydata(3)=ratio3(i,j);ydata(4)=ratio4(i,j);
ydata(5)=ratio5(i,j);ydata(6)=ratio6(i,j);ydata(7)=ratio7(i,j);ydata(8)=ratio8(i,j);
x0=[0.001 0.01]';
xmin=[0.0005 0];
xmax=[0.1 0.5];
x=lsqcurvefit(@MultB,x0,xdata,ydata,xmin,xmax);
%运用3sigema法则去掉偏差较大的点
sigema=0;
Vars=zeros(1,length(ydata));
Vars=1;
while(any(abs(Vars)>1.5*sigema))
ydataLsq=MultB(x,xdata);
VarMean=sum(abs(ydata-ydataLsq))/length(ydata);
Vars=abs(ydata-ydataLsq)-VarMean;
sigema=sqrt(sum(Vars.^2)/(length(ydata)-1));
if(any(abs(Vars)>1.5*sigema))
n=find(max(abs(Vars)));
ydata(n)=[];
xdata(n)=[];
end
x=lsqcurvefit(@MultB,x0,xdata,ydata,xmin,xmax)
end
%% MultB函数
function f=MultB(x,xdata)
n=length(xdata);
for i=1:n
f(i)=exp(-xdata(i)*x(1))*sqrt(pi/(4*xdata(i)*x(2)));
end