3 # calibcant - tools for thermally calibrating AFM cantilevers
5 # Copyright (C) 2007,2008, William Trevor King
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.
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.
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
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.
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.
31 The relevent physical quantities are :
32 Vphoto The photodiode vertical deflection voltage (what we measure)
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
43 PLOT_GUESSED_LORENTZIAN=False
45 def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
46 zeroVphoto_bits=config.zeroVphoto_bits) :
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
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.
56 Vphoto_std_bits = deflection_bits.std()
57 Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
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) :
64 Calculate the variance in the raw data due to the cantilever's
65 thermal oscillation and convert it to Volts**2.
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).
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
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.
80 Vphoto_t = numpy.zeros((len(deflection_bits),),
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])
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)
95 A,B,C = fit_psd(freq_axis, power, minFreq, maxFreq)
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)
101 vib_plot(deflection_bits, freq_axis, power, A, B, C,
102 minFreq, maxFreq, plotVerbose=plotVerbose)
104 # Our A is in uV**2, so convert back to Volts**2
105 return lorentzian_area(A,B,C) * 1e-12
107 def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
109 freq_axis : array of frequencies in Hz
110 psd_data : array of PSD amplitudes for the frequencies in freq_axis
112 Calculate the variance in the raw data due to the cantilever's
113 thermal oscillation and convert it to Volts**2.
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).
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().
124 # cut out the relevent frequency range
125 if config.TEXT_VERBOSE :
126 print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
128 while freq_axis[imin] < minFreq : imin += 1
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)
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 :
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
146 # C = (2 k_B T B) / (pi m)
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)
154 # guess Q = 1, and adjust from there
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)
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.
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)
185 # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
187 if False: # Gnuplot fit worse than scipy.optimize.leastsq fit.
188 # fit Lorentzian using Gnuplot's 'fit' command
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]))
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'))
205 os.remove(datafilename)
207 # fit Lorenzian using scipy.optimize.leastsq
208 def residual(p, y, x):
212 Y = C / ((A**2-x**2)**2 + (B*x)**2)
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
225 print "Solution converged"
227 print "Solution did not converge"
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)
240 def gnuplot_check_fit(datafilename, A, B, C) :
242 return a string containing a gnuplot script displaying the fit.
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
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
258 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
259 log_name="vib_%gHz" % freq)
260 log.write_dict_of_arrays(data)
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
269 dl = data_logger.data_load()
270 data = dl.read_dict_of_arrays(datafile)
273 def vib_plot(deflection_bits, freq_axis, power, A, B, C,
274 minFreq=None, maxFreq=None, plotVerbose=True) :
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).
281 if plotVerbose or config.PYLAB_VERBOSE :
282 common._import_pylab()
283 common._pylab.figure(config.BASE_FIGNUM+2)
285 if deflection_bits != None:
287 common._pylab.subplot(311)
288 common._pylab.hold(False)
289 common._pylab.plot(deflection_bits, 'r.')
290 common._pylab.title("free oscillation")
292 # plot histogram distribution and gaussian fit
293 common._pylab.subplot(312)
294 common._pylab.hold(False)
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()
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)
310 axes = common._pylab.subplot(313)
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()
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,
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])
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.
338 for i in range(imin, imax) :
340 fd.write("%g\t%g\n" % (freq_axis[i],
341 C / ((A**2-x**2)**2 + (x*B)**2) ) )
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'")
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)
356 g("pause mouse 'click with mouse'")
357 g.getvar("MOUSE_BUTTON")
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,
366 Read the vib files listed in tweak_file.
367 Return an array of Vphoto variances (due to the cantilever) in Volts**2
369 dl = data_logger.data_load()
371 for path in file(tweak_file, 'r') :
373 path = path.split('\n')[0]
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)
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)
386 # commandline interface functions
389 def read_data(ifile):
391 ifile can be a filename string or open (seekable) file object.
392 returns an array of data values (1 column)
394 #returns (column 1 array, column 2 array)
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])
400 if __name__ == '__main__' :
401 # command line interface
402 from optparse import OptionParser
404 usage_string = ('%prog <input-file>\n'
405 '2008, W. Trevor King.\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)',
427 parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
428 help='Print gnuplot fit check script to stderr',
430 parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
431 help='Produce pylab fit checks during execution',
433 parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
434 help='Produce plots of the pre-fit Lorentzian',
436 parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
437 help='Run in tweak-file mode',
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.',
442 parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
443 help='Print lots of debugging information',
446 options,args = parser.parse_args()
448 assert len(args) >= 1, "Need an input file"
452 if options.ofilename != None :
453 ofile = file(options.ofilename, 'w')
456 if options.verbose == True :
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
466 if options.tweakmode == False :
467 if options.datalogger_mode:
468 data = vib_load(ifilename)
469 deflection_bits = data['Deflection input']
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
477 Vphoto_vars = vib_load_analyze_tweaked(ifilename,
478 minFreq=options.min_freq,
479 maxFreq=options.max_freq)
480 if options.comma_out :
484 common.write_array(ofile, Vphoto_vars, sep)
486 if common._final_flush_plot != None:
487 common._final_flush_plot()
489 if options.ofilename != None :