"""
import os, time, numpy
-import GnuplotBiDir # used for fitting lorentzian
+import GnuplotBiDir # used for fitting lorentzian
+import scipy.optimize # alternative for fitting lorentzian
import common # common module for the calibcant package
import config # config module for the calibcant package
import data_logger
import FFT_tools
+from splittable_kwargs import splittableKwargsFunction, \
+ make_splittable_kwargs_function
+PLOT_GUESSED_LORENTZIAN=False
+
+@splittableKwargsFunction()
def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
zeroVphoto_bits=config.zeroVphoto_bits) :
"""
Vphoto_in2V is a function converting Vphoto values from bits to Volts.
This function is assumed linear, see bump_analyze().
- zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
"""
Vphoto_std_bits = deflection_bits.std()
- Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
+ Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zeroVphoto_bits)
return Vphoto_std**2
-
-def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000,
+
+#@splittableKwargsFunction((fit_psd, 'freq_axis', 'psd_data'),
+# (vib_plot, 'deflection_bits', 'freq_axis', 'power',
+# 'A', 'B', 'C', 'minFreq', 'maxFreq'))
+# Some of the child functions aren't yet defined, so postpone
+# make-splittable until later in the module.
+def vib_analyze(deflection_bits, freq,
chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V, plotVerbose=False) :
+ Vphoto_in2V=config.Vphoto_in2V, **kwargs) :
"""
Calculate the variance in the raw data due to the cantilever's
thermal oscillation and convert it to Volts**2.
- Improves on vib_analyze_naive() by first converting Vphoto(t) to
- frequency space, and fitting a Lorentzian in the relevant frequency
- range (see cantilever_calib.pdf for derivation).
+ Improves on vib_analyze_naive() by first converting Vphoto(t) to
+ frequency space, and fitting a Lorentzian in the relevant
+ frequency range (see cantilever_calib.pdf for derivation).
+ However, there may be cases where the fit is thrown off by noise
+ spikes in frequency space. To protect from errors, the fitted
+ variance is compared to the naive variance (where all noise is
+ included), and the minimum variance is returned.
- The conversion to frequency space generates an average power spectrum
- by breaking the data into windowed chunks and averaging the power
- spectrums for the chunks together.
- See FFT_tools.avg_power_spectrum
-
Vphoto_in2V is a function converting Vphoto values from bits to Volts.
This function is assumed linear, see bump_analyze().
- zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
"""
+ fit_psd_kwargs,vib_plot_kwargs = vib_analyze._splitargs(vib_analyze,kwargs)
+ if 'minFreq' in fit_psd_kwargs:
+ vib_plot_kwargs['minFreq'] = fit_psd_kwargs['minFreq']
+ if 'maxFreq' in fit_psd_kwargs:
+ vib_plot_kwargs['maxFreq'] = fit_psd_kwargs['maxFreq']
+
Vphoto_t = numpy.zeros((len(deflection_bits),),
dtype=numpy.float)
# convert the data from bits to volts
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, minFreq, maxFreq)
+ A,B,C = 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)
- vib_plot(deflection_bits, freq_axis, power, A, B, C, plotVerbose=plotVerbose)
+ vib_plot(deflection_bits, freq_axis, power, A, B, C, **kwargs)
# Our A is in uV**2, so convert back to Volts**2
- return lorentzian_area(A,B,C) * 1e-12
+ fitted_variance = lorentzian_area(A,B,C) * 1e-12
+ 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
+
+ return min(fitted_variance, naive_variance)
+
+@splittableKwargsFunction()
def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
"""
freq_axis : array of frequencies in Hz
Calculate the variance in the raw data due to the cantilever's
thermal oscillation and convert it to Volts**2.
- Improves on vib_analyze_naive() by working on frequency space data
- and fitting a Lorentzian in the relevant frequency range (see
- cantilever_calib.pdf for derivation).
-
The conversion to frequency space generates an average power spectrum
by breaking the data into windowed chunks and averaging the power
spectrums for the chunks together.
- See FFT_tools.unitary_avg_power_spectrum().
+ See FFT_tools.unitary_avg_power_spectrum().
"""
# cut out the relevent frequency range
if config.TEXT_VERBOSE :
assert imax >= imin + 10 , \
"Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
- # save about-to-be-fitted data to a temporary file
- datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time())
- fd = file(datafilename, 'w')
- for i in range(imin, imax) :
- fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i]))
- fd.close()
-
# generate guesses for Lorentzian parameters A,B,C
max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin
max_psd = psd_data[max_psd_index]
# peak height = C / (3 x_res^4 + A^4)
# Q = A/B
#
- # guess Q = 5, and adjust from there
- Q_guess = 5
+ # guess Q = 1, and adjust from there
+ Q_guess = 1
# so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
# B = x_res / sqrt(Q^2-1/2)
if config.TEXT_VERBOSE :
# (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
# so w is a disaster ;)
# For some values of A and B (non-underdamped), W is imaginary or negative.
+ #
+ # The mass m is given by m = k_B T / (pi A)
+ # The spring constant k is given by k = m (omega_0)**2
+ # The quality factor Q is given by Q = omega_0 m / gamma
if config.TEXT_VERBOSE :
print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
+ if PLOT_GUESSED_LORENTZIAN:
+ vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess,
+ minFreq, maxFreq, plotVerbose=True)
- # fit Lorentzian using Gnuplot's 'fit' command
- g = GnuplotBiDir.Gnuplot()
- # The Lorentzian is the predicted one-sided power spectral density
- # For an oscillator whose position x obeys
- # m x'' + gamma x' + kx = F_thermal(t)
- # A is the ?driving noise?
- # B is gamma, the damping term
- # C is omega_0, the resonant frequency
# Fitting the PSD of V = photoSensitivity*x just rescales the parameters
- g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)")
- ## The mass m is given by m = k_B T / (pi A)
- ## The spring constant k is given by k = m (omega_0)**2
- ## The quality factor Q is given by Q = omega_0 m / gamma
- g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
- g("fit f(x) '%s' via A,B,C" % datafilename)
- A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x)
- B=abs(float(g.getvar('B'))) # so ensure we get positive values
- C=float(g.getvar('C'))
- os.remove(datafilename)
+
+ if False: # Gnuplot fit worse than scipy.optimize.leastsq fit.
+ # fit Lorentzian using Gnuplot's 'fit' command
+
+ # save about-to-be-fitted data to a temporary file
+ datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time())
+ fd = file(datafilename, 'w')
+ for i in range(imin, imax) :
+ fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i]))
+ fd.close()
+
+ g = GnuplotBiDir.Gnuplot()
+ g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)")
+ g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
+ g("fit f(x) '%s' via A,B,C" % datafilename)
+ A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x)
+ B=abs(float(g.getvar('B'))) # so ensure we get positive values
+ C=float(g.getvar('C'))
+
+ os.remove(datafilename)
+ else:
+ # fit Lorenzian using scipy.optimize.leastsq
+ def residual(p, y, x):
+ A = p[0]
+ B = p[1]
+ C = p[2]
+ Y = C / ((A**2-x**2)**2 + (B*x)**2)
+ return Y-y
+ leastsq = scipy.optimize.leastsq
+ p,cov,info,mesg,ier = leastsq(residual, (A_guess, B_guess, C_guess),
+ args=(psd_data[imin:imax],
+ freq_axis[imin:imax]),
+ full_output=True, maxfev=10000)
+ if config.TEXT_VERBOSE:
+ print "Fitted params:",p
+ print "Covariance mx:",cov
+ #print "Info:", info # too much information :p
+ print "mesg:", mesg
+ if ier == 1 :
+ print "Solution converged"
+ else :
+ print "Solution did not converge"
+ A,B,C = p
+ 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)
def lorentzian_area(A, B, C) :
string += "plot '%s', f(x)\n" % datafilename
return string
+@splittableKwargsFunction()
def vib_save(data, log_dir=None) :
"""Save the dictionary data, using data_logger.data_log()
data should be a dictionary of arrays with the fields
'Deflection input'
"""
if log_dir :
+ freq = data['sample frequency Hz']
log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
log_name="vib_%gHz" % freq)
log.write_dict_of_arrays(data)
data = dl.read_dict_of_arrays(datafile)
return data
+@splittableKwargsFunction()
def vib_plot(deflection_bits, freq_axis, power, A, B, C,
- plotVerbose=True) :
+ minFreq=None, maxFreq=None, plotVerbose=False) :
"""
If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
Time series (Vphoto vs sample index) (show any major drift events),
"""
if plotVerbose or config.PYLAB_VERBOSE :
common._import_pylab()
+ common._pylab.figure(config.BASE_FIGNUM+2)
+
+ if deflection_bits != None:
+ # plot time series
+ common._pylab.subplot(311)
+ common._pylab.hold(False)
+ common._pylab.plot(deflection_bits, 'r.')
+ common._pylab.title("free oscillation")
+
+ # plot histogram distribution and gaussian fit
+ common._pylab.subplot(312)
+ common._pylab.hold(False)
+ n, bins, patches = \
+ common._pylab.hist(deflection_bits, bins=30,
+ normed=1, align='center')
+ gaus = numpy.zeros((len(bins),), dtype=numpy.float)
+ mean = deflection_bits.mean()
+ std = deflection_bits.std()
+ pi = numpy.pi
+ exp = numpy.exp
+ for i in range(len(bins)) :
+ gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
+ common._pylab.hold(True)
+ common._pylab.plot(bins, gaus, 'r-');
+ common._pylab.hold(False)
+
+ # plot FFTed data
+ axes = common._pylab.subplot(313)
+ else:
+ # use a nice big subplot just for the FFTed data
+ axes = common._pylab.subplot(111)
common._pylab.hold(False)
- common._pylab.figure(BASE_FIGNUM+2)
-
- # plot time series
- common._pylab.subplot(311)
- common._pylab.plot(data["Deflection input"], 'r.')
- common._pylab.title("free oscillation")
-
- # plot histogram distribution and gaussian fit
- common._pylab.subplot(312)
- n, bins, patches = \
- common._pylab.hist(data["Deflection input"], bins=30,
- normed=1, align='center')
- gaus = numpy.zeros((len(bins),))
- mean = data["Deflection input"].mean()
- std = data["Deflection input"].std()
- pi = numpy.pi
- exp = numpy.exp
- for i in range(len(bins)) :
- gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
+ common._pylab.semilogy(freq_axis, power, 'r.-')
common._pylab.hold(True)
- common._pylab.plot(bins, gaus, 'r-');
+ xmin,xmax = axes.get_xbound()
+ ymin,ymax = axes.get_ybound()
+
+ if minFreq is not None and maxFreq is not None:
+ # highlight the region we're fitting in
+ 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-');
common._pylab.hold(False)
+ axes.axis([xmin,xmax,ymin,ymax])
- # plot FFTed data
- common._pylab.subplot(313)
- common._pylab.semilogy(freq_axis, power, 'r.-')
common._flush_plot()
- if (plotVerbose or config.GNUPLOT_VERBOSE) and False : # TODO, clean up and double check...
+ if (plotVerbose or config.GNUPLOT_VERBOSE) and False: # TODO: cleanup and test
# write all the ft data now
fd = file(datafilename, 'w')
for i in range(len(freq_axis)) :
g.getvar("MOUSE_BUTTON")
del(g)
-def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
- chunk_size=2048, overlap=True,
- window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V,
- plotVerbose=False) :
+# Make vib_analyze splittable (was postponed until the child functions
+# were defined).
+make_splittable_kwargs_function(vib_analyze,
+ (fit_psd, 'freq_axis', 'psd_data'),
+ (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C',
+ 'minFreq', 'maxFreq'))
+
+@splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq'))
+def vib_load_analyze_tweaked(tweak_file, **kwargs) :
"""
Read the vib files listed in tweak_file.
Return an array of Vphoto variances (due to the cantilever) in Volts**2
"""
+ vib_analyze_kwargs, = \
+ vib_load_analyze_tweaked._splitargs(vib_load_analyze_tweaked, kwargs)
dl = data_logger.data_load()
Vphoto_var = []
for path in file(tweak_file, 'r') :
if config.TEXT_VERBOSE :
print "Analyzing '%s' at %g Hz" % (path, freq)
# analyze
- Vphoto_var.append(vib_analyze(deflection_bits, freq, minFreq, maxFreq,
- chunk_size, overlap, window,
- Vphoto_in2V, plotVerbose))
+ Vphoto_var.append(vib_analyze(deflection_bits, freq,
+ **vib_analyze_kwargs))
return numpy.array(Vphoto_var, dtype=numpy.float)
# commandline interface functions
'2008, W. Trevor King.\n'
'\n'
'Deflection power spectral density (PSD) data (X^2/Hz)\n'
- 'and returns the variance in X units (X^2)'
- '<input-file> should be whitespace-delimited, 2 column ASCII\n'
- 'without a header line.\n'
- 'e.g: "<frequency_Hz>\\t<deflection_psd X^2/Hz>\\n"\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')
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.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
help='Print gnuplot fit check script to stderr',
default=False)
+ parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
+ help='Produce pylab fit checks during execution',
+ default=False)
+ parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
+ help='Produce plots of the pre-fit Lorentzian',
+ default=False)
parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
help='Run in tweak-file mode',
default=False)
+ parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
+ help='Read 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('-v', '--verbose', dest='verbose', action='store_true',
help='Print lots of debugging information',
default=False)
ofile = file(options.ofilename, 'w')
else :
ofile = sys.stdout
- if options.verbose == True :
- vfile = sys.stderr
- else :
- vfile = None
config.TEXT_VERBOSE = options.verbose
- config.PYLAB_VERBOSE = False
+ config.PYLAB_INTERACTIVE = False
+ config.PYLAB_VERBOSE = options.pylab
config.GNUPLOT_VERBOSE = options.gnuplot
+ PLOT_GUESSED_LORENTZIAN = options.plot_guess
if options.tweakmode == False :
- data = read_data(ifilename)
- deflection_bits = data['Deflection input']
+ if options.datalogger_mode:
+ data = vib_load(ifilename)
+ deflection_bits = data['Deflection input']
+ else:
+ deflection_bits = read_data(ifilename)
Vphoto_var = vib_analyze(deflection_bits, freq=options.freq,
minFreq=options.min_freq,
maxFreq=options.max_freq)
sep = '\n'
common.write_array(ofile, Vphoto_vars, sep)
+ if common._final_flush_plot != None:
+ common._final_flush_plot()
+
if options.ofilename != None :
ofile.close()
+