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 21e13cb308c6281f626fc2009d00baf9dbe036d1..37839367d5b70420aba06595f9dc5518ace812cd 100755 (executable)
--- 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'
                     '<input-file> should be whitespace-delimited, 2 column ASCII\n'
                     'without a header line.  e.g: "<zp_DAC>\\t<deflection_ADC>\\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 a3eddf0d7d96900787a7cc921f1291f0015b7fa1..77b84020d24ee44d7172bb822a8ed675052d5480 100755 (executable)
--- 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)