diff --git a/Measurement and analysis routines/FreqNoiseCal/MyFreqNoiseCal.m b/Measurement and analysis routines/FreqNoiseCal/MyFreqNoiseCal.m index 9516c6e..ce5c3fc 100644 --- a/Measurement and analysis routines/FreqNoiseCal/MyFreqNoiseCal.m +++ b/Measurement and analysis routines/FreqNoiseCal/MyFreqNoiseCal.m @@ -1,202 +1,207 @@ -% Routine that calibrates noise spectrum in frequency units using a -% calibration tone with known phase modulation depth beta +% Routine for the calibration of frequency noise spectrum using a +% calibration tone with known phase modulation depth, beta, defined as +% follows: +% +% E_0(t) = A*Exp(-i\beta \cos(\Omega_{cal} t)) classdef MyFreqNoiseCal < MyAnalysisRoutine properties (Access = public, SetObservable = true) % Thermomechanical spectrum Data MyTrace % Spectrum, calibrated in frequency units. It can be % S_\omega % S_f = S_\omega/(2pi)^2 % sqrt(S_f) FreqSpectrum MyTrace % The type of calculated frequency spectrum, % S_\omega, S_f or sqrt(S_f) - spectrum_type = 'sqrt(S_f)' + spectrum_type = 'S_f' % Cursors for the selection of calibration tone CalCursors MyCursor cal_range = [0, 0] - % Phase modulation depth of calibration tone + % Phase modulation depth of the calibration tone, supplied + % externally beta = 0.1 end properties (GetAccess = public, SetAccess = protected, ... SetObservable = true) Axes Gui + % Frequency of the calibration tone, found in the cal_range cal_freq % Conversion factor between S_V and S_\omega defined such that % S_\omega = cf*S_V cf = 1 end methods (Access = public) function this = MyFreqNoiseCal(varargin) p = inputParser(); addParameter(p, 'Data', MyTrace()); addParameter(p, 'Axes', [], @isaxes); addParameter(p, 'enable_gui', true, @islogical); parse(p, varargin{:}); this.Data = p.Results.Data; this.Axes = p.Results.Axes; this.FreqSpectrum = MyTrace(); if ~isempty(this.Axes) % Add two sets of vertical cursors for the selection of % integration ranges xlim = this.Axes.XLim; x1 = xlim(1)+0.4*(xlim(2)-xlim(1)); x2 = xlim(1)+0.45*(xlim(2)-xlim(1)); this.CalCursors = ... [MyCursor(this.Axes, ... 'orientation', 'vertical', 'position', x1, ... 'Label','Cal 1', 'Color', [0, 0, 0.6]), ... MyCursor(this.Axes, ... 'orientation', 'vertical', 'position', x2, ... 'Label','Cal 2', 'Color', [0, 0, 0.6])]; end % Gui is created right before the construction of object % is over if p.Results.enable_gui this.Gui = GuiFreqNoiseCal(this); end end function delete(this) if ~isempty(this.CalCursors) delete(this.CalCursors); end end - % Calculates the frequency noise spectrum + % Calculates cf function calcConvFactor(this) if isempty(this.Data) || isDataEmpty(this.Data) warning('Data is empty'); return end cr = this.cal_range; % Find the frequency of calibration tone ind = (this.Data.x>cr(1) & this.Data.x<=cr(2)); freq = sum(this.Data.x(ind) .* this.Data.y(ind))/ ... sum(this.Data.y(ind)); this.cal_freq = freq; % Calculate area under the calibration tone peak area = integrate(this.Data, cr(1), cr(2)); % Average square of frequency excursions due to calibration % tone vSqCt = this.beta^2*(2*pi*freq)^2/2; this.cf = vSqCt/area; end function convertSpectrum(this) if isempty(this.Data) || isDataEmpty(this.Data) warning('Data is empty'); return end this.FreqSpectrum.x = this.Data.x; this.FreqSpectrum.name_x = this.Data.name_x; this.FreqSpectrum.unit_x = 'Hz'; switch this.spectrum_type case 'S_\omega' this.FreqSpectrum.name_y = '$S_{\omega}$'; this.FreqSpectrum.unit_y = 'rad$^2$/Hz'; this.FreqSpectrum.y = this.Data.y*this.cf; case 'S_f' this.FreqSpectrum.name_y = '$S_{\omega}/(2\pi)^2$'; this.FreqSpectrum.unit_y = 'Hz$^2$/Hz'; this.FreqSpectrum.y = ... this.Data.y*this.cf/(2*pi)^2; case 'sqrt(S_f)' this.FreqSpectrum.name_y ='$\sqrt{S_{\omega}}/(2\pi)$'; this.FreqSpectrum.unit_y ='Hz/$\sqrt{\mathrm{Hz}}$'; this.FreqSpectrum.y = ... sqrt(this.Data.y*this.cf)/(2*pi); otherwise error(['Unknown frequency spectrum type ' ... this.spectrum_type]) end % Update metadata this.FreqSpectrum.UserMetadata = createMetadata(this); triggerNewAnalysisTrace(this,'Trace',copy(this.FreqSpectrum)); end end methods (Access = protected) function Mdt = createMetadata(this) Mdt = MyMetadata('title', 'CalibrationParameters'); if ~isempty(this.Data.file_name) % If data file has name, indicate it addParam(Mdt, 'source', this.Data.file_name, ... 'comment', 'File containing raw data'); end addParam(Mdt, 'spectrum_type', this.spectrum_type); addParam(Mdt, 'beta', this.beta, ... 'comment', 'Phase modulation depth'); addParam(Mdt, 'cal_freq', this.cal_freq, ... 'comment', ['Calibration tone frequency (' ... this.Data.unit_x ')']); addParam(Mdt, 'cf', this.cf, ... 'comment', 'Conversion factor S_\omega = cf*S_V'); end end % Set and get methods methods % Get the integration range for the calibration tone function val = get.cal_range(this) if ~isempty(this.CalCursors) && all(isvalid(this.CalCursors)) % If cursors exist, return the range between the % cursors xmin = min(this.CalCursors.value); xmax = max(this.CalCursors.value); val = [xmin, xmax]; else % Otherwise the value can be set programmatically val = this.cal_range; end end function set.spectrum_type(this, val) val_list = {'S_\omega', 'S_f', 'sqrt(S_f)'}; assert(ismember(val, val_list), ['The value of ' ... 'spectrum_type must be one of the following: ' ... var2str(val_list)]) this.spectrum_type = val; end end end diff --git a/Measurement and analysis routines/PhaseModCal/MyPhaseModCal.m b/Measurement and analysis routines/PhaseModCal/MyPhaseModCal.m index 2f18af6..81b8b44 100644 --- a/Measurement and analysis routines/PhaseModCal/MyPhaseModCal.m +++ b/Measurement and analysis routines/PhaseModCal/MyPhaseModCal.m @@ -1,258 +1,257 @@ % Routine for the calibration of beta-factor, characterizing the phase % modulation of light, using heterodyne signal spectrum. % -% Beta in the following expression for phase-modulated complex amplitude of -% light: +% Beta is defined in the following expression for phase-modulated complex +% amplitude of light: % -% A*Exp(-i\beta \cos(\Omega_{cal} t)) +% E_0(t) = A*Exp(-i\beta \cos(\Omega_{cal} t)) classdef MyPhaseModCal < MyAnalysisRoutine properties (Access = public, SetObservable = true) Data MyTrace % Minimum thereshold for peak search. If MinHeightCursor exists, it % has priority over the programmatically set value. min_peak_height double mod_freq = 0 % modulation frequency (Hz) end properties (Access = public, Dependent = true, SetObservable = true) enable_cursor end properties (GetAccess = public, SetAccess = protected, ... SetObservable = true) Axes Gui MinHeightCursor MyCursor beta = 0 % Phase modulation depth end properties (Access = protected) % Line that displays the positions of peaks found PlottedPeaks end methods (Access = public) function this = MyPhaseModCal(varargin) p = inputParser(); addParameter(p, 'Data', MyTrace()); addParameter(p, 'Axes', [], @isaxes); addParameter(p, 'enable_cursor', true, @islogical); addParameter(p, 'enable_gui', true, @islogical); parse(p, varargin{:}); this.Data = p.Results.Data; this.Axes = p.Results.Axes; if ~isempty(this.Axes) && isvalid(this.Axes) ylim = this.Axes.YLim; pos = min(ylim(1)+0.1*(ylim(2)-ylim(1)), 10*ylim(1)); this.MinHeightCursor = MyCursor(this.Axes, ... 'orientation', 'horizontal', ... 'position', pos, ... 'Label', 'Peak threshold', ... 'Color', [0.6, 0, 0]); this.min_peak_height = pos; this.enable_cursor = p.Results.enable_cursor; else this.min_peak_height = 1e-12; end % Gui is created right before the construction of object % is over if p.Results.enable_gui this.Gui = GuiPhaseModCal(this); end end function delete(this) if ~isempty(this.PlottedPeaks) delete(this.PlottedPeaks) end if ~isempty(this.MinHeightCursor) delete(this.MinHeightCursor) end end % Calculate the depth of phase modulation from the hights of peaks % in the spectrum function calcBeta(this) min_y = this.min_peak_height; % Find peaks above the given threshold % Returned values: [y, x, widths, prominences] [peak_y, peak_x, peak_w, ~] = findpeaks( ... this.Data.y, this.Data.x, 'MinPeakHeight', min_y); n_peaks = length(peak_y); assert(n_peaks >= 3, ['Less than 3 peaks are found in '... 'the data with given threshold (' num2str(min_y) '). ' ... 'Phase modulation depth cannot be calculated.']) % Find the central peak, which is not necessarily the highest mean_freq = sum(peak_x.*peak_y)/sum(peak_y); [~, cent_ind] = min(abs(peak_x-mean_freq)); - % Take the integration width to a few times the width of the + % Take the integration width to be a few times the width of the % central peak. int_w = 6*peak_w(cent_ind); % Check if the peaks are rougly equally spaced harmonics, if % not, use the pre-specified value of modulation frequency to % select the right peaks. peak_x_diff = peak_x(2:n_peaks)-peak_x(1:n_peaks-1); mod_f = mean(peak_x_diff); % Specify the tolerance to the mismatch of frequencies between % the found modulation peaks mod_peak_tol = 0.1; if max(abs(peak_x_diff-mod_f))/mod_f > mod_peak_tol % Try using the approximate value of modulation frequency % that can be specified by the user. disp(['Distances between the found peaks are not ' ... 'regular, will use the pre-defined value of ' ... 'modulation frequency to post select peaks.']); if isempty(this.mod_freq) || this.mod_freq<=0 % Prompt user to specify approximate modulation % frequency. Show warning dialog if running in a gui % mode or output warning in the command line otherwise. if ~isempty(this.Gui) Wd = warndlg(['Cannot identify modulation ' ... 'sidebands automatically. Please input ' ... 'an approximate value of modulation ' ... 'frequency and try again.'], 'Warning'); centerFigure(Wd); else warning(['An approximate value on modulation ' ... 'frequency must be specified by setting ' ... 'mod_freq property. Please specify the ' ... 'frequency and try again.']); end return end mod_f = this.mod_freq; end % Delete the peaks that do not appear at the expected % frequencies of harmonics of the modulation frequency scaled_peak_x = (peak_x - peak_x(cent_ind))/mod_f; sb_ind = cent_ind; for i = 1:ceil(n_peaks/2) % Iterate over the sideband index i and find pairs of % sidebands [err_p, ind_p] = min(abs(scaled_peak_x - i)); [err_m, ind_m] = min(abs(scaled_peak_x + i)); if (err_p/i < mod_peak_tol) && (err_m/i < mod_peak_tol) % Add the found indices to the list of sideband % peak indices sb_ind = [ind_m, sb_ind, ind_p]; %#ok else break end end % Out of all peaks select sideband peaks that appear in pairs peak_y = peak_y(sb_ind); peak_x = peak_x(sb_ind); n_peaks = length(peak_y); assert(n_peaks >= 3, ['Less than 3 peaks are found. ' ... 'Phase modulation depth cannot be calculated.']) % Re-calculate the modulation frequency mod_f = (peak_x(end)-peak_x(1))/(n_peaks-1); % Display the found peaks if ~isempty(this.Axes) && isvalid(this.Axes) if ~isempty(this.PlottedPeaks)&&isvalid(this.PlottedPeaks) set(this.PlottedPeaks,'XData',peak_x,'YData',peak_y); else this.PlottedPeaks = line(this.Axes, ... 'XData', peak_x, 'YData', peak_y, 'Color', 'r', ... 'LineStyle', 'none', 'Marker', 'o'); end end % Calculate areas under the peaks peak_int = zeros(1, n_peaks); for i = 1:n_peaks peak_int(i) = integrate(this.Data, peak_x(i)-int_w/2, ... peak_x(i)+int_w/2); end % Scale by the maximum area for better fit convergence peak_int = peak_int/max(peak_int); % Find beta value by fitting Ft = fittype('a*besselj(n, beta)^2', 'independent', 'n', ... 'coefficients', {'a', 'beta'}); Opts = fitoptions('Method', 'NonLinearLeastSquares',... 'StartPoint', [1, 0.1],... 'MaxFunEvals', 2000,... 'MaxIter', 2000,... 'TolFun', 1e-10,... 'TolX', 1e-10); peak_ord = -floor(n_peaks/2):floor(n_peaks/2); FitResult = fit(peak_ord(:), peak_int(:), Ft, Opts); % Store the result in class variables this.beta = abs(FitResult.beta); this.mod_freq = mod_f; end function clearPeakDisp(this) if ~isempty(this.PlottedPeaks) delete(this.PlottedPeaks) end end end % Set and get methods methods function set.enable_cursor(this, val) if ~isempty(this.MinHeightCursor) && ... isvalid(this.MinHeightCursor) this.MinHeightCursor.Line.Visible = val; end end function val = get.enable_cursor(this) if ~isempty(this.MinHeightCursor) && ... isvalid(this.MinHeightCursor) val = strcmpi(this.MinHeightCursor.Line.Visible, 'on'); else val = false; end end function val = get.min_peak_height(this) - if ~isempty(this.MinHeightCursor) && ... - isvalid(this.MinHeightCursor) + if this.enable_cursor val = this.MinHeightCursor.value; else val = this.min_peak_height; end end end end diff --git a/Resources/Analysis routines/beta_def.png b/Measurement and analysis routines/PhaseModCal/beta_def.png similarity index 100% rename from Resources/Analysis routines/beta_def.png rename to Measurement and analysis routines/PhaseModCal/beta_def.png diff --git a/Resources/Analysis routines/g0_cal.png b/Measurement and analysis routines/g0Cal/g0_cal.png similarity index 100% rename from Resources/Analysis routines/g0_cal.png rename to Measurement and analysis routines/g0Cal/g0_cal.png