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
Ts = []
for line in file(tweak_file, 'r') :
parsed = line.split()
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
if textVerboseFile != None :
print >> textVerboseFile, "Reading data from %s" % (path)
# read the data
'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'
'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')
'A T file contains a sequence of 32 bit floats representing temperature in K.\n'
)
parser = OptionParser(usage=usage_string, version='%prog 0.1')
options,args = parser.parse_args()
parser.destroy()
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")
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])
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, \
if options.verbose :
print >> sys.stderr, \
photoSensitivity = []
for line in file(tweak_file, 'r') :
parsed = line.split()
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
if config.TEXT_VERBOSE :
print "Reading data from %s with ranges %s" % (path, parsed[1:])
# read the data
'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'
'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'
'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'
_pylab = None
def _dummy_fn(): pass
_pylab = None
def _dummy_fn(): pass
+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
"""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
LOG_DIR = '$DEFAULT$/calibrate_cantilever'
GNUFIT_DATA_BASE='./calibrate_cantilever_fitdata'
TEXT_VERBOSE = True # for debugging
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
PYLAB_INTERACTIVE = True # select between draw() and show() for flushing plots
BASE_FIGNUM = 20 # to avoid writing to already existing figures
FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
freq, chunk_size,overlap, window)
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"
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
# Our A is in uV**2, so convert back to Volts**2
fitted_variance = lorentzian_area(A,B,C) * 1e-12
naive_variance = vib_analyze_naive(deflection_bits,
Vphoto_in2V=Vphoto_in2V)
if config.TEXT_VERBOSE:
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)
return min(fitted_variance, naive_variance)
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)
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
#
# 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
else :
print "Solution did not converge"
A,B,C = p
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
A=abs(A) # A and B only show up as squares in f(x)
B=abs(B) # so ensure we get positive values
def lorentzian_area(A, B, C) :
# Integrating the the power spectral density per unit time (power) over the
def lorentzian_area(A, B, C) :
# Integrating the the power spectral density per unit time (power) over the
return data
@splittableKwargsFunction()
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:
minFreq=None, maxFreq=None, plotVerbose=False) :
"""
If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
common._pylab.axvspan(minFreq, maxFreq, facecolor='g', alpha=0.1,
zorder = -1)
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])
common._pylab.hold(False)
axes.axis([xmin,xmax,ymin,ymax])
# were defined).
make_splittable_kwargs_function(vib_analyze,
(fit_psd, 'freq_axis', 'psd_data'),
# 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'))
'minFreq', 'maxFreq'))
@splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq'))
dl = data_logger.data_load()
Vphoto_var = []
for path in file(tweak_file, 'r') :
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']
# read the data
data = vib_load(path)
deflection_bits = data['Deflection input']
'Deflection power spectral density (PSD) data (X^2/Hz)\n'
'and returns the variance in X units (X^2)\n'
'<input-file> should be 1 column ASCII without a header line.\n'
'Deflection power spectral density (PSD) data (X^2/Hz)\n'
'and returns the variance in X units (X^2)\n'
'<input-file> should be 1 column ASCII without a header line.\n'
- 'e.g: "<deflection bits>\\n..."\n')
+ 'e.g: "<deflection bits>\\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)',
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)',