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