From f8efb4c0ce7ee4e1ef9bc0b97fbd4ce14fd0ee25 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Tue, 27 Jan 2009 09:30:41 -0500 Subject: [PATCH] Various adjustments. I should commit more often ;). Added # comments to tweakfile syntax. Played around with adding a white-noise floor in the vibration fitting, but the fits didn't look all that convincing. Some of the white-noise code is still in place, but I think it's currently disabled ;). Fixed some typos in TEXT_VERBOSE output in vib_analyze.py --- calibcant/T_analyze.py | 5 ++++- calibcant/analyze.py | 13 +++++++++---- calibcant/bump_analyze.py | 5 ++++- calibcant/common.py | 2 +- calibcant/config.py | 4 ++-- calibcant/vib_analyze.py | 33 ++++++++++++++++++++------------- 6 files changed, 40 insertions(+), 22 deletions(-) diff --git a/calibcant/T_analyze.py b/calibcant/T_analyze.py index 6f0a8a1..d4ddf72 100755 --- a/calibcant/T_analyze.py +++ b/calibcant/T_analyze.py @@ -106,7 +106,9 @@ def T_load_analyze_tweaked(tweak_file, convert_to_K=C_to_K, textVerboseFile=None Ts = [] for line in file(tweak_file, 'r') : parsed = line.split() - path = parsed[0].split('\n')[0] + path = parsed[0].strip() + if path[0] == '#': # a comment + continue if textVerboseFile != None : print >> textVerboseFile, "Reading data from %s" % (path) # read the data @@ -142,6 +144,7 @@ if __name__ == '__main__' : 'Tweak file mode:\n' 'Runs the same analysis as in single file mode for each T file in\n' 'a tweak file. Each line in the tweak file specifies a single T file.\n' + 'Blank lines and those beginning with a pound sign (#) are ignored.\n' 'A T file contains a sequence of 32 bit floats representing temperature in K.\n' ) parser = OptionParser(usage=usage_string, version='%prog 0.1') diff --git a/calibcant/analyze.py b/calibcant/analyze.py index 7ed7564..c485009 100755 --- a/calibcant/analyze.py +++ b/calibcant/analyze.py @@ -269,19 +269,24 @@ if __name__ == '__main__' : options,args = parser.parse_args() parser.destroy() - + + config.TEXT_VERBOSE = options.verbose + config.PYLAB_INTERACTIVE = False + config.PYLAB_VERBOSE = options.plot + bumps = get_array(options.bump_string, options.bump_file, "bump") temps = get_array(options.temp_string, options.temp_file, "temp") vibs = get_array(options.vib_string, options.vib_file, "vib") - if options.plot : - calib_plot(bumps, temps, vibs) + #if options.plot : + # calib_plot(bumps, temps, vibs) if options.celsius : for i in range(len(temps)) : temps[i] = T_analyze.C_to_K(temps[i]) - km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = calib_analyze(bumps, temps, vibs) + km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = \ + calib_analyze(bumps, temps, vibs, plotVerbose=options.plot) if options.verbose : print >> sys.stderr, \ diff --git a/calibcant/bump_analyze.py b/calibcant/bump_analyze.py index 11a8f51..36a6b81 100755 --- a/calibcant/bump_analyze.py +++ b/calibcant/bump_analyze.py @@ -128,7 +128,9 @@ def bump_load_analyze_tweaked(tweak_file, **kwargs): photoSensitivity = [] for line in file(tweak_file, 'r') : parsed = line.split() - path = parsed[0].split('\n')[0] + path = parsed[0].strip() + if path[0] == '#' : # a comment + continue if config.TEXT_VERBOSE : print "Reading data from %s with ranges %s" % (path, parsed[1:]) # read the data @@ -192,6 +194,7 @@ if __name__ == '__main__' : 'Tweak file mode:\n' 'Runs the same analysis as in single file mode for each bump in\n' 'a tweak file. Each line in the tweak file specifies a single bump.\n' + 'Blank lines and those beginning with a pound sign (#) are ignored.\n' 'The format of a line is a series of whitespace-separated fields--\n' 'a base file path followed by optional point index ranges, e.g.:\n' '20080919/20080919132500_bump_surface 10:651 1413:2047\n' diff --git a/calibcant/common.py b/calibcant/common.py index c58e06b..1e8ff71 100644 --- a/calibcant/common.py +++ b/calibcant/common.py @@ -31,7 +31,7 @@ _final_flush_plot = None _pylab = None def _dummy_fn(): pass -def _import_pylab() : +def _import_pylab() : # TODO: auto detect no DISPLAY and abort """Import pylab plotting functions for when we need to plot. This function can be called multiple times, and ensures that the pylab setup is only imported once. It defines the functions diff --git a/calibcant/config.py b/calibcant/config.py index 8d7d511..4dbd1a2 100644 --- a/calibcant/config.py +++ b/calibcant/config.py @@ -31,8 +31,8 @@ LOG_DATA = True # quietly grab all real-world data and log to LOG_DIR LOG_DIR = '$DEFAULT$/calibrate_cantilever' GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata' TEXT_VERBOSE = True # for debugging -GNUPLOT_VERBOSE = True # turn on fit check plotting -PYLAB_VERBOSE = True # turn on plotting +GNUPLOT_VERBOSE = True # turn on fit check plotting +PYLAB_VERBOSE = False # turn on plotting PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots BASE_FIGNUM = 20 # to avoid writing to already existing figures diff --git a/calibcant/vib_analyze.py b/calibcant/vib_analyze.py index a98ea05..ff37f12 100755 --- a/calibcant/vib_analyze.py +++ b/calibcant/vib_analyze.py @@ -103,13 +103,13 @@ def vib_analyze(deflection_bits, freq, FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6, freq, chunk_size,overlap, window) - A,B,C = fit_psd(freq_axis, power, **kwargs) + A,B,C,D = fit_psd(freq_axis, power, **kwargs) 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" % (A, B, C) + print "A = %g, \t B = %g, \t C = %g, \t D = %g" % (A, B, C, D) - vib_plot(deflection_bits, freq_axis, power, A, B, C, **kwargs) + vib_plot(deflection_bits, freq_axis, power, A, B, C, D, **kwargs) # Our A is in uV**2, so convert back to Volts**2 fitted_variance = lorentzian_area(A,B,C) * 1e-12 @@ -117,8 +117,8 @@ def vib_analyze(deflection_bits, freq, naive_variance = vib_analyze_naive(deflection_bits, Vphoto_in2V=Vphoto_in2V) if config.TEXT_VERBOSE: - print "Fitted variance: %g V**2", fitted_variance - print "Naive variance : %g V**2", naive_variance + print "Fitted variance: %g V**2" % fitted_variance + print "Naive variance : %g V**2" % naive_variance return min(fitted_variance, naive_variance) @@ -175,6 +175,7 @@ def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) : B_guess = res_freq / numpy.sqrt(Q_guess**2-0.5) A_guess = Q_guess*B_guess C_guess = max_psd * (-res_freq**4 + A_guess**4) + D_guess = 0 # allow a noise floor offset (DISABLED in fitting below) # # Half width w on lower side when L(a-w) = L(a)/2 # (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2 @@ -241,9 +242,10 @@ def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) : else : print "Solution did not converge" A,B,C = p + D = 0 # do not fit the noise floor (fit doesn't look very convincing) A=abs(A) # A and B only show up as squares in f(x) B=abs(B) # so ensure we get positive values - return (A, B, C) + return (A, B, C, D) def lorentzian_area(A, B, C) : # Integrating the the power spectral density per unit time (power) over the @@ -290,7 +292,7 @@ def vib_load(datafile=None) : return data @splittableKwargsFunction() -def vib_plot(deflection_bits, freq_axis, power, A, B, C, +def vib_plot(deflection_bits, freq_axis, power, A, B, C, D, minFreq=None, maxFreq=None, plotVerbose=False) : """ If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures: @@ -342,8 +344,10 @@ def vib_plot(deflection_bits, freq_axis, power, A, B, C, common._pylab.axvspan(minFreq, maxFreq, facecolor='g', alpha=0.1, zorder = -1) - fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2) - common._pylab.plot(freq_axis, fitdata, 'b-'); + fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2) + D + common._pylab.plot(freq_axis, fitdata, 'b-') + noisefloor = D + 0*freq_axis; + common._pylab.plot(freq_axis, noisefloor) common._pylab.hold(False) axes.axis([xmin,xmax,ymin,ymax]) @@ -381,7 +385,7 @@ def vib_plot(deflection_bits, freq_axis, power, A, B, C, # were defined). make_splittable_kwargs_function(vib_analyze, (fit_psd, 'freq_axis', 'psd_data'), - (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C', + (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C', 'D', 'minFreq', 'maxFreq')) @splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq')) @@ -395,8 +399,9 @@ def vib_load_analyze_tweaked(tweak_file, **kwargs) : dl = data_logger.data_load() Vphoto_var = [] for path in file(tweak_file, 'r') : - if path[-1] == '\n': - path = path.split('\n')[0] + path = path.strip() + if path[0] == '#': # a comment + continue # read the data data = vib_load(path) deflection_bits = data['Deflection input'] @@ -432,7 +437,9 @@ if __name__ == '__main__' : 'Deflection power spectral density (PSD) data (X^2/Hz)\n' 'and returns the variance in X units (X^2)\n' ' should be 1 column ASCII without a header line.\n' - 'e.g: "\\n..."\n') + 'e.g: "\\n..."\n' + #TODO, better overview of input modes, e.g. a la T_analyze + ) parser = OptionParser(usage=usage_string, version='%prog 0.1') parser.add_option('-f', '--sample-freq', dest='freq', default=100e3, help='sample frequency in Hz (default %default)', -- 2.26.2