how to speed up gaussian fitting
19 次查看(过去 30 天)
显示 更早的评论
Greetings,
I am working on a piece of code that takes in some data and fits a gaussian with a constant offset to it. I first wrote a quick way to fit the gaussian that wasn't a true least squares method but it ran very fast ( could fit all my data in a matter of 10s of seconds). I now would like to go back and take the parameters that the quick fit gavem and use them to initialize a true least squares fit. What I've got it as follows :
fo=fitoptions('Method','NonlinearLeastSquares',...
'Lower',[0 0 0 0],...
'Upper',[Inf Inf Inf Inf],...
'StartPoint', [ Ai Bi Ci Di]); % - setup the fitting options
ft=fittype( 'a+b*exp(-(x-c)*(x-c)/(2*d*d))','options',fo); % make the fit type
[CURVE,GOF] = fit(x,y,ft); % - do the fit!
and I do this for each set of data points I want to fit. I have roughly 10000-20000 fits that I would like to do and this code is taking forever ( as in hours!) am I missing a way to speed this up? I could constrainthe upper and lower bounds slightly more, would that make a big difference?
For completeness I have the entire function below :
function [xc, resid, yfit, A,diff] = recenterg(x, y, xc, width)
% Fit y = c + m* x + A * exp(- (x-xc)^2 / (2 * width^2)) )
% c -> average background intensity
% m -> slope of background intensity
% A -> Amplitude of peak
% xc -> peak center
% width -> width of peak (essential psf)
%
% Input -> x,y -> vector of positions
% xc -> initial guess for peak position
% width -> width of psf
%
%
% Outputs ->
% xc -> refined peak center position
% resid -> rms of fitted intensities
% yfit -> fitted values of y
% A -> fitted peak amplitude
%
%
% Pick initial fit values so loop excutes at least once
Niter = 0;
Niter_max = 5; % Limit to 5 iterations
xc_cut = 0.01; % Limit to fit precision of 0.01 pixels
xc_min = min(x) + width; % Limits on center position
xc_max = max(x) - width;
while (Niter<=Niter_max)
Niter = Niter + 1;
% Perform linear regression which should be fast
% of y = c + m* x + A * f + (A * xc_err) * df/dxc
% c = cbest(1); m = cbest(2) ; A = cbest(3); A * xc_err = cbest(4)
f = exp( - (x - xc).^2 / (2*width^2))/ width; % Peak function
M = [ones(size(x)), x, f, f.*(x-xc)/width^2];
cbest = (M'*M) \ M' * y;
rmsd = std(y - M*cbest); % RMS error
xc_err = cbest(4)/cbest(3);
if (cbest(3)> 2*rmsd)&(abs(xc_err)<width)
% Peak signal to noise > 2 sigma
% Correction in peak position < width of peak
xc = xc + xc_err; % Update center position
yfit = M*cbest;
resid = sqrt(mean((y - yfit).^2));
A = cbest(3);
% End fit if convert to within precision of xc_cut
if (abs(xc_err)<xc_cut)
Niter = Niter_max + 1;
end
else
%xc = xc + 0.3 * xc_err;
Niter = Niter_max +1 ; % Abort fit
xc = NaN;
resid = NaN;
yfit = NaN;
A = cbest(3);
end
% ok ,ran with the quick and dirty f it, now take those quick
%values and use them to initialize a true Gaussian fit
%****** START HERE GAUSSIAN FITTING PROCEDURE
% first check if the value was NaN or 0 if it was then leave it as
% it
if (isnan(xc)) % when the quick fit doesnt work the gaussin fit won't either so give them NaN
xc = NaN;
resid = NaN;
yfit = NaN;
A = cbest(3);
diff=NaN;
else
Ai= (sum(y(1:4))+sum(y(end-3:end)))./8; % a rough guess of the background
Bi= A; % - guess for the amplitude
Ci= xc; % - from the quick fit
Di= 1.9 ; % the width as a first guess
fo=fitoptions('Method','NonlinearLeastSquares',...
'Lower',[0 0 0 0],...
'Upper',[Inf Inf Inf Inf],...
'StartPoint', [ Ai Bi Ci Di]); % - setup the fitting options
ft=fittype( 'a+b*exp(-(x-c)*(x-c)/(2*d*d))','options',fo); % make the fit type
[CURVE,GOF] = fit(x,y,ft); % - do the fit!
diff=abs( xc - CURVE.c); % - to keep track of the difference between quick method and real fit
xc = CURVE.c; % - the center position of the fit
yfit = CURVE.b; %- don't really need this but might be useful later
resid = CURVE.d; % - the std of fit
A = CURVE.b ; %- the fit amplitude
end
0 个评论
回答(1 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Problem-Based Optimization Setup 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!