From: W. Trevor King Date: Tue, 16 Jun 2009 17:05:34 +0000 (-0400) Subject: Restored linear-fitting option to surface bumps. X-Git-Tag: 0.4~2 X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=6369a5dcb5592450db8a5fa5f323c69fabb7ae2b;p=calibcant.git Restored linear-fitting option to surface bumps. Now you don't have to use quadratic fitting if you don't want to. This mitigates problems associated with poorly defined contact-kinks, since it's not a good idea to focus on the post-kink slope if you don't have a good idea of where the kink actually is. The limited_linear model retains the non-contact and high-voltage-rail flat-line portions, with a linear contact regime between the two. Additional minor changes include: * Average deflection printout from vib_analyze in TEXT_VERBOSE mode. To make it easier to demonstrate that variance increases as the mean deflection deviates further from zero. * Corrected usage string in bump_analyze.py which had been out of date before. --- diff --git a/calibcant/bump_analyze.py b/calibcant/bump_analyze.py index 21e13cb..3783936 100755 --- a/calibcant/bump_analyze.py +++ b/calibcant/bump_analyze.py @@ -99,7 +99,7 @@ def bump_analyze(data, **kwargs) : With the current implementation, the data is regressed in DAC/ADC bits and THEN converted, so we're assuming that both conversions are LINEAR. - if they aren't, rewrite to convert before the regression. + If they aren't, rewrite to convert before the regression. """ bump_fit_kwargs,slope_bitspbit2Vpnm_kwargs = \ bump_analyze._splitargs(bump_analyze, kwargs) @@ -108,27 +108,26 @@ def bump_analyze(data, **kwargs) : **bump_fit_kwargs) return slope_bitspbit2Vpnm(Vphoto2Vzp_out_bit, **slope_bitspbit2Vpnm_kwargs) -def limited_quadratic(x, params): +def limited_linear(x, params): """ Model the bump as: flat region (off-surface) - quadratic region (in-contact) + linear region (in-contact) flat region (high-voltage-rail) Parameters: x_contact (x value for the surface-contact kink) y_contact (y value for the surface-contact kink) slope (dy/dx at the surface-contact kink) - quad (d**2 y / dx**2, allow decreasing sensitivity with increased x) """ high_voltage_rail = 2**16 - 1 # bits - x_contact,y_contact,slope,quad = params - y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact + x_contact,y_contact,slope = params + y = slope*(x-x_contact) + y_contact y = numpy.clip(y, y_contact, high_voltage_rail) return y -def limited_quadratic_param_guess(x, y) : +def limited_linear_param_guess(x, y) : """ - Guess rough parameters for a limited_quadratic model. Assumes the + Guess rough parameters for a limited_linear model. Assumes the bump approaches (raising the deflection as it does so) first. Retracting after the approach is optional. Approximates the contact position and an on-surface (high) position by finding first crossings @@ -149,6 +148,44 @@ def limited_quadratic_param_guess(x, y) : x_contact = float(x[i_low]) x_high = float(x[i_high]) slope = (y_high - y_contact) / (x_high - x_contact) + return (x_contact, y_contact, slope) + +def limited_linear_sensitivity(params): + """ + Return the estimated sensitivity to small deflections according to + limited_linear fit parameters. + """ + slope = params[2] + return slope + +def limited_quadratic(x, params): + """ + Model the bump as: + flat region (off-surface) + quadratic region (in-contact) + flat region (high-voltage-rail) + Parameters: + x_contact (x value for the surface-contact kink) + y_contact (y value for the surface-contact kink) + slope (dy/dx at the surface-contact kink) + quad (d**2 y / dx**2, allow decreasing sensitivity with increased x) + """ + high_voltage_rail = 2**16 - 1 # bits + x_contact,y_contact,slope,quad = params + y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact + y = numpy.clip(y, y_contact, high_voltage_rail) + return y + +def limited_quadratic_param_guess(x, y) : + """ + Guess rough parameters for a limited_quadratic model. Assumes the + bump approaches (raising the deflection as it does so) first. + Retracting after the approach is optional. Approximates the contact + position and an on-surface (high) position by finding first crossings + of thresholds 0.3 and 0.7 of the y value's total range. Not the + most efficient algorithm, but it seems fairly robust. + """ + x_contact,y_contact,slope = limited_linear_param_guess(x,y) quad = 0 return (x_contact, y_contact, slope, quad) @@ -162,7 +199,7 @@ def limited_quadratic_sensitivity(params): @splittableKwargsFunction() def bump_fit(zpiezo_output_bits, deflection_input_bits, - paramGuesser=limited_quadratic_param_guess, + param_guesser=limited_quadratic_param_guess, model=limited_quadratic, sensitivity_from_fit_params=limited_quadratic_sensitivity, plotVerbose=False) : @@ -170,9 +207,9 @@ def bump_fit(zpiezo_output_bits, deflection_input_bits, y = deflection_input_bits def residual(p, y, x) : return model(x, p) - y - paramGuess = paramGuesser(x, y) + param_guess = param_guesser(x, y) p,cov,info,mesg,ier = \ - scipy.optimize.leastsq(residual, paramGuess, args=(y, x), + scipy.optimize.leastsq(residual, param_guess, args=(y, x), full_output=True, maxfev=int(10e3)) if config.TEXT_VERBOSE : print "Fitted params:",p @@ -184,7 +221,7 @@ def bump_fit(zpiezo_output_bits, deflection_input_bits, else : print "Solution did not converge" if plotVerbose or config.PYLAB_VERBOSE : - yguess = model(x, paramGuess) + yguess = model(x, param_guess) #yguess = None # Don't print the guess, since I'm convinced it's ok ;). yfit = model(x, p) bump_plot(data={"Z piezo output":x, "Deflection input":y}, @@ -310,8 +347,9 @@ if __name__ == '__main__' : 'and one to analyze tweak files.\n' '\n' 'Single file mode (the default) :\n' - 'Scales raw DAC/ADC bit data and fits a straight line.\n' - 'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm.\n' + 'Scales raw DAC/ADC bit data and fits a bounded quadratic.\n' + 'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm, determined by.\n' + 'the slope at the kink between the non-contact region and the contact region.\n' ' should be whitespace-delimited, 2 column ASCII\n' 'without a header line. e.g: "\\t\\n"\n' '\n' @@ -341,6 +379,9 @@ if __name__ == '__main__' : parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true', help='Run input files with datalogger.read_dict_of_arrays(). This is useful, for example, to test a single line from a tweakfile.', default=False) + parser.add_option('-q', '--disable-quadratic', dest='quadratic', action='store_false', + help='Disable quadratic term in fitting (i.e. use bounded linear fits).', + default=True) parser.add_option('-v', '--verbose', dest='verbose', action='store_true', help='Print lots of debugging information', default=False) @@ -359,17 +400,31 @@ if __name__ == '__main__' : config.PYLAB_INTERACTIVE = False config.PYLAB_VERBOSE = options.pylab config.GNUPLOT_VERBOSE = False + if options.quadratic == True: + param_guesser = limited_quadratic_param_guess + model = limited_quadratic + sensitivity_from_fit_params = limited_quadratic_sensitivity + else: + param_guesser = limited_linear_param_guess + model = limited_linear + sensitivity_from_fit_params = limited_linear_sensitivity if options.tweakmode == False : if options.datalogger_mode: data = bump_load(ifilename) else: data = read_data(ifilename) - photoSensitivity = bump_analyze(data) + photoSensitivity = bump_analyze(data, + param_guesser=param_guesser, + model=model, + sensitivity_from_fit_params=sensitivity_from_fit_params) print >> ofile, photoSensitivity else : # tweak file mode - slopes = bump_load_analyze_tweaked(ifilename) + slopes = bump_load_analyze_tweaked(ifilename, + param_guesser=param_guesser, + model=model, + sensitivity_from_fit_params=sensitivity_from_fit_params) if options.comma_out : sep = ',' else : diff --git a/calibcant/vib_analyze.py b/calibcant/vib_analyze.py index a3eddf0..77b8402 100755 --- a/calibcant/vib_analyze.py +++ b/calibcant/vib_analyze.py @@ -95,9 +95,11 @@ def vib_analyze(deflection_bits, freq, print "Converting %d bit values to voltages" % len(Vphoto_t) for i in range(len(deflection_bits)) : Vphoto_t[i] = Vphoto_in2V(deflection_bits[i]) - + if config.TEXT_VERBOSE : + print "Average deflection (Volts):", Vphoto_t.mean() + # Compute the average power spectral density per unit time (in uV**2/Hz) - if config.TEXT_VERBOSE : + if config.TEXT_VERBOSE : print "Compute the averaged power spectral density in uV**2/Hz" (freq_axis, power) = \ FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6, @@ -105,7 +107,7 @@ def vib_analyze(deflection_bits, freq, A,B,C,D = fit_psd(freq_axis, power, **kwargs) - if config.TEXT_VERBOSE : + if config.TEXT_VERBOSE : print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with" print "A = %g, \t B = %g, \t C = %g, \t D = %g" % (A, B, C, D)