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)&gt1.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)&gt1.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