diff --git a/@MyFit/MyFit.m b/@MyFit/MyFit.m index 0f27ed8..1738bcc 100644 --- a/@MyFit/MyFit.m +++ b/@MyFit/MyFit.m @@ -1,262 +1,274 @@ classdef MyFit < handle properties Gui Data; Fit; Parser; fit_name='linear' init_params=[]; scale_init=[]; lim_lower; lim_upper; FitStruct; Fitdata; coeffs; enable_gui=1; enable_plot; plot_handle; hline_init; end properties (Dependent=true) fit_function; fit_tex; fit_params; fit_param_names; valid_fit_names; n_params; scaled_params; init_param_fun; end methods function this=MyFit(varargin) createFitStruct(this); createParser(this); parse(this.Parser,varargin{:}); parseInputs(this); if ismember('Data',this.Parser.UsingDefaults) &&... ~ismember('x',this.Parser.UsingDefaults) &&... ~ismember('y',this.Parser.UsingDefaults) this.Data.x=this.Parser.Results.x; this.Data.y=this.Parser.Results.y; end - genInitParams(this); + if ~isempty(this.Data.x) || ~isempty(this.Data.y) + genInitParams(this); + else + this.init_params=ones(1,this.n_params); + end + this.scale_init=ones(1,this.n_params); if this.enable_gui createGui(this); end end %Creates the GUI of MyFit createGui(this); function delete(this) if this.enable_gui set(this.Gui.Window,'CloseRequestFcn',''); %Deletes the figure delete(this.Gui.Window); %Removes the figure handle to prevent memory leaks this.Gui=[]; end end %Close figure callback simply calls delete function for class function closeFigure(this,~,~) delete(this); end function createParser(this) p=inputParser; addParameter(p,'fit_name','linear',@ischar) addParameter(p,'Data',MyTrace()); addParameter(p,'Fit',MyTrace()); addParameter(p,'x',[]); addParameter(p,'y',[]); addParameter(p,'enable_gui',1); addParameter(p,'enable_plot',0); addParameter(p,'plot_handle',[]); this.Parser=p; end %Sets the class variables to the inputs from the inputParser. function parseInputs(this) for i=1:length(this.Parser.Parameters) %Takes the value from the inputParser to the appropriate %property. if isprop(this,this.Parser.Parameters{i}) this.(this.Parser.Parameters{i})=... this.Parser.Results.(this.Parser.Parameters{i}); end end end function fitTrace(this) this.Fit.x=linspace(min(this.Data.x),max(this.Data.x),1e3); switch this.fit_name case 'linear' this.coeffs=polyfit(this.Data.x,this.Data.y,1); + this.Fit.y=polyval(this.coeffs,this.Fit.x); case 'quadratic' this.coeffs=polyfit(this.Data.x,this.Data.y,2); + this.Fit.y=polyval(this.coeffs,this.Fit.x); case {'exponential','gaussian'} this.doFit this.coeffs=coeffvalues(this.Fitdata); this.Fit.y=this.Fitdata(this.Fit.x); otherwise this.doFit; this.Fit.y=this.Fitdata(this.Fit.x); this.coeffs=coeffvalues(this.Fitdata); end this.init_params=this.coeffs; this.scale_init=ones(1,this.n_params); - updateGui(this); + if this.enable_gui; updateGui(this); end + if this.enable_plot; plotFit(this); end end function doFit(this) this.Fitdata=fit(this.Data.x,this.Data.y,this.fit_function,... 'Lower',this.lim_lower,'Upper',this.lim_upper,... 'StartPoint',this.init_params); end + function plotFit(this) + plot(this.plot_handle,this.Fit.x,this.Fit.y); + end + function createFitStruct(this) %Adds fits - addFit(this,'linear','a*x_b','$$ax+b$$',{'a','b'},... + addFit(this,'linear','a*x+b','$$ax+b$$',{'a','b'},... {'Gradient','Offset'}) addFit(this,'quadratic','a*x^2+b*x+c','$$ax^2+bx+c$$',... {'a','b','c'},{'Quadratic coeff.','Linear coeff.','Offset'}); addFit(this,'gaussian','a*exp(-((x-c)/b)^2/2)+d',... '$$ae^{-\frac{(x-c)^2}{2b^2}}+d$$',{'a','b','c','d'},... {'Amplitude','Width','Center','Offset'}); addFit(this,'lorentzian','a/(pi)*(b/((x-c)^2+b^2))',... '$$\frac{a}{1+\frac{(x-c)^2}{b^2}}+d$$',{'a','b','c','d'},... {'Amplitude','Width','Center','Offset'}); addFit(this,'exponential','a*exp(b*x)+c',... '$$ae^{bx}+c$$',{'a','b','c'},... {'Amplitude','Rate','Offset'}); end function updateGui(this) %Converts the scale variable to the value between 0 and 100 %necessary for the slider slider_vals=25*log10(this.scale_init)+50; for i=1:this.n_params set(this.Gui.(sprintf('edit_%s',this.fit_params{i})),... 'String',sprintf('%3.3e',this.scaled_params(i))); set(this.Gui.(sprintf('slider_%s',this.fit_params{i})),... 'Value',slider_vals(i)); end end %Adds a fit to the list of fits function addFit(this,fit_name,fit_function,fit_tex,fit_params,... fit_param_names) this.FitStruct.(fit_name).fit_function=fit_function; this.FitStruct.(fit_name).fit_tex=fit_tex; this.FitStruct.(fit_name).fit_params=fit_params; this.FitStruct.(fit_name).fit_param_names=fit_param_names; %Generates the anonymous fit function from the above args=['@(x,', strjoin(fit_params,','),')']; anon_fit_fun=str2func(vectorize([args,fit_function])); this.FitStruct.(fit_name).anon_fit_fun=anon_fit_fun; end function genInitParams(this) switch this.fit_name case 'exponential' [this.init_params,this.lim_lower,this.lim_upper]=... initParamExponential(this.Data.x,this.Data.y); case 'gaussian' [this.init_params,this.lim_lower,this.lim_upper]=... initParamGaussian(this.Data.x,this.Data.y); otherwise this.init_params=ones(1,this.n_params); end end function slider_Callback(this, param_ind, hObject, ~) %Gets the value from the slider scale=get(hObject,'Value'); %Updates the scale with a new value this.scale_init(param_ind)=10^((scale-50)/25); %Updates the edit box with the new value from the slider set(this.Gui.(sprintf('edit_%s',this.fit_params{param_ind})),... 'String',sprintf('%3.3e',this.scaled_params(param_ind))); if this.enable_plot; plotInitFun(this); end end function edit_Callback(this, hObject, ~) init_param=str2double(get(hObject,'String')); tag=get(hObject,'Tag'); %Finds the index where the fit_param name begins (convention is %after the underscore) fit_param=tag((strfind(tag,'_')+1):end); param_ind=strcmp(fit_param,this.fit_params); %Updates the slider to be such that the scaling is 1 set(this.Gui.(sprintf('slider_%s',fit_param)),... 'Value',50); %Updates the correct initial parameter this.init_params(param_ind)=init_param; if this.enable_plot; plotInitFun(this); end end function plotInitFun(this) %Substantially faster than any alternative - generating %anonymous functions is very cpu intensive. x_vec=linspace(min(this.Data.x),max(this.Data.x),1000); input_cell=num2cell(this.scaled_params); y_vec=feval(this.FitStruct.(this.fit_name).anon_fit_fun,x_vec,... input_cell{:}); if isempty(this.hline_init) this.hline_init=plot(this.plot_handle,x_vec,y_vec); else set(this.hline_init,'XData',x_vec,'YData',y_vec); end end function set.fit_name(this,fit_name) assert(ischar(fit_name),'The fit name must be a string'); this.fit_name=lower(fit_name); end function valid_fit_names=get.valid_fit_names(this) valid_fit_names=fieldnames(this.FitStruct); end function fit_function=get.fit_function(this) assert(ismember(this.fit_name,this.valid_fit_names),... '%s is not a supported fit name',this.fit_name); fit_function=this.FitStruct.(this.fit_name).fit_function; end function fit_tex=get.fit_tex(this) assert(ismember(this.fit_name,this.valid_fit_names),... '%s is not a supported fit name',this.fit_name); fit_tex=this.FitStruct.(this.fit_name).fit_tex; end function fit_params=get.fit_params(this) assert(ismember(this.fit_name,this.valid_fit_names),... '%s is not a supported fit name',this.fit_name); fit_params=this.FitStruct.(this.fit_name).fit_params; end function fit_param_names=get.fit_param_names(this) assert(ismember(this.fit_name,this.valid_fit_names),... '%s is not a supported fit name',this.fit_name); fit_param_names=this.FitStruct.(this.fit_name).fit_param_names; end function scaled_params=get.scaled_params(this) scaled_params=this.scale_init.*this.init_params; end function n_params=get.n_params(this) n_params=length(this.fit_params); end end end \ No newline at end of file diff --git a/Utility functions/fitExponential.m b/Utility functions/fitExponential.m index 09c71f2..c7e5e67 100644 --- a/Utility functions/fitExponential.m +++ b/Utility functions/fitExponential.m @@ -1,30 +1,17 @@ -function fitdata=fitExponential(x,y,p_in) -x=x(:); -y=y(:); -ffun='a*exp(b*x)+c'; +function fitdata=fitExponential(x,y,varargin) + +p=createFitParser(3); -%Setting upper and lower limits -[amp_max,ind_max]=max(y); -[amp_min,ind_min]=min(y); +parse(p,x,y,varargin{:}); -lim_upper=[Inf,Inf,Inf]; -lim_lower=-lim_upper; +%Converts to column vectors +x=p.Results.x(:); +y=p.Results.y(:); -if abs(amp_max)>abs(amp_min) - lim_upper(1)=Inf; - lim_lower(1)=0; -else - lim_upper(1)=0; - lim_lower(1)=-Inf; -end +assert(length(x)==length(y),'The length of x and y must be equal'); -if ind_max>ind_min - lim_upper(2)=Inf; - lim_lower(2)=0; -else - lim_upper(2)=0; - lim_lower(2)=-Inf; -end +ffun='a*exp(b*x)+c'; -fitdata=fit(x,y,ffun,'Lower',lim_lower,'Upper',lim_upper,'StartPoint',p_in); +fitdata=fit(x,y,ffun,'Lower',p.Results.Lower,'Upper',p.Results.Upper,... + 'StartPoint',p.Results.StartPoint); end \ No newline at end of file diff --git a/Utility functions/initParamExponential.m b/Utility functions/initParamExponential.m index 5cca33a..d6fa3f0 100644 --- a/Utility functions/initParamExponential.m +++ b/Utility functions/initParamExponential.m @@ -1,44 +1,46 @@ function [p_in,lim_lower,lim_upper]=initParamExponential(x,y) %Assumes form a*exp(-bx)+c %Setting upper and lower limits [amp_max,ind_max]=max(y); [amp_min,ind_min]=min(y); lim_upper=[Inf,Inf,Inf]; lim_lower=-lim_upper; if abs(amp_max)>abs(amp_min) lim_upper(1)=Inf; lim_lower(1)=0; else lim_upper(1)=0; lim_lower(1)=-Inf; end if ind_max>ind_min lim_upper(2)=Inf; lim_lower(2)=0; else lim_upper(2)=0; lim_lower(2)=-Inf; end %Method for estimating initial parameters taken from -%http://www.matrixlab-examples.com/exponential-regression.html +%http://www.matrixlab-examples.com/exponential-regression.html. Some +%modifications to account for negative y values + n=length(x); y2=log(y); j=sum(x); k=sum(y2); l=sum(x.^2); r2=sum(x .* y2); p_in(2)=(n * r2 - k * j)/(n * l - j^2); p_in(1)=exp((k-p_in(2)*j)/n); if abs(amp_max)>abs(amp_min) p_in(3)=min(y); else p_in(3)=max(y); end end \ No newline at end of file diff --git a/Utility functions/initParamGaussian.m b/Utility functions/initParamGaussian.m index 0ea4532..379bcba 100644 --- a/Utility functions/initParamGaussian.m +++ b/Utility functions/initParamGaussian.m @@ -1,32 +1,38 @@ function [p_in,lim_lower,lim_upper]=initParamGaussian(x,y) %Assumes a*exp(-((x-c)/b)^2/2)+d - remember matlab orders the fit %parameters alphabetically -[amp,ind_max]=max(y); -center=x(ind_max); -width=sqrt(sum(y.*(x-center).^2)/sum(y)); -bg=median(y); -p_in=[amp,width,center, bg]; +[amp_max,ind_max]=max(y); +[amp_min,ind_min]=min(y); lim_upper=[Inf,Inf,Inf,Inf]; lim_lower=-lim_upper; if abs(amp_max)>abs(amp_min) + amp=amp_max; + center=x(ind_max); lim_upper(1)=Inf; lim_lower(1)=0; else + amp=amp_min; + center=x(ind_min); lim_upper(1)=0; lim_lower(1)=-Inf; end +width=sqrt(sum(y.*(x-center).^2)/sum(y)); +bg=median(y); + %Sets the lower limit on width to zero lim_lower(2)=0; %Sets the upper limit on width to 100 times the range of the data lim_upper(2)=100*range(x); %Sets upper and lower limit on the center lim_lower(3)=min(x)/2; lim_upper(3)=max(x)*2; +p_in=[amp_max,width,center, bg]; + end \ No newline at end of file diff --git a/Utility functions/testfit.m b/Utility functions/testfit.m new file mode 100644 index 0000000..45dc8cd --- /dev/null +++ b/Utility functions/testfit.m @@ -0,0 +1,23 @@ +%Testing tool for MyFit +clear +x_vec=linspace(0,10,100); + +testFit=MyFit('fit_name','exponential','enable_gui',0); +params=cell(1,testFit.n_params); +for i=1:testFit.n_params + params{i}=rand; +end +params{2}=-rand; + +y_vec=testFit.FitStruct.(testFit.fit_name).anon_fit_fun(x_vec+rand(size(x_vec)),params{:}); +figure(1) +ax=gca; +plot(x_vec,y_vec,'x'); +hold on +testFit.plot_handle=ax; +testFit.enable_plot=1; +testFit.Data.x=x_vec; +testFit.Data.y=y_vec; +testFit.genInitParams; +testFit.init_params +testFit.fitTrace;