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