calibcant/vib_analyze.py

   1 #!/usr/bin/python
   2 #
   3 # calibcant - tools for thermally calibrating AFM cantilevers
   4 #
   5 # Copyright (C) 2007,2008, William Trevor King
   6 #
   7 # This program is free software; you can redistribute it and/or
   8 # modify it under the terms of the GNU General Public License as
   9 # published by the Free Software Foundation; either version 3 of the
  10 # License, or (at your option) any later version.
  11 #
  12 # This program is distributed in the hope that it will be useful, but
  13 # WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  15 # See the GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with this program; if not, write to the Free Software
  19 # Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  20 # 02111-1307, USA.
  21 #
  22 # The author may be contacted at <wking@drexel.edu> on the Internet, or
  23 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
  24 # Philadelphia PA 19104, USA.
  25
  26 """
  27 Separate the more general vib_analyze() from the other vib_*()
  28 functions in calibcant.  Also provide a command line interface
  29 for analyzing data acquired through other workflows.
  30
  31 The relevent physical quantities are :
  32  Vphoto   The photodiode vertical deflection voltage (what we measure)
  33 """
  34
  35 import os, time, numpy
  36 import GnuplotBiDir   # used for fitting lorentzian
  37 import scipy.optimize # alternative for fitting lorentzian
  38 import common # common module for the calibcant package
  39 import config # config module for the calibcant package
  40 import data_logger
  41 import FFT_tools
  42
  43 PLOT_GUESSED_LORENTZIAN=False
  44
  45 def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
  46                       zeroVphoto_bits=config.zeroVphoto_bits) :
  47     """
  48     Calculate the variance of the raw data, and convert it to Volts**2.
  49     This method is simple and easy to understand, but it highly succeptible
  50     to noise, drift, etc.
  51
  52     Vphoto_in2V is a function converting Vphoto values from bits to Volts.
  53     This function is assumed linear, see bump_analyze().
  54     """
  55     Vphoto_std_bits = deflection_bits.std()
  56     Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zeroVphoto_bits)
  57     return Vphoto_std**2
  58
  59 def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000,
  60                 chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
  61                 Vphoto_in2V=config.Vphoto_in2V, plotVerbose=False) :
  62     """
  63     Calculate the variance in the raw data due to the cantilever's
  64     thermal oscillation and convert it to Volts**2.
  65
  66     Improves on vib_analyze_naive() by first converting Vphoto(t) to
  67     frequency space, and fitting a Lorentzian in the relevant
  68     frequency range (see cantilever_calib.pdf for derivation).
  69     However, there may be cases where the fit is thrown off by noise
  70     spikes in frequency space.  To protect from errors, the fitted
  71     variance is compared to the naive variance (where all noise is
  72     included), and the minimum variance is returned.
  73
  74     Vphoto_in2V is a function converting Vphoto values from bits to Volts.
  75     This function is assumed linear, see bump_analyze().
  76     """
  77     Vphoto_t = numpy.zeros((len(deflection_bits),),
  78                            dtype=numpy.float)
  79     # convert the data from bits to volts
  80     if config.TEXT_VERBOSE :
  81         print "Converting %d bit values to voltages" % len(Vphoto_t)
  82     for i in range(len(deflection_bits)) :
  83         Vphoto_t[i] = Vphoto_in2V(deflection_bits[i])
  84
  85     # Compute the average power spectral density per unit time (in uV**2/Hz)
  86     if config.TEXT_VERBOSE :
  87         print "Compute the averaged power spectral density in uV**2/Hz"
  88     (freq_axis, power) = \
  89         FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
  90                                              freq, chunk_size,overlap, window)
  91
  92     A,B,C = fit_psd(freq_axis, power, minFreq, maxFreq)
  93
  94     if config.TEXT_VERBOSE :
  95         print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
  96         print "A = %g, \t B = %g, \t C = %g" % (A, B, C)
  97
  98     vib_plot(deflection_bits, freq_axis, power, A, B, C,
  99              minFreq, maxFreq, plotVerbose=plotVerbose)
 100
 101     # Our A is in uV**2, so convert back to Volts**2
 102     fitted_variance = lorentzian_area(A,B,C) * 1e-12
 103
 104     naive_variance = vib_analyze_naive(deflection_bits,
 105                                        Vphoto_in2V=Vphoto_in2V)
 106     if config.TEXT_VERBOSE:
 107         print "Fitted variance: %g V**2", fitted_variance
 108         print "Naive variance : %g V**2", naive_variance
 109
 110     return min(fitted_variance, naive_variance)
 111
 112 def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
 113     """
 114     freq_axis : array of frequencies in Hz
 115     psd_data  : array of PSD amplitudes for the frequencies in freq_axis
 116
 117     Calculate the variance in the raw data due to the cantilever's
 118     thermal oscillation and convert it to Volts**2.
 119
 120     The conversion to frequency space generates an average power spectrum
 121     by breaking the data into windowed chunks and averaging the power
 122     spectrums for the chunks together.
 123     See FFT_tools.unitary_avg_power_spectrum().
 124     """
 125     # cut out the relevent frequency range
 126     if config.TEXT_VERBOSE :
 127         print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
 128     imin = 0
 129     while freq_axis[imin] < minFreq : imin += 1
 130     imax = imin
 131     while freq_axis[imax] < maxFreq : imax += 1
 132     assert imax >= imin + 10 , \
 133         "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
 134
 135     # generate guesses for Lorentzian parameters A,B,C
 136     max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin
 137     max_psd = psd_data[max_psd_index]
 138     half_max_index = imin
 139     while psd_data[half_max_index] < max_psd/2 :
 140         half_max_index += 1
 141     res_freq = freq_axis[max_psd_index]
 142     half_width = freq_axis[max_psd_index] - freq_axis[half_max_index]
 143     # Lorentzian L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
 144     # is expected power spectrum for
 145     # x'' + B x' + A^2 x'' = F_external(t)/m
 146     # (A = omega_0)
 147     # C = (2 k_B T B) / (pi m)
 148     #
 149     # A = resonant frequency
 150     # peak at  x_res = sqrt(A^2 - B^2/2)
 151     #  which could be complex if there isn't a peak (overdamped)
 152     # peak height    = C / (3 x_res^4 + A^4)
 153     # Q = A/B
 154     #
 155     # guess Q = 1, and adjust from there
 156     Q_guess = 1
 157     # so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
 158     #    B = x_res / sqrt(Q^2-1/2)
 159     if config.TEXT_VERBOSE :
 160         print "params : %g, %g" % (res_freq, max_psd)
 161     B_guess = res_freq / numpy.sqrt(Q_guess**2-0.5)
 162     A_guess = Q_guess*B_guess
 163     C_guess = max_psd * (-res_freq**4 + A_guess**4)
 164     #
 165     # Half width w on lower side when L(a-w) = L(a)/2
 166     #  (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
 167     # Let W=(a-w)**2, A=a**2, and B=b**2
 168     #  (A - W)**2 + BW = 2*AB
 169     #  W**2 - 2AW + A**2 + BW = 2AB
 170     #  W**2 + (B-2A)W + (A**2-2AB) = 0
 171     #  W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
 172     #    = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 173     #  (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 174     #  so w is a disaster ;)
 175     # For some values of A and B (non-underdamped), W is imaginary or negative.
 176     #
 177     # The mass m is given by m = k_B T / (pi A)
 178     # The spring constant k is given by k = m (omega_0)**2
 179     # The quality factor Q is given by Q = omega_0 m / gamma
 180     if config.TEXT_VERBOSE :
 181         print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
 182     if PLOT_GUESSED_LORENTZIAN:
 183         vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess,
 184                  minFreq, maxFreq, plotVerbose=True)
 185
 186     # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
 187
 188     if False: # Gnuplot fit worse than scipy.optimize.leastsq fit.
 189         # fit Lorentzian using Gnuplot's 'fit' command
 190
 191         # save about-to-be-fitted data to a temporary file
 192         datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time())
 193         fd = file(datafilename, 'w')
 194         for i in range(imin, imax) :
 195             fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i]))
 196         fd.close()
 197
 198         g = GnuplotBiDir.Gnuplot()
 199         g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)")
 200         g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
 201         g("fit f(x) '%s' via A,B,C" % datafilename)
 202         A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x)
 203         B=abs(float(g.getvar('B'))) # so ensure we get positive values
 204         C=float(g.getvar('C'))
 205
 206         os.remove(datafilename)
 207     else:
 208         # fit Lorenzian using scipy.optimize.leastsq
 209         def residual(p, y, x):
 210             A = p[0]
 211             B = p[1]
 212             C = p[2]
 213             Y = C / ((A**2-x**2)**2 + (B*x)**2)
 214             return Y-y
 215         leastsq = scipy.optimize.leastsq
 216         p,cov,info,mesg,ier = leastsq(residual, (A_guess, B_guess, C_guess),
 217                                       args=(psd_data[imin:imax],
 218                                             freq_axis[imin:imax]),
 219                                       full_output=True, maxfev=10000)
 220         if config.TEXT_VERBOSE :
 221             print "Fitted params:",p
 222             print "Covariance mx:",cov
 223             print "Info:", info
 224             print "mesg:", mesg
 225             if ier == 1 :
 226                 print "Solution converged"
 227             else :
 228                 print "Solution did not converge"
 229         A,B,C = p
 230         A=abs(A) # A and B only show up as squares in f(x)
 231         B=abs(B) # so ensure we get positive values
 232     return (A, B, C)
 233
 234 def lorentzian_area(A, B, C) :
 235     # Integrating the the power spectral density per unit time (power) over the
 236     # frequency axis [0, fN] returns the total signal power per unit time
 237     #  int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
 238     #                      = <V(t)**2>, the variance for our AC signal.
 239     # The variance from our fitted Lorentzian is the area under the Lorentzian
 240     #  <V(t)**2> = (pi*C) / (2*B*A**2)
 241     return (numpy.pi * C) / (2 * B * A**2)
 242
 243 def gnuplot_check_fit(datafilename, A, B, C) :
 244     """
 245     return a string containing a gnuplot script displaying the fit.
 246     """
 247     string = 'f(x) = C / ((A**2-x**2)**2 + (x*B)**2)\n'
 248     string += 'A= %g ; B= %g ; C= %g\n' % (A, B, C)
 249     string += 'set logscale y\n'
 250     string += "plot '%s', f(x)\n" % datafilename
 251     return string
 252
 253 def vib_save(data, log_dir=None) :
 254     """Save the dictionary data, using data_logger.data_log()
 255     data should be a dictionary of arrays with the fields
 256       'Z piezo output'
 257       'Z piezo input'
 258       'Deflection input'
 259     """
 260     if log_dir :
 261         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 262                                    log_name="vib_%gHz" % freq)
 263         log.write_dict_of_arrays(data)
 264
 265 def vib_load(datafile=None) :
 266     """Load the dictionary data, using data_logger.date_load()"
 267     data should be a dictionary of arrays with the fields
 268       'Z piezo output'
 269       'Z piezo input'
 270       'Deflection input'
 271     """
 272     dl = data_logger.data_load()
 273     data = dl.read_dict_of_arrays(datafile)
 274     return data
 275
 276 def vib_plot(deflection_bits, freq_axis, power, A, B, C,
 277              minFreq=None, maxFreq=None, plotVerbose=True) :
 278     """
 279     If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
 280      Time series (Vphoto vs sample index) (show any major drift events),
 281      A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
 282      FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
 283     """
 284     if plotVerbose or config.PYLAB_VERBOSE :
 285         common._import_pylab()
 286         common._pylab.figure(config.BASE_FIGNUM+2)
 287
 288         if deflection_bits != None:
 289             # plot time series
 290             common._pylab.subplot(311)
 291             common._pylab.hold(False)
 292             common._pylab.plot(deflection_bits, 'r.')
 293             common._pylab.title("free oscillation")
 294
 295             # plot histogram distribution and gaussian fit
 296             common._pylab.subplot(312)
 297             common._pylab.hold(False)
 298             n, bins, patches = \
 299                 common._pylab.hist(deflection_bits, bins=30,
 300                                    normed=1, align='center')
 301             gaus = numpy.zeros((len(bins),), dtype=numpy.float)
 302             mean = deflection_bits.mean()
 303             std = deflection_bits.std()
 304             pi = numpy.pi
 305             exp = numpy.exp
 306             for i in range(len(bins)) :
 307                 gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
 308             common._pylab.hold(True)
 309             common._pylab.plot(bins, gaus, 'r-');
 310             common._pylab.hold(False)
 311
 312             # plot FFTed data
 313             axes = common._pylab.subplot(313)
 314         else:
 315             # use a nice big subplot just for the FFTed data
 316             axes = common._pylab.subplot(111)
 317         common._pylab.hold(False)
 318         common._pylab.semilogy(freq_axis, power, 'r.-')
 319         common._pylab.hold(True)
 320         xmin,xmax = axes.get_xbound()
 321         ymin,ymax = axes.get_ybound()
 322
 323         if minFreq is not None and maxFreq is not None:
 324             # highlight the region we're fitting in
 325             common._pylab.axvspan(minFreq, maxFreq, facecolor='g', alpha=0.1,
 326                                   zorder = -1)
 327
 328         fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2)
 329         common._pylab.plot(freq_axis, fitdata, 'b-');
 330         common._pylab.hold(False)
 331         axes.axis([xmin,xmax,ymin,ymax])
 332
 333         common._flush_plot()
 334     if (plotVerbose or config.GNUPLOT_VERBOSE) and False: # TODO: cleanup and test
 335         # write all the ft data now
 336         fd = file(datafilename, 'w')
 337         for i in range(len(freq_axis)) :
 338             fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
 339         fd.write("\n") # blank line to separate fit data.
 340         # write the fit data
 341         for i in range(imin, imax) :
 342             x = freq_axis[i]
 343             fd.write("%g\t%g\n" % (freq_axis[i],
 344                                    C / ((A**2-x**2)**2 + (x*B)**2) ) )
 345         fd.close()
 346         print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
 347         g("set terminal epslatex color solid")
 348         g("set output 'calibration.tex'")
 349         g("set size 2,2") # multipliers on default 5"x3.5"
 350         #g("set title 'Thermal calibration'")
 351         g("set logscale y")
 352         g("set xlabel 'Frequency (Hz)'")
 353         g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
 354         g("set xrange [0:20000]")
 355         g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
 356         g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
 357
 358         g("set mouse")
 359         g("pause mouse 'click with mouse'")
 360         g.getvar("MOUSE_BUTTON")
 361         del(g)
 362
 363 def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
 364                              chunk_size=2048, overlap=True,
 365                              window=FFT_tools.window_hann,
 366                              Vphoto_in2V=config.Vphoto_in2V,
 367                              plotVerbose=False) :
 368     """
 369     Read the vib files listed in tweak_file.
 370     Return an array of Vphoto variances (due to the cantilever) in Volts**2
 371     """
 372     dl = data_logger.data_load()
 373     Vphoto_var = []
 374     for path in file(tweak_file, 'r') :
 375         if path[-1] == '\n':
 376             path = path.split('\n')[0]
 377         # read the data
 378         data = vib_load(path)
 379         deflection_bits = data['Deflection input']
 380         freq = float(path.split('_')[-1].split('Hz')[0])
 381         if config.TEXT_VERBOSE :
 382             print "Analyzing '%s' at %g Hz" % (path, freq)
 383         # analyze
 384         Vphoto_var.append(vib_analyze(deflection_bits, freq, minFreq, maxFreq,
 385                                       chunk_size, overlap, window,
 386                                       Vphoto_in2V, plotVerbose))
 387     return numpy.array(Vphoto_var, dtype=numpy.float)
 388
 389 # commandline interface functions
 390 import scipy.io, sys
 391
 392 def read_data(ifile):
 393     """
 394     ifile can be a filename string or open (seekable) file object.
 395     returns an array of data values (1 column)
 396     """
 397     #returns (column 1 array, column 2 array)
 398     #"""
 399     if ifile == None :  ifile = sys.stdin
 400     unlabeled_data=scipy.io.read_array(ifile)
 401     return unlabeled_data #(unlabeled_data[:,0], unlabeled_data[:,1])
 402
 403 if __name__ == '__main__' :
 404     # command line interface
 405     from optparse import OptionParser
 406
 407     usage_string = ('%prog <input-file>\n'
 408                     '2008, W. Trevor King.\n'
 409                     '\n'
 410                     'Deflection power spectral density (PSD) data (X^2/Hz)\n'
 411                     'and returns the variance in X units (X^2)\n'
 412                     '<input-file> should be 1 column ASCII without a header line.\n'
 413                     'e.g: "<deflection bits>\\n..."\n')
 414     parser = OptionParser(usage=usage_string, version='%prog 0.1')
 415     parser.add_option('-f', '--sample-freq', dest='freq', default=100e3,
 416                       help='sample frequency in Hz (default %default)',
 417                       type='float', metavar='FREQ')
 418     parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
 419                       help='minimum frequency in Hz (default %default)',
 420                       type='float', metavar='FREQ')
 421     parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
 422                       help='maximum frequency in Hz (default %default)',
 423                       type='float', metavar='FREQ')
 424     parser.add_option('-o', '--output-file', dest='ofilename',
 425                       help='write output to FILE (default stdout)',
 426                       type='string', metavar='FILE')
 427     parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
 428                       help='Output comma-seperated values (default %default)',
 429                       default=False)
 430     parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
 431                       help='Print gnuplot fit check script to stderr',
 432                       default=False)
 433     parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
 434                       help='Produce pylab fit checks during execution',
 435                       default=False)
 436     parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
 437                       help='Produce plots of the pre-fit Lorentzian',
 438                       default=False)
 439     parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
 440                       help='Run in tweak-file mode',
 441                       default=False)
 442     parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
 443                       help='Read input files with datalogger.read_dict_of_arrays().  This is useful, for example, to test a single line from a tweakfile.',
 444                       default=False)
 445     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
 446                       help='Print lots of debugging information',
 447                       default=False)
 448
 449     options,args = parser.parse_args()
 450     parser.destroy()
 451     assert len(args) >= 1, "Need an input file"
 452
 453     ifilename = args[0]
 454
 455     if options.ofilename != None :
 456         ofile = file(options.ofilename, 'w')
 457     else :
 458         ofile = sys.stdout
 459     if options.verbose == True :
 460         vfile = sys.stderr
 461     else :
 462         vfile = None
 463     config.TEXT_VERBOSE = options.verbose
 464     config.PYLAB_INTERACTIVE = False
 465     config.PYLAB_VERBOSE = options.pylab
 466     config.GNUPLOT_VERBOSE = options.gnuplot
 467     PLOT_GUESSED_LORENTZIAN = options.plot_guess
 468
 469     if options.tweakmode == False :
 470         if options.datalogger_mode:
 471             data = vib_load(ifilename)
 472             deflection_bits = data['Deflection input']
 473         else:
 474             deflection_bits = read_data(ifilename)
 475         Vphoto_var = vib_analyze(deflection_bits, freq=options.freq,
 476                                  minFreq=options.min_freq,
 477                                  maxFreq=options.max_freq)
 478         print >> ofile, Vphoto_var
 479     else : # tweak mode
 480         Vphoto_vars = vib_load_analyze_tweaked(ifilename,
 481                                                minFreq=options.min_freq,
 482                                                maxFreq=options.max_freq)
 483         if options.comma_out :
 484             sep = ','
 485         else :
 486             sep = '\n'
 487         common.write_array(ofile, Vphoto_vars, sep)
 488
 489     if common._final_flush_plot != None:
 490         common._final_flush_plot()
 491
 492     if options.ofilename != None :
 493         ofile.close()
 494