fitgaussian

 GAUSSIAN_PARAM FITERR GAUSSFIT SIGERROR]= FITGAUSSIAN(DATA,[property_value_pair]);
 
 DATA is a 1D vector to which we want to fit a gaussian profile. The
 function will return a structure with various fitted parameters.
 SCALEFACTOR is optional and is applied in the horizontal axis.

 Property value pairs:

  scale    : 1 (default)
  plot     : 0 no plots (default), 1 plots fits against data

 Initial Fit parameters are fit automatically unless specified.
  integral : estimated sum integral of the function
  mean     : estimate of the mean of the gaussian 
  sigma    : estimated sigma of the gaussian
  DC       : known DC component to remove (set to zero to not fit)
  bg_grad  : known background gradient to fit (set to zero to not fit)
  Assym    : known gaussian assymetry (set to zero to not fit)

 FITERR    Final fit error returned by the internal error fuction 
 GAUSSFIT  Final fit function
 SIGERROR  Error estimate of the sigma fit in relative terms. Multiply by
           100 to get in percent.

 Original script by R. Dowd
 In functional format by E. Tan 31/07/2009

0001 function varargout = fitgaussian(varargin)
0002 
0003 % GAUSSIAN_PARAM FITERR GAUSSFIT SIGERROR]= FITGAUSSIAN(DATA,[property_value_pair]);
0004 %
0005 % DATA is a 1D vector to which we want to fit a gaussian profile. The
0006 % function will return a structure with various fitted parameters.
0007 % SCALEFACTOR is optional and is applied in the horizontal axis.
0008 %
0009 % Property value pairs:
0010 %
0011 %  scale    : 1 (default)
0012 %  plot     : 0 no plots (default), 1 plots fits against data
0013 %
0014 % Initial Fit parameters are fit automatically unless specified.
0015 %  integral : estimated sum integral of the function
0016 %  mean     : estimate of the mean of the gaussian
0017 %  sigma    : estimated sigma of the gaussian
0018 %  DC       : known DC component to remove (set to zero to not fit)
0019 %  bg_grad  : known background gradient to fit (set to zero to not fit)
0020 %  Assym    : known gaussian assymetry (set to zero to not fit)
0021 %
0022 % FITERR    Final fit error returned by the internal error fuction
0023 % GAUSSFIT  Final fit function
0024 % SIGERROR  Error estimate of the sigma fit in relative terms. Multiply by
0025 %           100 to get in percent.
0026 %
0027 % Original script by R. Dowd
0028 % In functional format by E. Tan 31/07/2009
0029 
0030 DEBUG = 0;
0031 ASYM_FLAG = 0;
0032 
0033 [reg, prop] = parseparams(varargin);
0034 
0035 if nargin > 0
0036     data = reg{1};
0037 else
0038     disp('No data to work with');
0039     return
0040 end
0041 data = data(:)';
0042 
0043 % defaults
0044 scalefac = 1;
0045 plotfig = 0;
0046 autofitDC = 1;
0047 autofitGrad = 1;
0048 autofitAssym = 1;
0049 userx = [];
0050 for i=1:length(prop)/2
0051     ind = (i-1)*2+1;
0052     switch lower(prop{ind})
0053         case 'scale'
0054             scalefac = prop{ind+1};
0055         case 'plot'
0056             plotfig = prop{ind+1};
0057         case 'integral'
0058             guess_area = prop{ind+1};
0059         case 'mean'
0060             guess_mu = prop{ind+1};
0061         case 'sigma'
0062             guess_sigma = prop{ind+1};
0063         case 'dc'
0064             DCfit = prop{ind+1};
0065             autofitDC = 0;
0066         case 'bg_grad'
0067             linefit = prop{ind+1};
0068             autofitGrad = 0;
0069         case 'assym'
0070             Assym = prop{ind+1};
0071             autofitAssym = 0;
0072         case 'x'
0073             userx = prop{ind+1};
0074     end
0075 end
0076 
0077 % percentage of at the start and end of the data set assumed to be
0078 % representative of the background noise.
0079 pcbackground = 0.1;
0080 
0081 datasize = length(data);
0082 
0083 if isempty(userx)
0084     x=[1:datasize]*scalefac;
0085 else
0086     x = userx;
0087 end
0088 
0089 % Fit a gaussian; first find some starting parameters
0090 
0091 % Used to guess the DC component. Assume flat for the first 10% of data
0092 % points.
0093 startDC = mean(data(1:fix(datasize*pcbackground)));
0094 endDC =  mean(data(end-fix(datasize*pcbackground):end));
0095 
0096 
0097 if ~exist('guess_area','var')
0098     % Guess area AUTOMATICALLY
0099 %     guess_area = sum(data) - 0.5*(startDC+endDC)*datasize;
0100 %     guess_area = guess_area*scalefac;
0101     guess_area = sum((data(2:end) + data(1:end-1)).*diff(x)/2) - 0.5*(startDC+endDC)*(x(end)-x(1));
0102     guess_area = guess_area;
0103 end
0104 
0105 if ~exist('guess_mu','var')
0106     % Guess the center of mass AUTOMATICALLY
0107     [~, maxind] = max(data);
0108     guess_mu = x(maxind);
0109 end
0110 
0111 if ~exist('guess_sigma','var')
0112     % Guess sigma in pixels AUTOMATICALLY
0113     maxval = max(data);
0114     indices = find(data > (maxval+((startDC+endDC)/2) )/2);
0115     guess_sigma = (x(indices(end)) - x(indices(1)))/2.3;
0116     guess_sigma = guess_sigma;
0117 end
0118 
0119 % So far everything has been calculated in units of data points. Apply
0120 % scaling factor here.
0121 fixedvals = [NaN NaN NaN];
0122 
0123 Starting(1) = guess_area;
0124 Starting(2) = guess_mu;
0125 Starting(3) = guess_sigma;
0126 if autofitDC
0127     Starting(end+1) = startDC;
0128 else
0129     fixedvals(1) = DCfit;
0130 end
0131 if autofitGrad
0132     % Guess if there is a background gradient. Again assuming first and last
0133     % 10% of data set is "background".
0134     Starting(end+1) = -(endDC-startDC)/datasize*scalefac;
0135 else
0136     fixedvals(2) = linefit;
0137 end
0138 if autofitAssym
0139     % Initial Assymetry factor
0140     Starting(end+1) = 0;
0141 else
0142     fixedvals(3) = Assym;
0143 end
0144 
0145 
0146 if DEBUG
0147     options = optimset('Display','iter','MaxIter',1500,'TolX',1e-6,'TolFun',1e-10);
0148 else        
0149     options = optimset('Display','off','MaxIter',1500,'TolX',1e-6,'TolFun',1e-10);
0150 end
0151 
0152 [Estimates fval] = fminsearch(@myfit,Starting,options,x,data,fixedvals);
0153 
0154 
0155 fitparam.xdata = x;
0156 fitparam.rawdata = data;
0157 fitparam.area = Estimates(1);
0158 fitparam.mu = Estimates(2);
0159 fitparam.sigma = Estimates(3);
0160 i = 1;
0161 if autofitDC
0162     fitparam.DC = Estimates(3+i);
0163     i = i + 1;
0164 else
0165     fitparam.DC = fixedvals(1);
0166 end
0167 if autofitGrad
0168     fitparam.bg_gradient = Estimates(3+i);
0169     i = i + 1;
0170 else
0171     fitparam.bg_gradient = fixedvals(2);
0172 end
0173 if autofitAssym
0174     fitparam.Assym_factor = Estimates(3+i);
0175 else
0176     fitparam.Assym_factor = fixedvals(3);
0177 end
0178 fitparam.final_fit_val = fval;
0179 
0180 
0181 
0182 varargout{1} = fitparam;
0183 if nargout > 1
0184     varargout{2} = fval;
0185 end
0186 if nargout > 2
0187     gaussianfit = ones(size(data));
0188     for i = 1:datasize
0189         c = x(i);
0190         gaussianfit(i) = fitparam.area * exp(-0.5*((c-fitparam.mu)./((1+sign(c-fitparam.mu)*fitparam.Assym_factor)*fitparam.sigma)).^2) / sqrt(2*pi*fitparam.sigma^2) + fitparam.bg_gradient*c + fitparam.DC;
0191     end
0192     varargout{3} = gaussianfit;
0193 end
0194 if nargout > 3
0195     % Calculate error in sigma
0196     er=[];
0197     errscale = ones(size(Estimates));
0198     for perc=0.90:0.001:1.1;
0199         errscale(3) = perc; % change the sigma value and see the fit.
0200         er(end+1) = myfit(Estimates.*errscale,x,data,fixedvals);
0201     end
0202     % Normalise
0203     er = er./min(er);
0204     % threshold a 5% change in the error function;
0205     ind = find(er<1.05);
0206     perc = 0.90:0.001:1.1;
0207     % percerror
0208     sigmaerror = (perc(ind(end))-1);
0209     
0210     varargout{4} = sigmaerror;
0211 end
0212 
0213 if DEBUG || plotfig
0214     gaussianfit = ones(size(data));
0215     for i = 1:datasize
0216         c = x(i);
0217         gaussianfit(i) = fitparam.area * exp(-0.5*((c-fitparam.mu)./((1+sign(c-fitparam.mu)*fitparam.Assym_factor)*fitparam.sigma)).^2) / sqrt(2*pi*fitparam.sigma^2) + fitparam.bg_gradient*c + fitparam.DC;
0218     end
0219     figure(233);
0220     plot(x,data, '.-r');
0221     hold on;
0222     plot(x,gaussianfit, '-b');
0223     hold off;
0224 %     title(sprintf('Fitting STD error %g (\\sigma = %f)',...
0225 %         fittingerror(end),fitsigmas(end)));
0226 end
0227 
0228 
0229 
0230 
0231 function sse=myfit(params, x, Dist, fixedvals)
0232 
0233 % if length(params) > 0
0234     Afit = params(1);
0235 % else
0236 %     Afit = 1;
0237 % end
0238 % if length(params) > 1
0239     mufit = params(2);
0240 % else
0241 %     mufit = length(x)/2;
0242 % end
0243 % if length(params) > 2
0244     sigmafit = params(3);
0245 % else
0246 %     sigmafit = length(x)/10;
0247 % end
0248 
0249 i = 1;
0250 if isnan(fixedvals(1))
0251     DCfit = params(3+i);
0252     i = i + 1;
0253 else
0254     DCfit = fixedvals(1);
0255 end
0256 
0257 if isnan(fixedvals(2))
0258     linefit = params(3+i);
0259     i = i + 1;
0260 else
0261     linefit = fixedvals(2);
0262 end
0263 
0264 if isnan(fixedvals(3))
0265     Asym = params(3+i);
0266 else
0267     Asym = fixedvals(3);
0268 end
0269 
0270 fittedcurve = Afit * exp(-0.5*((x-mufit)./((1+sign(x-mufit)*Asym)*sigmafit)).^2) / sqrt(2*pi*sigmafit^2) + (linefit*x) + DCfit;
0271 
0272 sse = sum((fittedcurve - Dist).^2)/length(Dist);