how to make this file work "Error using unique"

조회 수: 11 (최근 30일)
Baozhuo Zhang
Baozhuo Zhang 2015년 2월 9일
댓글: per isakson 2015년 2월 9일
single_peaks=1:8; %single peaks to fit
double_peaks=[];
plot_fits=1; %set to 1 to plot raw and fitted data
azi_bin=360; %number of pre-existing eta bins to put in a single bin for fitting (1 <=azi_bin<= eta.N)
dx=8; %number of radial points on either side of guessed peak position to evaluate for fit
bx=1; %number of background points on either side to use in fitting poly background
fwhmguess=2.5; %guess for peak fwhm, in pixels
maxdrift=3; %maximum range in pixels to find the y max (keep tight as possible!)
peak_tol=0.5; %minimum peak intensity (background subtracted, evaluated over range 2*maxdrift)at which peaks will be fit - otherwise set as 1e-6
below_tol_val=999; %value to give all parameters for peaks which are below intensity tolerance
delta=0; %detector offset parameter, keep at 0 unless you have a good reason not to!
de=0; %strain amplitude (guess)
% crystallographic parameters, hkls to fit and resulting peak position guesses
a1=3.054; % refined from NiTi (austenite_B2)lattice parameters from pdf #18-0899
%a1=3.054; % refined from NiTi (austenite_B2)lattice parameters from pdf #18-0899
% % % NiTiCu B19 lattice parameters from pdf Potatov et. al.
a2=2.919;
b2=4.288;
c2=4.500;
%alpha=1.5708; %90 degrees
%beta=1.5708; %90 degrees
%gamma=1.5708; %90 degrees
% % % NiTiPd B19 lattice parameters from pdf Potatov et. al.
% a3=2.811;
% b3=4.421;
% c3=4.694;
% alpha=1.5708; % 90 degrees
% beta=1.5708; % 90 degrees
% gamma=1.5708; % 90 degrees
% a3=2.885; b3=4.622; c3=4.120; alpha=1.5708; beta=1.68948; gamma=1.5708; % NiTi (mart) B19 lattice parameters from pdf #35-1281
a3=2.898; b3=4.108; c3=4.646; %alpha=1.5708; beta=1.7066; gamma=1.5708; % NiTi (mart) B19' lattice parameters from pdf #03-065-0145
h=[1 0 1 0 1 1 1 1];
k=[0 1 1 0 1 1 2 1];
l=[0 1 0 2 1 2 0 2];
%h=[1 1 0 1 2 2];
%k=[0 1 1 1 0 1];
%l=[0 0 2 1 0 1];
for i=1:1:4; cen0(i)=round(cen0calc('cubic',h(i),k(i),l(i),a2,a2,a2,center.z,energy,pix.size,delta)); end; % NiTi (austenite) cubic peaks
%for i=5:1:11; cen0(i)=round(cen0calcmono('mono',h(i),k(i),l(i),a3,b3,c3,alpha,beta,gamma,center.z,energy,pix.size)); end; % NiTi (martensite) monoclinic peaks
%for i=1:1:6; cen0(i)=round(cen0calc('cubic',h(i),k(i),l(i),a1,a1,a1,center.z,energy,pix.size,delta)); end; % NiTi (austenite) cubic peaks
%for i=2:1:8; cen0(i)=round(cen0calcmono('mono',h(i),k(i),l(i),a3,b3,c3,alpha,beta,gamma,center.z,energy,pix.size)); end; % NiTi (martensite) monoclinic peaks
%%calculate arrays based on user defined parameters
nPeak = length(cen0); nDouble=length(double_peaks);
nAzi=size(imm,1)/azi_bin; %bin data using azi_bin
azi_deg=[0.5:nAzi-0.5]*(eta.max-eta.min)/eta.N*azi_bin; %mean values of azimuths for fit, in degrees
%%fit singlet peaks
opt = optimset('disp','off','large','on','jacobi','on','tolx',1e-4,'tolf',1e-4,'maxi',120);
fn = {'backg1' 'psv1'}; npar = [2 4]; %define peakfitting choice (PV here)
low = [-inf 0 0 0 0 0]; up = [inf inf inf inf inf 1]; %set limits for fitting parameters
pfix = [nan nan nan nan nan nan]; % fix fitting parameters (currently none are fixed)
if ~isempty(single_peaks); fprintf('\nFitting single peaks for following azimuth bins:\n'); end;
for azi = nAzi:-1:1
if ~isempty(single_peaks); fprintf('%4i',azi*azi_bin);
for peak = single_peaks;
cen0a(peak,azi)=round(cen0(peak)*(1+de/2*cos(2.*azi_deg(azi)/57.3))); %azi-dependent starting guess
x = cen0a(peak,azi) - dx: cen0a(peak,azi) + dx; Nx=length(x);
yraw = mean(imm([azi*azi_bin:-1:azi*azi_bin-azi_bin+1], x-rho.min));
xback=x([1:2 Nx-1:Nx]); yback=yraw([1:2 Nx-1:Nx]); %x,y values to use for background determination
[pback_pars,s,mu]=polyfit(xback,yback,1); % 1st order polynomial (linear) background fit
pback=polyval(pback_pars,x,[],mu); %mu are centering and scaling parameters to help polyfit
y=yraw-pback+1000;
ymax= max(max(y)); xmax = x(find(round(y)==round(ymax))); xmax = xmax(1); ymin=min(min(y));
clear x pback_pars pback y yraw xback yback;
x = xmax - dx: xmax + dx; Nx=length(x); % shift data so xmax is center point
x = cen0a(peak,azi) - dx: cen0a(peak,azi) + dx; Nx=length(x); %do not shift data
yraw = mean(imm((azi*azi_bin:-1:azi*azi_bin-azi_bin+1), x-rho.min));ymin=min(min(yraw));
xback=x([1:2 Nx-1:Nx]); yback=yraw([1:2 Nx-1:Nx]); %x,y values to use for background determination
[pback_pars,s,mu]=polyfit(xback,yback,1); % 1st order polynomial (linear) background fit
pback=polyval(pback_pars,x,[],mu); %mu are centering and scaling parameters to help polyfit
y=yraw-pback+1000;
ymax = max(max(y((length(y)-1)/2-maxdrift:(length(y)-1)/2+maxdrift)));
xmax = x(find(round(y)==round(ymax))); xmax = xmax(1); ymin=min(min(y));
%check if peak intensity is above threshold, to decide if this particular azimuth will be fitted
fom=ymax-1000;
if plot_fits
figure(6);clf;
plot(x,yraw,'b+-',x,pback,'r-');
xlabel('radial position, pixels');ylabel('intensity');
title(sprintf('Peak %i and azimuth %i of %s; blue=data, red =background',peak,azi,filename)); grid on;
if (fom>=peak_tol); text(0.05,0.95,sprintf('Will fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause;
else text(0.05,0.95,sprintf('Will NOT fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause; end;
end; %plotting raw data and background
if (fom>=peak_tol); % do the fit as intensity is above threshold
p0 = [0, ymin, (ymax-ymin), xmax, fwhmguess, .5];
y0 = sumfun1(p0,pfix,npar,fn,x);
[pfit,rn,rs,ex,out,lam,jac] = lsqnonlin('sumfun1',p0,low,up,opt,pfix,npar,fn,x,y,ones(size(y)));
[~, var] = confint(pfit,rs,jac); % variance of fitted params
var=full(var);
yfit = sumfun1(pfit,pfix,npar,fn,x);
[~, dum, intint]=psv1([pfit(3:end)],x);
if plot_fits
figure(7);clf; subplot('position',[0.1 0.3 0.85 0.6]);
plot(x,y0,'-r',x,y,'+',x,yfit,'-b'); ylabel('Intensity+100(offset)');
title(sprintf('Fit of peak %i and azimuth %i of %s; symbols=data, blue=fit, red =guess',peak,azi,filename)); grid on;
text(0.05,0.9,sprintf('peak int= %f +- %f',pfit(3),sqrt(var(3))),'sc'); text(0.05,0.85,sprintf('center= %f +- %f',pfit(4),sqrt(var(4))),'sc');
text(0.05,0.8,sprintf('fwhm= %f +- %f',pfit(5),sqrt(var(5))),'sc'); text(0.05,0.75,sprintf('fgauss= %f +- %f',1-pfit(6),sqrt(var(6))),'sc');
text(0.05,0.7,sprintf('integ int= %f ',intint),'sc'); text(0.05,0.65,sprintf('min int(tol)= %f ',peak_tol),'sc');
subplot('position',[0.1 0.1 0.85 0.15]); plot(x,y-yfit'); ylabel('Residual'); xlabel('Radial position'); grid on; pause;
end; %if plot_fits==1;
fit(ifile,peak).int(azi) = pfit(3);fit(ifile,peak).cen(azi) = pfit(4); fit(ifile,peak).fwhm(azi) = pfit(5); fit(ifile,peak).shape(azi)=pfit(6);
fit(ifile,peak).dint(azi) = sqrt(var(3));fit(ifile,peak).dcen(azi) = sqrt(var(4));fit(ifile,peak).dfwhm(azi) = sqrt(var(5));fit(ifile,peak).dshape(azi)=sqrt(var(6));
else %peak intensity too weak to fit
fit(ifile,peak).int(azi) = below_tol_val; fit(ifile,peak).cen(azi) = below_tol_val; fit(ifile,peak).fwhm(azi) = below_tol_val; fit(ifile,peak).shape(azi)=below_tol_val;
fit(ifile,peak).dint(azi) = below_tol_val;fit(ifile,peak).dcen(azi) = below_tol_val;fit(ifile,peak).dfwhm(azi) = below_tol_val;fit(ifile,peak).dshape(azi)=below_tol_val;
end; %if statement for thresholding
end; %looping over azimuths
end; end; %looping over single peaks and checking if there are any singlepeaks
%%fit 2 PV simultaneously
opt = optimset('disp','off','large','on','jacobi','off','tolx',1e-2,'tolf',1e-2,'maxi',80);
slguess=0; fwhmguess1 = 3; fwhmguess2=3; %some guesses for peaks
fn = {'backg1' 'psv1' 'psv1'}; npar = [2 4 4];
low = [-inf 0 0 0 0 0 0 0 0 0]; up = [inf inf inf inf inf 1 inf inf inf 1];
for igroup=1:nDouble;
fprintf('\nFitting two peaks simultaneously for following azimuth bins:\n');
peaks=double_peaks{igroup}; dpeak1=peaks(1);dpeak2=peaks(2);
for azi = nAzi:-1:1
fprintf('%4i',azi*azi_bin);
x = cen0(dpeak2) - dx2: cen0(dpeak1) + dx2; Nx=length(x);
yraw = median(imm([azi*azi_bin:-1:azi*azi_bin-azi_bin+1], x-rho.min));
xback=x([1:2 Nx-1:Nx]); yback=yraw([1:2 Nx-1:Nx]); %x,y values to use for background determination
[pback_pars,s,mu]=polyfit(xback,yback,1); % 1st order polynomial (linear) background fit
pback=polyval(pback_pars,x,[],mu); %mu are centering and scaling parameters to help polyfit
y=yraw-pback+1000;
ymax= max(max(y)); xmax = x(round(y)==round(ymax)); xmax = xmax(1);
ymin = min(y); ymin = ymin(1);
ymax1 = y(x==cen0(dpeak1)); ymax2=y(x==cen0(dpeak2));
fom=ymax-1000;
if plot_fits;
figure(6);clf;
plot(x,yraw,'b+-',x,pback,'r-');
xlabel('radial position, pixels');ylabel('intensity');
title(sprintf('Peaks %i and %i and azimuth %i of %s; blue=data, red =background',dpeak1,dpeak2,azi,filename)); grid on;
if (fom>=peak_tol); text(0.05,0.95,sprintf('Will fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause;
else text(0.05,0.95,sprintf('Will NOT fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause; end;
end; %plotting raw data and background
if (fom>=peak_tol); % do the fit as intensity is above threshold
p0 = [slguess, ymin, (ymax1-ymin), cen0(dpeak1), fwhmguess1, .5,(ymax2-ymin), cen0(dpeak2), fwhmguess2, .5];
y0 = sumfun1(p0,[],npar,fn,x);
[pfit,rn,rs,ex,out,lam,jac] = lsqnonlin('sumfun1',p0,low,up,opt,[],npar,fn,x,y,ones(size(y)));
[dum, var] = confint(pfit,rs,jac); % variance of fitted params
var=full(var);
yfit = sumfun1(pfit,[],npar,fn,x);
yfit1=sumfun1(pfit(1:6),[],[2 4],{'backg1' 'psv1'},x); yfit2=sumfun1(pfit([1:2 7:10]),[],[2 4],{'backg1' 'psv1'},x);
if plot_fits
figure(8);clf; subplot('position',[0.1 0.3 0.85 0.6]);
plot(x,y,'+',x,y0,'-r',x,yfit1,'-b',x,yfit2,'--b',x,yfit,'-k'); ylabel('Intensity+100(offset)');
title(sprintf('Peaks %i and %i & azi %i of %s; symb=data, red=guess, blue/black=fit',dpeak1,dpeak2,azi,filename)); grid on;
text(0.05,0.9,sprintf('peak ints= %f & %f',pfit(3),pfit(7)),'sc'); text(0.05,0.85,sprintf('centers= %f & %f',pfit(4),pfit(8)),'sc');
text(0.05,0.8,sprintf('fwhms= %f & %f',pfit(5),pfit(9)),'sc'); text(0.05,0.75,sprintf('fgauss= %f &%f',1-pfit(6),1-pfit(10)),'sc');
subplot('position',[0.1 0.1 0.85 0.15]); plot(x,y-yfit'); ylabel('Residual'); xlabel('Radial position'); grid on; pause;
end; %if plot_fits==1;
fit(ifile,dpeak1).int(azi) = pfit(3);fit(ifile,dpeak1).cen(azi) = pfit(4); fit(ifile,dpeak1).fwhm(azi) = pfit(5);
fit(ifile,dpeak1).shape(azi) = pfit(6);fit(ifile,dpeak1).dint(azi) = sqrt(var(3));
fit(ifile,dpeak1).dcen(azi) = sqrt(var(4));fit(ifile,dpeak1).dfwhm(azi) = sqrt(var(5));
fit(ifile,dpeak2).int(azi) = pfit(7); fit(ifile,dpeak2).cen(azi) = pfit(8); fit(ifile,dpeak2).fwhm(azi) = pfit(9);
fit(ifile,dpeak2).shape(azi) = pfit(10);fit(ifile,dpeak2).dint(azi) = sqrt(var(7));
fit(ifile,dpeak2).dcen(azi) = sqrt(var(8));fit(ifile,dpeak2).dfwhm(azi) = sqrt(var(9));
else %peak too weak to fit
fit(ifile,dpeak1).int(azi) = below_tol_val; fit(ifile,dpeak1).cen(azi) = below_tol_val; fit(ifile,dpeak1).fwhm(azi) = below_tol_val; fit(ifile,dpeak1).shape(azi)=below_tol_val;
fit(ifile,dpeak1).dint(azi) = below_tol_val;fit(ifile,dpeak1).dcen(azi) = below_tol_val;fit(ifile,dpeak1).dfwhm(azi) = below_tol_val;fit(ifile,dpeak1).dshape(azi)=below_tol_val;
fit(ifile,dpeak2).int(azi) = below_tol_val; fit(ifile,dpeak2).cen(azi) = below_tol_val; fit(ifile,dpeak2).fwhm(azi) = below_tol_val; fit(ifile,dpeak2).shape(azi)=below_tol_val;
fit(ifile,dpeak2).dint(azi) = below_tol_val;fit(ifile,dpeak2).dcen(azi) = below_tol_val;fit(ifile,dpeak2).dfwhm(azi) = below_tol_val;fit(ifile,dpeak2).dshape(azi)=below_tol_val;
end;
end;
end;
fprintf('\n');
  댓글 수: 1
per isakson
per isakson 2015년 2월 9일
  • I don't think unique is called by your script.
  • The script cannot be run because of missing variables, e.g. center
  • Descibe your problem in more detail and upload the file with the "paper-clip-button"

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Descriptive Statistics에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by