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