diff --git a/@MyFit/MyFit.m b/@MyFit/MyFit.m index 8fb4127..f39d8ab 100644 --- a/@MyFit/MyFit.m +++ b/@MyFit/MyFit.m @@ -1,1023 +1,1032 @@ classdef MyFit < dynamicprops %Note that dynamicprops classes are handle classes. properties (Access=public) Data; %MyTrace object contains the data to be fitted to init_params=[]; %Contains the initial parameters lim_lower; %Lower limits for fit parameters lim_upper; %Upper limits for fit parameters enable_plot; %If enabled, plots initial parameters in the plot_handle plot_handle; %The handle which fits and init params are plotted in %Calibration values supplied externally CalVals=struct(); init_color='c'; end properties (GetAccess=public, SetAccess=private) Fit; %MyTrace object containing the fit Gui; %Gui handles %Output structures from fit: Fitdata; Gof; FitInfo; %Contains information about all the fits and their parameters FitStruct; coeffs; %The name of the fit used. fit_name='Linear' end properties (Access=private) %Structure used for initializing GUI of userpanel UserGui; Parser; %Input parser for constructor enable_gui=1; hline_init; %Handle for the plotted init values %Private struct used for saving file information when there is no %gui SaveInfo slider_vecs; %Vectors for varying the range of the sliders for different fits end %Dependent variables with no set methods properties (Dependent=true, SetAccess=private) %These grab and set the appropriate field of the FitStruct fit_function; anon_fit_fun; fit_tex; fit_params; fit_param_names; valid_fit_names; n_params; %These are used to create the usergui n_user_fields; user_field_tags; user_field_names; user_field_vals; %Vector used for plotting, depends on the data trace x_vec; %Variables used for saving fullpath; save_path; end %Dependent variables with associated set methods properties (Dependent=true) filename; base_dir; session_name; end %Events for communicating with outside entities events NewFit; NewInitVal; end %Parser function methods (Access=private) %Creates parser for constructor 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',[]); addParameter(p,'base_dir',this.SaveInfo.filename); addParameter(p,'session_name',this.SaveInfo.session_name); addParameter(p,'filename',this.SaveInfo.base_dir); this.Parser=p; end end %Public methods methods (Access=public) %Constructor function function this=MyFit(varargin) %Sets the default parameters for the save directory and %filename. this.SaveInfo.filename='placeholder'; this.SaveInfo.session_name='placeholder'; this.SaveInfo.base_dir=getLocalSettings('measurement_base_dir'); %We first create the FitStruct, which contains all the %information about the available fits. createFitStruct(this); %We now create the parser for parsing the arguments to the %constructor, and parse the variables. createParser(this); parse(this.Parser,varargin{:}); parseInputs(this); %Loads values into the CalVals struct depending on the type of %fit initCalVals(this); %Allows us to load either x/y data or a MyTrace object directly 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 %Sets dummy values for the GUI this.init_params=ones(1,this.n_params); this.lim_lower=-Inf(1,this.n_params); this.lim_upper=Inf(1,this.n_params); %Creates the structure that contains variables for calibration %of fit results createUserGuiStruct(this); %Creates the gui if the flag is enabled. This function is in a %separate file. if this.enable_gui createGui(this) %Generates the slider lookup table genSliderVecs(this); end %If the data is appropriate, generates initial %parameters if validateData(this); genInitParams(this); end end %Deletion function of object function delete(this) if this.enable_gui %Avoids loops set(this.Gui.Window,'CloseRequestFcn',''); %Deletes the figure delete(this.Gui.Window); %Removes the figure handle to prevent memory leaks this.Gui=[]; end if ~isempty(this.hline_init); delete(this.hline_init); end if ~isempty(this.Fit.hlines); delete(this.Fit.hlines{:}); end end %Close figure callback simply calls delete function for class function closeFigure(this,~,~) delete(this); end %Saves the metadata function saveParams(this,varargin) p=inputParser; addParameter(p,'save_user_params',true); addParameter(p,'save_gof',true); parse(p,varargin{:}); %Flags for saving the user parameters or goodness of fit save_user_params=p.Results.save_user_params; save_gof=p.Results.save_gof; assert(~isempty(this.coeffs) && ... length(this.coeffs)==this.n_params,... ['The number of calculated coefficients (%i) is not',... ' equal to the number of parameters (%i).', ... ' Perform a fit before trying to save parameters.'],... length(this.coeffs),this.n_params); %Creates combined strings of form: Linewidth (b), where %Linewidth is the parameter name and b is the parameter tag headers=cellfun(@(x,y) sprintf('%s (%s)',x,y),... this.fit_param_names, this.fit_params,'UniformOutput',0); save_data=this.coeffs; if save_user_params %Creates headers for the user fields user_field_headers=cellfun(@(x,y) ... sprintf('%s. %s',this.UserGui.Fields.(x).parent,y),... this.user_field_tags,this.user_field_names,... 'UniformOutput',0)'; %Appends the user headers and data to the save data headers=[headers, user_field_headers]; save_data=[save_data,this.user_field_vals']; end if save_gof %Appends GOF headers and data to the save data headers=[headers,fieldnames(this.Gof)']; save_data=[save_data,struct2array(this.Gof)]; end %Find out at the end how many columns we have n_columns=length(headers); %Sets the column width. Pads 2 for legibility. col_width=cellfun(@(x) length(x), headers)+2; %Min column width of 24 col_width(col_width<24)=24; %Create the right directories if ~exist(this.base_dir,'dir') mkdir(this.base_dir) end if ~exist(this.save_path,'dir') mkdir(this.save_path) end %We automatically append to the file if it already exists, %otherwise create a new file if exist(this.fullpath,'file') fileID=fopen(this.fullpath,'a'); fprintf('Appending data to %s \n',this.fullpath); else fileID=fopen(this.fullpath,'w'); pre_fmt_str=repmat('%%%is\\t',1,n_columns); fmt_str=sprintf([pre_fmt_str,'\r\n'],col_width); fprintf(fileID,fmt_str,headers{:}); end pre_fmt_str_nmb=repmat('%%%i.15e\\t',1,n_columns); nmb_fmt_str=sprintf([pre_fmt_str_nmb,'\r\n'],col_width); fprintf(fileID,nmb_fmt_str,save_data); fclose(fileID); end %We can load a fit from a file with appropriately formatted columns %We simply load the coefficients from the file into the fit. function loadFit(this,fullfilename,varargin) p=inputParser; addParameter(p,'line_no',1); parse(p,varargin{:}) n=p.Results.line_no; load_table=readtable(fullfilename); load_names=fieldnames(load_table); for i=1:this.n_params this.coeffs(i)=load_table.(load_names{i})(n); end end %This function is used to set the coefficients, to avoid setting it %to a number not equal to the number of parameters function setFitParams(this,coeffs) assert(length(coeffs)==this.n_params,... ['The length of the coefficient vector (currently %i) ',... 'must be equal to the number of parameters (%i)'],... length(this.coeffs),this.n_params) this.coeffs=coeffs; end %Initializes the CalVals structure. function initCalVals(this) switch this.fit_name case 'Lorentzian' %Line spacing is the spacing between all the lines, %i.e. number of lines times the spacing between each %one this.CalVals.line_spacing=1; case 'DoubleLorentzian' this.CalVals.line_spacing=1; end end %Fits the trace using currently set parameters, depending on the %model. function fitTrace(this) %Checks for valid data. assert(validateData(this),... ['The length of x is %d and the length of y is',... ' %d. The lengths must be equal and greater than ',... 'the number of fit parameters to perform a fit'],... length(this.Data.x),length(this.Data.y)) %Checks for valid limits. lim_check=this.lim_upper>this.lim_lower; assert(all(lim_check),... sprintf(['All upper limits must exceed lower limits. ',... 'Check limit %i, fit parameter %s'],find(~lim_check,1),... this.fit_params{find(~lim_check,1)})); switch this.fit_name case 'Linear' %Fits polynomial of order 1 this.coeffs=polyfit(this.Data.x,this.Data.y,1); case 'Quadratic' %Fits polynomial of order 2 this.coeffs=polyfit(this.Data.x,this.Data.y,2); this.Fit.y=polyval(this.coeffs,this.Fit.x); - case {'Lorentzian','LorentzianGrad','Gaussian',... + case 'Lorentzian' + scale_factor=max(this.Data.y); + this.Data.y=this.Data.y/scale_factor; + this.init_params(1)=this.init_params(1)/scale_factor; + this.init_params(4)=this.init_params(4)/scale_factor; + doFit(this); + this.coeffs(1)=this.coeffs(1)*scale_factor; + this.coeffs(4)=this.coeffs(4)*scale_factor; + this.Data.y=this.Data.y*scale_factor; + case {'LorentzianGrad','Gaussian',... 'DoubleLorentzian','DoubleLorentzianGrad',... 'Exponential','Gorodetsky2000',... 'Gorodetsky2000plus'} doFit(this) otherwise error('Selected fit is invalid'); end %This function calculates the fit trace, using this.x_vec as %the x axis calcFit(this); %This function calculates user-defined parameters calcUserParams(this); %Sets the new initial parameters to be the fitted parameters this.init_params=this.coeffs; %Updates the gui if it is enabled if this.enable_gui genSliderVecs(this); updateGui(this); end %Plots the fit if the flag is on if this.enable_plot; plotFit(this); end %Triggers new fit event triggerNewFit(this); end %This function calculates all the user-defined parameters shown in %the GUI. To add more parameters, add them in createUserGuiStruct, %then add them here when you wish to calculate them. function calcUserParams(this) switch this.fit_name case 'Lorentzian' this.mech_lw=this.coeffs(2); %#ok this.mech_freq=this.coeffs(3); %#ok this.Q=this.mech_freq/this.mech_lw; %#ok this.opt_lw=convOptFreq(this,this.coeffs(2)); %#ok this.Qf=this.mech_freq*this.Q; %#ok case 'DoubleLorentzian' this.opt_lw1=convOptFreq(this,this.coeffs(2)); %#ok this.opt_lw2=convOptFreq(this,this.coeffs(5)); %#ok splitting=abs(this.coeffs(6)-this.coeffs(3)); this.mode_split=convOptFreq(this,splitting); %#ok case 'Exponential' if ~isempty(this.coeffs) this.tau=abs(1/this.coeffs(2)); %#ok this.lw=abs(this.coeffs(2)/pi); %#ok end this.Q=pi*this.freq*this.tau; %#ok this.Qf=this.Q*this.freq; %#ok otherwise %If fit_name is not listed, do nothing end end %This function is used to convert the x-axis to frequency. function real_freq=convOptFreq(this,freq) real_freq=freq*this.spacing*this.line_no/this.CalVals.line_spacing; end %This struct is used to generate the UserGUI. Fields are seen under %tabs in the GUI. To create a new tab, you have to enter it under %this.UserGui.Tabs. A tab must have a tab_title and a field to add %Children. To add a field, use the addUserField function. function createUserGuiStruct(this) this.UserGui=struct('Fields',struct(),'Tabs',struct()); switch this.fit_name case 'Lorentzian' %Parameters for the tab relating to mechanics this.UserGui.Tabs.Mech.tab_title='Mech.'; this.UserGui.Tabs.Mech.Children={}; addUserField(this,'Mech','mech_lw','Linewidth (Hz)',1,... 'enable_flag','off') addUserField(this,'Mech','Q',... 'Qualify Factor (x10^6)',1e6,... 'enable_flag','off','conv_factor',1e6) addUserField(this,'Mech','mech_freq','Frequency (MHz)',1e6,... 'conv_factor',1e6, 'enable_flag','off') addUserField(this,'Mech','Qf','Q\times f (10^{14} Hz)',1e14,... 'conv_factor',1e14,'enable_flag','off'); %Parameters for the tab relating to optics this.UserGui.Tabs.Opt.tab_title='Optical'; this.UserGui.Tabs.Opt.Children={}; addUserField(this,'Opt','spacing',... 'Line Spacing (MHz)',1e6,'conv_factor',1e6,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','line_no','Number of lines',10,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','opt_lw','Linewidth (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); case 'DoubleLorentzian' this.UserGui.Tabs.Opt.tab_title='Optical'; this.UserGui.Tabs.Opt.Children={}; addUserField(this,'Opt','spacing',... 'Line Spacing (MHz)',1e6,'conv_factor',1e6,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','line_no','Number of lines',10,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','opt_lw1','Linewidth 1 (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); addUserField(this,'Opt','opt_lw2','Linewidth 2 (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); addUserField(this,'Opt','mode_split',... 'Modal splitting (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); case 'Exponential' this.UserGui.Tabs.Q.tab_title='Q'; this.UserGui.Tabs.Q.Children={}; addUserField(this,'Q','tau','\tau (s)',1,... 'enable_flag','off') addUserField(this,'Q','lw','Linewidth (Hz)',1,... 'enable_flag','off') addUserField(this,'Q','Q',... 'Qualify Factor (x10^6)',1e6,... 'enable_flag','off','conv_factor',1e6) addUserField(this,'Q','freq','Frequency (MHz)',1e6,... 'conv_factor',1e6, 'enable_flag','on',... 'Callback',@(~,~) calcUserParams(this)); addUserField(this,'Q','Qf','Q\times f (10^{14} Hz)',1e14,... 'conv_factor',1e14,'enable_flag','off'); otherwise %Do nothing if there is no defined user parameters end end %Parent is the parent tab for the userfield, tag is the tag given %to the GUI element, title is the text written next to the field, %initial value is the initial value of the property and change_flag %determines whether the gui element is enabled for writing or not. %conv_factor is used to have different units in the field. In the %program, the value is always saved as the bare value. function addUserField(this, parent, tag, title, ... init_val,varargin) %Parsing inputs p=inputParser(); addRequired(p,'Parent'); addRequired(p,'Tag'); addRequired(p,'Title'); addRequired(p,'init_val'); addParameter(p,'enable_flag','on'); addParameter(p,'Callback',''); addParameter(p,'conv_factor',1); parse(p,parent,tag,title,init_val,varargin{:}); tag=p.Results.Tag; %Populates the UserGui struct this.UserGui.Fields.(tag).parent=p.Results.Parent; this.UserGui.Fields.(tag).title=p.Results.Title; this.UserGui.Fields.(tag).init_val=p.Results.init_val; this.UserGui.Fields.(tag).enable_flag=... p.Results.enable_flag; this.UserGui.Fields.(tag).conv_factor=p.Results.conv_factor; this.UserGui.Fields.(tag).Callback=... p.Results.Callback; this.UserGui.Tabs.(p.Results.Parent).Children{end+1}=tag; %Adds the new property to the class addUserProp(this, tag); end %Every user field has an associated property, which is added by %this function. The get and set functions are set to use the GUI %through the getUserVal and setUserVal functions if the GUI is %enabled. function addUserProp(this,tag) prop=addprop(this,tag); if this.enable_gui prop.GetMethod=@(this) getUserVal(this,tag); prop.SetMethod=@(this, val) setUserVal(this, val, tag); prop.Dependent=true; end end %This function gets the value of the userprop from the GUI. The GUI %is the single point of truth function val=getUserVal(this, tag) conv_factor=this.UserGui.Fields.(tag).conv_factor; val=str2double(this.Gui.([tag,'Edit']).String)*conv_factor; end %As above, but instead we set the GUI through setting the property function setUserVal(this, val, tag) conv_factor=this.UserGui.Fields.(tag).conv_factor; this.Gui.([tag,'Edit']).String=num2str(val/conv_factor); end %Generates model-dependent initial parameters, lower and upper %boundaries. function genInitParams(this) assert(validateData(this), ['The data must be vectors of',... ' equal length greater than the number of fit parameters.',... ' Currently the number of fit parameters is %d, the',... ' length of x is %d and the length of y is %d'],... this.n_params,length(this.Data.x),length(this.Data.y)); %Cell for putting parameters in to be interpreted in the %parser. Element 1 contains the init params, Element 2 contains %the lower limits and Element 3 contains the upper limits. params={}; switch this.fit_name case 'Exponential' [params{1},params{2},params{3}]=... initParamExponential(this.Data.x,this.Data.y); case 'Gaussian' [params{1},params{2},params{3}]=... initParamGaussian(this.Data.x,this.Data.y); case 'Lorentzian' [params{1},params{2},params{3}]=... initParamLorentzian(this.Data.x,this.Data.y); case 'LorentzianGrad' [params{1},params{2},params{3}]=... initParamLorentzianGrad(this.Data.x,this.Data.y); case 'DoubleLorentzian' [params{1},params{2},params{3}]=... initParamDblLorentzian(this.Data.x,this.Data.y); case 'DoubleLorentzianGrad' [params{1},params{2},params{3}]=... initParamDblLorentzianGrad(this.Data.x,this.Data.y); case 'Gorodetsky2000' [params{1},params{2},params{3}]=... initParamGorodetsky2000(this.Data.x,this.Data.y); case 'Gorodetsky2000plus' [params{1},params{2},params{3}]=... initParamGorodetsky2000Plus(this.Data.x,this.Data.y); end %Validates the initial parameters p=createFitParser(this.n_params); parse(p,params{:}); %Loads the parsed results into the class variables this.init_params=p.Results.init_params; this.lim_lower=p.Results.lower; this.lim_upper=p.Results.upper; %Plots the fit function with the new initial parameters if this.enable_plot; plotInitFun(this); end %Updates the GUI and creates new lookup tables for the init %param sliders if this.enable_gui genSliderVecs(this); updateGui(this); end end %Calculates the trace object for the fit function calcFit(this) this.Fit.x=this.x_vec; input_coeffs=num2cell(this.coeffs); this.Fit.y=this.anon_fit_fun(this.Fit.x,input_coeffs{:}); end %Plots the trace contained in the Fit MyTrace object after %calculating the new values function plotFit(this,varargin) calcFit(this); assert((isa(this.plot_handle,'matlab.graphics.axis.Axes')||... isa(this.plot_handle,'matlab.ui.control.UIAxes')),... 'plot_handle property must be defined to valid axis in order to plot') this.Fit.plot(this.plot_handle,varargin{:}); clearInitFun(this); end %Clears the plots function clearFit(this) cellfun(@(x) delete(x), this.Fit.hlines); clearInitFun(this); this.Fit.hlines={}; end function clearInitFun(this) delete(this.hline_init); this.hline_init=[]; end %Function for plotting fit model with current initial parameters. function plotInitFun(this) %Substantially faster than any alternative - generating %anonymous functions is very cpu intensive. input_cell=num2cell(this.init_params); y_vec=feval(this.FitStruct.(this.fit_name).anon_fit_fun,... this.x_vec,input_cell{:}); if isempty(this.hline_init) this.hline_init=plot(this.plot_handle,this.x_vec,y_vec,... 'Color',this.init_color); else set(this.hline_init,'XData',this.x_vec,'YData',y_vec); end end end %Callbacks methods %Save fit function callback function saveFitCallback(this,~,~) assert(~isempty(this.base_dir),'Save directory is not specified'); assert(ischar(this.base_dir),... ['Save directory is not specified.',... ' Should be of type char but is %s.'], ... class(this.base_dir)) this.Fit.save('name',this.filename,... 'base_dir',this.save_path) end %Callback for saving parameters function saveParamCallback(this,~,~) saveParams(this); end function genSliderVecs(this) %Return values of the slider slider_vals=1:101; %Default scaling vector def_vec=10.^((slider_vals-51)/50); %Sets the cell to the default value for i=1:this.n_params this.slider_vecs{i}=def_vec*this.init_params(i); set(this.Gui.(sprintf('Slider_%s',this.fit_params{i})),... 'Value',50); end %Here we can put special cases such that the sliders can behave %differently for different fits. if validateData(this) switch this.fit_name case 'Lorentzian' %We choose to have the slider go over the range of %the x-values of the plot for the center of the %Lorentzian. this.slider_vecs{3}=... linspace(this.x_vec(1),this.x_vec(end),101); %Find the index closest to the init parameter [~,ind]=... min(abs(this.init_params(3)-this.slider_vecs{3})); %Set to ind-1 as the slider goes from 0 to 100 set(this.Gui.(sprintf('Slider_%s',... this.fit_params{3})),'Value',ind-1); end end end %Callback functions for sliders in GUI. Uses param_ind to find out %which slider the call is coming from, this was implemented to %speed up the callback. Note that this gets triggered whenever the %value of the slider is changed. function sliderCallback(this, param_ind, hObject, ~) %Gets the value from the slider val=get(hObject,'Value'); %Find out if the current slider value is correct for the %current init param value. If so, do not change anything. This %is required as the callback also gets called when the slider %values are changed programmatically [~,ind]=... min(abs(this.init_params(param_ind)-this.slider_vecs{param_ind})); if ind~=(val+1) %Updates the scale with a new value from the lookup table this.init_params(param_ind)=... this.slider_vecs{param_ind}(val+1); %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.init_params(param_ind))); if this.enable_plot; plotInitFun(this); end end end %Callback function for edit boxes in GUI function editCallback(this, hObject, ~) init_param=str2double(hObject.String); param_ind=str2double(hObject.Tag); %Centers the slider set(this.Gui.(sprintf('Slider_%s',this.fit_params{param_ind})),... 'Value',50); %Updates the correct initial parameter this.init_params(param_ind)=init_param; if this.enable_plot; plotInitFun(this); end %Triggers event for new init values triggerNewInitVal(this); %Generate the new slider vectors genSliderVecs(this); end %Callback function for editing limits in the GUI function limEditCallback(this, hObject,~) lim = str2double(hObject.String); %Regexp finds type (lower or upper bound) and index expr = '(?Upper|Lower)(?\d+)'; s=regexp(hObject.Tag,expr,'names'); ind=str2double(s.ind); switch s.type case 'Lower' this.lim_lower(ind)=lim; case 'Upper' this.lim_upper(ind)=lim; otherwise error('%s is not properly named for assignment of limits',... hObject.Tag); end end %Callback function for analyze button in GUI. Checks if the data is %ready for fitting. function analyzeCallback(this, ~, ~) fitTrace(this); end %Callback for clearing the fits on the axis. function clearFitCallback(this,~,~) clearFit(this); end %Callback function for generate init parameters button. function initParamCallback(this,~,~) genInitParams(this); end end %Private methods methods(Access=private) %Creates the GUI of MyFit, in separate file. createGui(this); %Creates a panel for the GUI, in separate file createTab(this, tab_tag, bg_color, button_h); %Creats two vboxes (from GUI layouts) to display values of %quantities createUnitBox(this, bg_color, h_parent, name); %Creates an edit box inside a UnitDisp for showing label and value of %a quantity. Used in conjunction with createUnitBox createUnitDisp(this,varargin); %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 %Does the fit with the currently set parameters function doFit(this) ft=fittype(this.fit_function,'coefficients',this.fit_params); opts=fitoptions('Method','NonLinearLeastSquares',... 'Lower',this.lim_lower,... 'Upper',this.lim_upper,... 'StartPoint',this.init_params,... 'MaxFunEvals',2000,... 'MaxIter',2000,... 'TolFun',1e-9,... 'TolX',1e-9); %Fits with the below properties. Chosen for maximum accuracy. [this.Fitdata,this.Gof,this.FitInfo]=... fit(this.Data.x,this.Data.y,ft,opts); %Puts the coeffs into the class variable. this.coeffs=coeffvalues(this.Fitdata); end %Triggers the NewFit event such that other objects can use this to %e.g. plot new fits function triggerNewFit(this) notify(this,'NewFit'); end %Triggers the NewInitVal event function triggerNewInitVal(this) notify(this,'NewInitVal'); end %Creates the struct used to get all things relevant to the fit %model. Ensure that fit parameters are listed alphabetically, as %otherwise the anon_fit_fun will not work properly. function createFitStruct(this) %Adds fits 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','1/pi*a*b/2/((x-c)^2+(b/2)^2)+d',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+d$$',... {'a','b','c','d'},... {'Amplitude','Width','Center','Offset'}); addFit(this,'LorentzianGrad','1/pi*a*b/2/((x-c)^2+(b/2)^2)+d*(x-c)+e',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+d*(x-b)+e$$',... {'a','b','c','d','e'},... {'Amplitude','Width','Center','Gradient','Offset'}); addFit(this,'Exponential','a*exp(b*x)+c',... '$$ae^{bx}+c$$',{'a','b','c'},... {'Amplitude','Rate','Offset'}); addFit(this,'DoubleLorentzian',... '1/pi*b/2*a/((x-c)^2+(b/2)^2)+1/pi*e/2*d/((x-f)^2+(e/2)^2)+g',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+\frac{d}{\pi}\frac{e/2}{(x-f)^2+(e/2)^2}+g$$',... {'a','b','c','d','e','f','g'},... {'Amplitude 1','Width 1','Center 1','Amplitude 2',... 'Width 2','Center 2','Offset'}); addFit(this,'DoubleLorentzianGrad',... '1/pi*b/2*a/((x-c)^2+(b/2)^2)+1/pi*e/2*d/((x-f)^2+(e/2)^2)+g*(x-c)+h',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+\frac{d}{\pi}\frac{e/2}{(x-f)^2+(e/2)^2}+g*x+h$$',... {'a','b','c','d','e','f','g','h'},... {'Amplitude 1','Width 1','Center 1','Amplitude 2',... 'Width 2','Center 2','Gradient','Offset'}); addFit(this,'Gorodetsky2000',... ['a*abs( (k0^2/4 - kex^2/4 + gamma^2/4 - (x-b).^2 + i*k0.*(x-b))',... './( (k0 + kex)^2/4 + gamma^2/4 - (x-b).^2 + i.*(x-b)*(k0 + kex) )).^2+c*(x-b)'],... ['$$a\left|\frac{\kappa_0^2/4-\kappa_{ex}^2/4+\gamma^2/4-(x-b)^2+i\kappa_0(x-b)/2}',... '{(\kappa_0+\kappa_{ex})^2/4+\gamma^2/4-(x-b)^2+i(x-b)(\kappa_0+\kappa_{ex})}\right|^2$$+c(x-b)'],... { 'a','b','c','gamma','k0', 'kex'},... {'Background','Center','BG Slope','Mode splitting',... 'Intrinsic','Extrinsic'}); addFit(this,'Gorodetsky2000plus',... ['(a+c*x+d*x.^2+e*x.^2)*abs( (k0^2/4 - kex^2/4 + gamma^2/4 - (x-b).^2 + i*k0.*(x-b))',... './( (k0 + kex)^2/4 + gamma^2/4 - (x-b).^2 + i.*(x-b)*(k0 + kex) )).^2+f'],... ['$$\left[a+cx+dx^2\right]\left|\frac{\kappa_0^2/4-\kappa_{ex}^2/4+\gamma^2/4-(x-b)^2+i\kappa_0(x-b)/2}',... '{(\kappa_0+\kappa_{ex})^2/4+\gamma^2/4-(x-b)^2+i(x-b)(\kappa_0+\kappa_{ex})}\right|^2$$'],... { 'a','b','c','d','e','f','gamma','k0', 'kex'},... {'Baseline','Center','Lin. Coeff','Quad. Coeff','Cubic Coeff.',... 'Background','Mode splitting','Intrinsic','Extrinsic'}); end %Adds a fit to the list of fits. See above for real examples %fit_function: the function used to fit to in MATLAB form %fit_tex: the fit function written in tex for display in the GUI %fit_params: the fit parameters %fit_param_names: longer names of fit parameters for GUI 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,','),')']; this.FitStruct.(fit_name).anon_fit_fun=... str2func(vectorize([args,fit_function])); end %Updates the GUI if the edit or slider boxes are changed from %elsewhere. function updateGui(this) for i=1:this.n_params str=this.fit_params{i}; set(this.Gui.(sprintf('Edit_%s',str)),... 'String',sprintf('%3.3e',this.init_params(i))); set(this.Gui.(sprintf('Lim_%s_upper',str)),... 'String',sprintf('%3.3e',this.lim_upper(i))) set(this.Gui.(sprintf('Lim_%s_lower',str)),... 'String',sprintf('%3.3e',this.lim_lower(i))) end end %Checks if the class is ready to perform a fit function bool=validateData(this) bool=~isempty(this.Data.x) && ~isempty(this.Data.y) && ... length(this.Data.x)==length(this.Data.y) && ... length(this.Data.x)>=this.n_params; end end % Set function for nondependent variable methods %Set function for fit_name. function set.fit_name(this,fit_name) assert(ischar(fit_name),'The fit name must be a string'); %Capitalizes the first letter fit_name=[upper(fit_name(1)),lower(fit_name(2:end))]; %Checks it is a valid fit name ind=strcmpi(fit_name,this.valid_fit_names);%#ok assert(any(ind),'%s is not a supported fit name',fit_name); this.fit_name=this.valid_fit_names{ind}; %#ok end end % Get functions for dependent variables methods %Generates the valid fit names function valid_fit_names=get.valid_fit_names(this) valid_fit_names=fieldnames(this.FitStruct); end %Grabs the correct fit function from FitStruct function fit_function=get.fit_function(this) fit_function=this.FitStruct.(this.fit_name).fit_function; end %Grabs the correct tex string from FitStruct function fit_tex=get.fit_tex(this) fit_tex=this.FitStruct.(this.fit_name).fit_tex; end %Grabs the correct anon fit function from FitStruct function anon_fit_fun=get.anon_fit_fun(this) anon_fit_fun=this.FitStruct.(this.fit_name).anon_fit_fun; end %Grabs the correct fit parameters from FitStruct function fit_params=get.fit_params(this) fit_params=this.FitStruct.(this.fit_name).fit_params; end %Grabs the correct fit parameter names from FitStruct function fit_param_names=get.fit_param_names(this) fit_param_names=this.FitStruct.(this.fit_name).fit_param_names; end %Calculates the number of parameters in the fit function function n_params=get.n_params(this) n_params=length(this.fit_params); end %Generates a vector of x values for plotting function x_vec=get.x_vec(this) x_vec=linspace(min(this.Data.x),max(this.Data.x),1e3); end %Used when creating the UserGUI, finds the number of user fields. function n_user_fields=get.n_user_fields(this) n_user_fields=length(this.user_field_tags); end %Finds all the user field tags function user_field_tags=get.user_field_tags(this) user_field_tags=fieldnames(this.UserGui.Fields); end %Finds all the titles of the user field tags function user_field_names=get.user_field_names(this) user_field_names=cellfun(@(x) this.UserGui.Fields.(x).title,... this.user_field_tags,'UniformOutput',0); end %Finds all the values of the user fields function user_field_vals=get.user_field_vals(this) user_field_vals=cellfun(@(x) this.(x), this.user_field_tags); end %Generates a full path for saving function fullpath=get.fullpath(this) fullpath=[this.save_path,this.filename,'.txt']; end %Generates a base path for saving function save_path=get.save_path(this) save_path=createSessionPath(this.base_dir,this.session_name); end end %Set and get functions for dependent variables with SetAccess methods %Gets the base dir from the gui function base_dir=get.base_dir(this) if this.enable_gui base_dir=this.Gui.BaseDir.String; else base_dir=this.SaveInfo.base_dir; end end function set.base_dir(this,base_dir) if this.enable_gui this.Gui.BaseDir.String=base_dir; else this.SaveInfo.base_dir=base_dir; end end function session_name=get.session_name(this) if this.enable_gui session_name=this.Gui.SessionName.String; else session_name=this.SaveInfo.session_name; end end function set.session_name(this,session_name) if this.enable_gui this.Gui.SessionName.String=session_name; else this.SaveInfo.session_name=session_name; end end function filename=get.filename(this) if this.enable_gui filename=this.Gui.FileName.String; else filename=this.SaveInfo.filename; end end function set.filename(this,filename) if this.enable_gui this.Gui.FileName.String=filename; else this.SaveInfo.filename=filename; end end end end \ No newline at end of file diff --git a/@MyZiRingdown/MyZiRingdown.m b/@MyZiRingdown/MyZiRingdown.m index 89adefb..c9a4e2e 100644 --- a/@MyZiRingdown/MyZiRingdown.m +++ b/@MyZiRingdown/MyZiRingdown.m @@ -1,869 +1,861 @@ % Class for performing ringdown measurements of mechanical oscillators % using Zurich Instruments UHF or MF lock-in. % % Operation: sweep the driving tone (drive_osc) using the sweep module % in LabOne web user interface, when the magnitude of the demodulator % signal exceeds trig_threshold the driving tone is switched off and % the recording of demodulated signal is started, the signal is recorded % for the duration of record_time. % % Features: % % Adaptive measurement oscillator frequency % % Averaging % % Auto saving % % Auxiliary output signal: If enable_aux_out=true % then after a ringdown is started a sequence of pulses is applied % to the output consisting of itermittent on and off periods % starting from on. classdef MyZiRingdown < MyZiLi & MyDataSource properties (Access=public) % Ringdown is recorded if the signal in the triggering demodulation % channel exceeds this value trig_threshold=1e-3 % V % Duration of the recorded ringdown record_time=1 % (s) % If enable_acq is true, then the drive is on and the acquisition % of record is triggered when signal exceeds trig_threshold enable_acq=false % Auxiliary output signal during ringdown. enable_aux_out=false % If auxiliary output is applied % time during which the output is in aux_out_on_lev state aux_out_on_t=1 % (s) % time during which the output is in aux_out_off_lev state aux_out_off_t=1 % (s) aux_out_on_lev=1 % (V), output trigger on level aux_out_off_lev=0 % (V), output trigger off level % Average the trace over n points to reduce amount of stored data % while keeping the demodulator bandwidth large downsample_n=1 fft_length=128 auto_save=false % if all ringdowns should be automatically saved % In adaptive measurement oscillator mode the oscillator frequency % is continuously changed to follow the signal frequency during % ringdown acquisition. This helps against the oscillator frequency % drift. adaptive_meas_osc=false end % The properties which are read or set only once during the class % initialization properties (GetAccess=public, SetAccess={?MyClassParser}) % enumeration for demodulators, oscillators and output starts from 1 demod=1 % demodulator used for both triggering and measurement % Enumeration in the node structure starts from 0, so, for example, % the default path to the trigger demodulator refers to the % demodulator #1 demod_path='/dev4090/demods/0' drive_osc=1 meas_osc=2 % Signal input, integers above 1 correspond to main inputs, aux % input etc. (see the user interface for device-specific details) signal_in=1 drive_out=1 % signal output used for driving % Number of an auxiliary channel used for the output of triggering % signal, primarily intended to switch the measurement apparatus % off during a part of the ringdown and thus allow for free % evolution of the oscillator during that period. aux_out=1 % Poll duration of 1 ms practically means that ziDAQ('poll', ... % returns immediately with the data accumulated since the % previous function call. poll_duration=0.001 % s poll_timeout=50 % ms % Margin for adaptive oscillator frequency adjustment - oscillator % follows the signal if the dispersion of frequency in the % demodulator band is below ad_osc_margin times the demodulation % bandwidth (under the condition that adaptive_meas_osc=true) ad_osc_margin=0.1 end % Internal variables properties (GetAccess=public, SetAccess=protected) recording=false % true if a ringdown is being recorded % true if adaptive measurement oscillator mode is on and if the % measurement oscillator is actually actively following. ad_osc_following=false % Reference timestamp at the beginning of measurement record. % Stored as uint64. t0 elapsed_t=0 % Time elapsed since the last recording was started DemodSpectrum % MyTrace object to store FFT of the demodulator data end % Setting or reading the properties below automatically invokes % communication with the device properties (Dependent=true) drive_osc_freq meas_osc_freq drive_on % true when the dirive output is on % demodulator sampling rate (as transferred to the computer) demod_rate % The properties below are only used within the program to display % the information about device state. drive_amp % (V), peak-to-peak amplitude of the driving tone lowpass_order % low-pass filter order lowpass_bw % low-pass filter bandwidth end % Other dependent variables that are dont device properties properties (Dependent=true) % Downsample the measurement record to reduce the amount of data % while keeping the large demodulation bandwidth. % (samples/s), sampling rate of the trace after avraging downsampled_rate % number of the oscillator presently in use with the demodulator current_osc % true/false, true if the signal output from aux out is in on state aux_out_on % Provides public access to properties of private AvgTrace n_avg % number of ringdowns to be averaged avg_count % the average counter fft_rbw % resolution bandwidth of fft end % Keeping handle objects fully private is the only way to restrict set % access to their properties properties (Access=private) PollTimer AuxOutOffTimer % Timer responsible for switching the aux out off AuxOutOnTimer % Timer responsible for switching the aux out on % Demodulator samples z(t) stored to continuously calculate % spectrum, values of z are complex here, z=x+iy. % osc_freq is the demodulation frequency DemodRecord=struct('t',[],'z',[],'osc_freq',[]) AvgTrace % MyAvgTrace object used for averaging ringdowns end events NewDemodSample % New demodulator samples received RecordingStarted % Acquisition of a new trace triggered end methods (Access=public) %% Constructor and destructor function this = MyZiRingdown(dev_serial, varargin) this=this@MyZiLi(dev_serial); P=MyClassParser(this); addRequired(P, dev_serial, @ischar) % Poll timer period addParameter(P,'poll_period',0.042,@isnumeric); processInputs(P, this, dev_serial, varargin{:}); % Create and configure trace objects % Trace is inherited from the superclass this.Trace=MyTrace(... 'name_x','Time',... 'unit_x','s',... 'name_y','Magnitude r',... 'unit_y','V'); this.DemodSpectrum=MyTrace(... 'name_x','Frequency',... 'unit_x','Hz',... 'name_y','PSD',... 'unit_y','V^2/Hz'); this.AvgTrace=MyAvgTrace(); % Set up the poll timer. Using a timer for anyncronous % data readout allows to use the wait time for execution % of other programs. % Fixed spacing is preferred as it is the most robust mode of % operation when keeping the intervals between callbacks % precisely defined is not the biggest concern. this.PollTimer=timer(... 'ExecutionMode','fixedSpacing',... 'Period',P.Results.poll_period,... 'TimerFcn',@(~,~)pollTimerCallback(this)); % Aux out timers use fixedRate mode for more precise timing. % The two timers are executed periodically with a time lag. % The first timer switches the auxiliary output off this.AuxOutOffTimer=timer(... 'ExecutionMode','fixedRate',... 'TimerFcn',@(~,~)auxOutOffTimerCallback(this)); % The second timer switches the auxiliary output on this.AuxOutOnTimer=timer(... 'ExecutionMode','fixedRate',... 'TimerFcn',@(~,~)auxOutOnTimerCallback(this)); this.demod_path = sprintf('/%s/demods/%i', this.dev_id, ... this.demod-1); end function delete(this) % delete function should never throw errors, so protect % statements with try-catch try stopPoll(this) catch warning(['Could not usubscribe from the demodulator ', ... 'or stop the poll timer.']) end % Delete timers to prevent them from running indefinitely in % the case of program crash try delete(this.PollTimer) catch warning('Could not delete the poll timer.') end try stop(this.AuxOutOffTimer); delete(this.AuxOutOffTimer); catch warning('Could not stop and delete AuxOutOff timer.') end try stop(this.AuxOutOnTimer); delete(this.AuxOutOnTimer); catch warning('Could not stop and delete AuxOutOn timer.') end end %% Other methods function startPoll(this) % Configure the oscillators, demodulator and driving output % -1 accounts for the difference in enumeration conventions % in the software names (starting from 1) and node numbers % (starting from 0). % First, update the demodulator path this.demod_path = sprintf('/%s/demods/%i', ... this.dev_id, this.demod-1); % Set the data transfer rate so that it satisfies the Nyquist % criterion (>x2 the bandwidth of interest) this.demod_rate=4*this.lowpass_bw; % Configure the demodulator. Signal input: ziDAQ('setInt', ... [this.demod_path,'/adcselect'], this.signal_in-1); % Oscillator: ziDAQ('setInt', ... [this.demod_path,'/oscselect'], this.drive_osc-1); % Enable data transfer from the demodulator to the computer ziDAQ('setInt', [this.demod_path,'/enable'], 1); % Configure the signal output - disable all the oscillator % contributions excluding the driving tone path = sprintf('/%s/sigouts/%i/enables/*', ... this.dev_id, this.drive_out-1); ziDAQ('setInt', path, 0); path = sprintf('/%s/sigouts/%i/enables/%i', ... this.dev_id, this.drive_out-1, this.drive_osc-1); ziDAQ('setInt', path, 1); % Enable output this.drive_on=true; % By convention, we start form 'enable_acq=false' state this.enable_acq=false; % Configure the auxiliary trigger output - put it in the manual % mode so it does not output demodulator readings path=sprintf('/%s/auxouts/%i/outputselect', ... this.dev_id, this.aux_out-1); ziDAQ('setInt', path, -1); % The convention is that aux out is on by default this.aux_out_on=true; % Subscribe to continuously receive samples from the % demodulator. Samples accumulated between timer callbacks % will be read out using ziDAQ('poll', ... ziDAQ('subscribe',[this.demod_path,'/sample']); % Start continuous polling start(this.PollTimer) end function stopPoll(this) stop(this.PollTimer) ziDAQ('unsubscribe',[this.demod_path,'/sample']); end % Main function that polls data from the device demodulator function pollTimerCallback(this) % ziDAQ('poll', ... with short poll_duration returns % immediately with the data accumulated since the last timer % callback Data = ziDAQ('poll', this.poll_duration, this.poll_timeout); if ziCheckPathInData(Data, [this.demod_path,'/sample']) % Demodulator returns data DemodSample= ... Data.(this.dev_id).demods(this.demod).sample; % Append new samples to the record and recalculate spectrum appendSamplesToBuff(this, DemodSample); calcfft(this); if this.recording % If recording is under way, append the new samples to % the trace rec_finished = appendSamplesToTrace(this, DemodSample); % Recording can be manually stopped by setting % enable_acq=false if ~this.enable_acq rec_finished=true; end % Update elapsed time this.elapsed_t=this.Trace.x(end); % If the adaptive measurement frequency mode is on, % update the measurement oscillator frequency. % Make sure that the demodulator record actually % contains signal by comparing the dispersion of % frequency to demodulator bandwidth. if this.adaptive_meas_osc [df_avg, df_dev]=calcfreq(this); if df_dev < this.ad_osc_margin*this.lowpass_bw this.meas_osc_freq=df_avg; % Change indicator this.ad_osc_following=true; else this.ad_osc_following=false; end else this.ad_osc_following=false; end else r=sqrt(DemodSample.x.^2+DemodSample.y.^2); if this.enable_acq && max(r)>this.trig_threshold % Start acquisition of a new trace if the maximum % of the signal exceeds threshold this.recording=true; % Find index at which the threshold was % exceeded ind0=find(r>this.trig_threshold,1,'first'); this.t0=DemodSample.timestamp(ind0); this.elapsed_t=0; % Switch the drive off this.drive_on=false; % Set the measurement oscillator frequency to be % the frequency at which triggering occurred this.meas_osc_freq=this.drive_osc_freq; % Switch the oscillator this.current_osc=this.meas_osc; - % Clear the buffer on ZI data server from existing - % demodulator samples, as these samples were - % recorded with drive on - ziDAQ('poll',this.poll_duration,this.poll_timeout); - % Optionally start the auxiliary output timers if this.enable_aux_out % Configure measurement periods and delays T=this.aux_out_on_t+this.aux_out_off_t; this.AuxOutOffTimer.Period=T; this.AuxOutOnTimer.Period=T; this.AuxOutOffTimer.startDelay=... this.aux_out_on_t; this.AuxOutOnTimer.startDelay=T; % Start timers start(this.AuxOutOffTimer) start(this.AuxOutOnTimer) end % Clear trace and append new data starting from the - % index, at which triggering occurred + % index, at which triggering occurred. + % Theoretically, a record can be finished with + % this one portion if the record time is set small. clearData(this.Trace); - % Append the first portion of samples to the - % record, starting from the index at which - % triggering occurred. Theoretically, a record can - % be finished with this one portion if the record - % time is set small. rec_finished = ... appendSamplesToTrace(this, DemodSample, ind0); notify(this, 'RecordingStarted'); else rec_finished=false; end % Indicator for adaptive measurement is off, since % recording is not under way this.ad_osc_following=false; end notify(this,'NewDemodSample'); % Stop recording if a record was completed if rec_finished % stop recording this.recording=false; this.ad_osc_following=false; % Stop auxiliary timers stop(this.AuxOutOffTimer); stop(this.AuxOutOnTimer); % Return the drive and aux out to the default state this.aux_out_on=true; this.current_osc=this.drive_osc; this.drive_on=true; % Downsample the trace to reduce the amount of data downsample(this.Trace, this.downsample_n, 'avg'); % Do trace averaging. If the new data length is not of % the same size as the length of the existing data % (which should happen only when the record period was % changed during recording or when recording was % manually stopped), truncate to the minimum length if ~isempty(this.AvgTrace) && ... (length(this.AvgTrace.y)~=length(this.Trace.y)) l=min(length(this.AvgTrace.y),length(this.Trace.y)); this.AvgTrace.y=this.AvgTrace.y(1:l); this.AvgTrace.x=this.AvgTrace.x(1:l); this.Trace.y=this.Trace.y(1:l); this.Trace.x=this.Trace.x(1:l); disp('Ringdown record was truncated') end avg_compl=addAverage(this.AvgTrace, this.Trace); % Trigger NewData if this.n_avg>1 end_str=sprintf('_%i', this.AvgTrace.avg_count); else end_str=''; end triggerNewData(this, 'save', this.auto_save, ... 'filename_ending', end_str); % If the ringdown averaging is complete, disable % further triggering to exclude data overwriting if avg_compl this.enable_acq=false; if this.n_avg>1 end_str='_avg'; % Trigger one more time to transfer the average % trace. % A new measurement header is not necessary % as the delay since the last triggering is % minimum. triggerNewData(this, ... 'Trace', copy(this.AvgTrace), ... 'save', this.auto_save, ... 'filename_ending', end_str); end end end end end % Append timestamps vs r=sqrt(x^2+y^2) to the measurement record. % Starting index can be supplied as varargin. % The output variable tells if the record is finished. function isfin = appendSamplesToTrace(this, DemodSample, varargin) if isempty(varargin) startind=1; else startind=varargin{1}; end r=sqrt(DemodSample.x(startind:end).^2 + ... DemodSample.y(startind:end).^2); % Subtract the reference time, convert timestamps to seconds ts=double(DemodSample.timestamp(startind:end) -... this.t0)/this.clockbase; % Check if recording should be stopped isfin=(ts(end)>=this.record_time); if isfin % Remove excess data points from the new data ind=(tsflen this.DemodRecord.t = this.DemodRecord.t(end-flen+1:end); this.DemodRecord.z = this.DemodRecord.z(end-flen+1:end); this.DemodRecord.osc_freq = ... this.DemodRecord.osc_freq(end-flen+1:end); end end function calcfft(this) flen=min(this.fft_length, length(this.DemodRecord.t)); [freq, spectr]=xyFourier( ... this.DemodRecord.t(end-flen+1:end), ... this.DemodRecord.z(end-flen+1:end)); this.DemodSpectrum.x=freq; this.DemodSpectrum.y=abs(spectr).^2; end % Calculate the average frequency and dispersion of the demodulator % signal function [f_avg, f_dev]=calcfreq(this) if ~isempty(this.DemodSpectrum) norm=sum(this.DemodSpectrum.y); % Calculate the center frequency of the spectrum f_avg=dot(this.DemodSpectrum.x, this.DemodSpectrum.y)/norm; f_dev=sqrt(dot(this.DemodSpectrum.x.^2, ... this.DemodSpectrum.y)/norm-f_avg^2); % Shift the FFT center by the demodulation frequency to % output absolute value f_avg=f_avg+mean(this.DemodRecord.osc_freq); else f_avg=[]; f_dev=[]; end end % Provide restricted access to private AvgTrace function resetAveraging(this) % Clear data and reset counter clearData(this.AvgTrace); notify(this,'NewSetting'); end function auxOutOffTimerCallback(this) this.aux_out_on=false; end function auxOutOnTimerCallback(this) this.aux_out_on=true; end %% measurement header functionality function Hdr=readHeader(this) Hdr=readHeader@MyZiLi(this); % Demodulator parameters addObjProp(Hdr, this, 'demod', 'comment', ... 'Number of the demodulator in use (starting from 1)'); addObjProp(Hdr, this, 'demod_rate', 'comment', ... '(samples/s), demodulator data transfer rate'); addObjProp(Hdr, this, 'lowpass_order', 'comment', ... 'Order of the demodulator low-pass filter'); addObjProp(Hdr, this, 'lowpass_bw', 'comment', ... ['(Hz), 3 dB bandwidth of the demodulator low-pass ', ... 'filter']); addObjProp(Hdr, this, 'meas_osc', 'comment', ... 'Measurement oscillator number'); addObjProp(Hdr, this, 'meas_osc_freq', 'comment', '(Hz)'); % Signal input addObjProp(Hdr, this, 'signal_in', 'comment', ... 'Singnal input number'); % Drive parameters addObjProp(Hdr, this, 'drive_out', 'comment', ... 'Driving output number'); addObjProp(Hdr, this, 'drive_osc', 'comment', ... 'Swept oscillator number'); addObjProp(Hdr, this, 'drive_amp', 'comment', ... '(V) peak to peak'); % Parameters of the auxiliary output addObjProp(Hdr, this, 'aux_out', 'comment', ... 'Auxiliary output number'); addObjProp(Hdr, this, 'enable_aux_out', 'comment', ... 'Auxiliary output is applied during ringdown'); addObjProp(Hdr, this, 'aux_out_on_lev', 'comment', '(V)'); addObjProp(Hdr, this, 'aux_out_off_lev', 'comment', '(V)'); addObjProp(Hdr, this, 'aux_out_on_t', 'comment', '(s)'); addObjProp(Hdr, this, 'aux_out_off_t', 'comment', '(s)'); % Software parameters addObjProp(Hdr, this, 'trig_threshold', 'comment', ... '(V), threshold for starting a ringdown record'); addObjProp(Hdr, this, 'record_time', 'comment', '(s)'); addObjProp(Hdr, this, 'n_avg', 'comment', ... 'Number of ringdowns to be averaged'); addObjProp(Hdr, this, 'downsampled_rate', 'comment', ... ['(samples/s), rate to which a ringown trace is ', ... 'downsampled with averaging after acquisition']); addObjProp(Hdr, this, 'auto_save', 'comment', '(s)'); % Adaptive measurement oscillator addObjProp(Hdr, this, 'adaptive_meas_osc', 'comment', ... ['If true the measurement oscillator frequency is ', ... 'adjusted during ringdown']); addObjProp(Hdr, this, 'ad_osc_margin'); addObjProp(Hdr, this, 'fft_length', 'comment', '(points)'); % Timer poll parameters addParam(Hdr, this.name, 'poll_period', ... this.PollTimer.Period, 'comment', '(s)'); addObjProp(Hdr, this, 'poll_duration', 'comment', '(s)'); addObjProp(Hdr, this, 'poll_timeout', 'comment', '(ms)'); end end %% Set and get methods. methods function freq=get.drive_osc_freq(this) path=sprintf('/%s/oscs/%i/freq', this.dev_id, this.drive_osc-1); freq=ziDAQ('getDouble', path); end function set.drive_osc_freq(this, val) assert(isfloat(val), ... 'Oscillator frequency must be a floating point number') path=sprintf('/%s/oscs/%i/freq', this.dev_id, this.drive_osc-1); ziDAQ('setDouble', path, val); notify(this,'NewSetting'); end function freq=get.meas_osc_freq(this) path=sprintf('/%s/oscs/%i/freq', this.dev_id, this.meas_osc-1); freq=ziDAQ('getDouble', path); end function set.meas_osc_freq(this, val) assert(isfloat(val), ... 'Oscillator frequency must be a floating point number') path=sprintf('/%s/oscs/%i/freq', this.dev_id, this.meas_osc-1); ziDAQ('setDouble', path, val); notify(this,'NewSetting'); end function set.drive_on(this, val) path=sprintf('/%s/sigouts/%i/on',this.dev_id,this.drive_out-1); % Use double() to convert from logical type ziDAQ('setInt', path, double(val)); notify(this,'NewSetting'); end function bool=get.drive_on(this) path=sprintf('/%s/sigouts/%i/on',this.dev_id,this.drive_out-1); bool=logical(ziDAQ('getInt', path)); end function set.current_osc(this, val) assert((val==this.drive_osc) || (val==this.meas_osc), ... ['The number of current oscillator must be that of ', ... 'the drive or measurement oscillator, not ', num2str(val)]) ziDAQ('setInt', [this.demod_path,'/oscselect'], val-1); notify(this,'NewSetting') end function osc_num=get.current_osc(this) osc_num=double(ziDAQ('getInt', ... [this.demod_path,'/oscselect']))+1; end function amp=get.drive_amp(this) path=sprintf('/%s/sigouts/%i/amplitudes/%i', ... this.dev_id, this.drive_out-1, this.drive_osc-1); amp=ziDAQ('getDouble', path); end function set.drive_amp(this, val) path=sprintf('/%s/sigouts/%i/amplitudes/%i', ... this.dev_id, this.drive_out-1, this.drive_osc-1); ziDAQ('setDouble', path, val); notify(this,'NewSetting'); end function set.lowpass_order(this, val) assert(any(val==[1,2,3,4,5,6,7,8]), ['Low-pass filter ', ... 'order must be an integer between 1 and 8']) ziDAQ('setInt', [this.demod_path,'/order'], val); notify(this,'NewSetting'); end function n=get.lowpass_order(this) n=ziDAQ('getInt', [this.demod_path,'/order']); end function bw=get.lowpass_bw(this) tc=ziDAQ('getDouble', [this.demod_path,'/timeconstant']); bw=ziBW2TC(tc, this.lowpass_order); end function set.lowpass_bw(this, val) tc=ziBW2TC(val, this.lowpass_order); ziDAQ('setDouble', [this.demod_path,'/timeconstant'], tc); notify(this,'NewSetting'); end function rate=get.demod_rate(this) rate=ziDAQ('getDouble', [this.demod_path,'/rate']); end function set.demod_rate(this, val) ziDAQ('setDouble', [this.demod_path,'/rate'], val); notify(this,'NewSetting'); end function set.downsample_n(this, val) n=round(val); assert(n>=1, ['Number of points for trace averaging must ', ... 'be greater than 1']) this.downsample_n=n; notify(this,'NewSetting'); end function set.aux_out_on(this, bool) path=sprintf('/%s/auxouts/%i/offset', ... this.dev_id, this.aux_out-1); if bool out_offset=this.aux_out_on_lev; else out_offset=this.aux_out_off_lev; end ziDAQ('setDouble', path, out_offset); notify(this,'NewSetting'); end function bool=get.aux_out_on(this) path=sprintf('/%s/auxouts/%i/offset', ... this.dev_id, this.aux_out-1); val=ziDAQ('getDouble', path); % Signal from the auxiliary output is continuous, we make the % binary decision about the output state depending on if % the signal is closer to the ON or OFF level bool=(abs(val-this.aux_out_on_lev) < ... abs(val-this.aux_out_off_lev)); end function set.downsampled_rate(this, val) dr=this.demod_rate; if val>dr % Downsampled rate should not exceed the data transfer rate val=dr; end % Round so that the averaging is done over an integer number of % points this.downsample_n=round(dr/val); notify(this,'NewSetting'); end function val=get.downsampled_rate(this) val=this.demod_rate/this.downsample_n; end function set.fft_length(this, val) if val<1 val=1; end % Round val to the nearest 2^n to make the calculation of % Fourier transform efficient n=round(log2(val)); this.fft_length=2^n; notify(this,'NewSetting'); end function val=get.fft_rbw(this) val=this.demod_rate/this.fft_length; end function set.fft_rbw(this, val) assert(val>0,'FFT resolution bandwidth must be greater than 0') % Rounding of fft_length to the nearest integer is handled by % its own set method this.fft_length=this.demod_rate/val; notify(this,'NewSetting'); end function set.n_avg(this, val) this.AvgTrace.n_avg=val; notify(this,'NewSetting'); end function val=get.n_avg(this) val=this.AvgTrace.n_avg; end function val=get.avg_count(this) val=this.AvgTrace.avg_count; end function set.aux_out_on_t(this, val) assert(val>0.001, ... 'Aux out on time must be greater than 0.001 s.') this.aux_out_on_t=val; notify(this,'NewSetting'); end function set.aux_out_off_t(this, val) assert(val>0.001, ... 'Aux out off time must be greater than 0.001 s.') this.aux_out_off_t=val; notify(this,'NewSetting'); end function set.enable_acq(this, val) this.enable_acq=logical(val); notify(this,'NewSetting'); end end end