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 common # common module for the calibcant package
38 import config # config module for the calibcant package
42 def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
43 zeroVphoto_bits=config.zeroVphoto_bits) :
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
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.
53 Vphoto_std_bits = deflection_bits.std()
54 Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
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) :
61 Calculate the variance in the raw data due to the cantilever's
62 thermal oscillation and convert it to Volts**2.
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).
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
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.
77 Vphoto_t = numpy.zeros((len(deflection_bits),),
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])
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)
92 A,B,C = fit_psd(freq_axis, power, minFreq, maxFreq)
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)
98 vib_plot(deflection_bits, freq_axis, power, A, B, C, plotVerbose=plotVerbose)
100 # Our A is in uV**2, so convert back to Volts**2
101 return lorentzian_area(A,B,C) * 1e-12
103 def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
105 freq_axis : array of frequencies in Hz
106 psd_data : array of PSD amplitudes for the frequencies in freq_axis
108 Calculate the variance in the raw data due to the cantilever's
109 thermal oscillation and convert it to Volts**2.
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).
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().
120 # cut out the relevent frequency range
121 if config.TEXT_VERBOSE :
122 print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
124 while freq_axis[imin] < minFreq : imin += 1
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)
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]))
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 :
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
149 # C = (2 k_B T B) / (pi m)
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)
157 # guess Q = 5, and adjust from there
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)
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)
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)
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)
211 def gnuplot_check_fit(datafilename, A, B, C) :
213 return a string containing a gnuplot script displaying the fit.
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
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
229 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
230 log_name="vib_%gHz" % freq)
231 log.write_dict_of_arrays(data)
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
240 dl = data_logger.data_load()
241 data = dl.read_dict_of_arrays(datafile)
244 def vib_plot(deflection_bits, freq_axis, power, A, B, C,
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).
252 if plotVerbose or config.PYLAB_VERBOSE :
254 common._import_pylab()
255 common._pylab.figure(config.BASE_FIGNUM+2)
256 common._pylab.hold(False)
259 common._pylab.subplot(311)
260 common._pylab.plot(deflection_bits, 'r.')
261 common._pylab.title("free oscillation")
263 # plot histogram distribution and gaussian fit
264 common._pylab.subplot(312)
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()
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)
280 common._pylab.subplot(313)
281 common._pylab.semilogy(freq_axis, power, 'r.-')
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.
290 for i in range(imin, imax) :
292 fd.write("%g\t%g\n" % (freq_axis[i],
293 C / ((A**2-x**2)**2 + (x*B)**2) ) )
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'")
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)
308 g("pause mouse 'click with mouse'")
309 g.getvar("MOUSE_BUTTON")
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,
318 Read the vib files listed in tweak_file.
319 Return an array of Vphoto variances (due to the cantilever) in Volts**2
321 dl = data_logger.data_load()
323 for path in file(tweak_file, 'r') :
325 path = path.split('\n')[0]
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)
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)
338 # commandline interface functions
341 def read_data(ifile):
343 ifile can be a filename string or open (seekable) file object.
344 returns an array of data values (1 column)
346 #returns (column 1 array, column 2 array)
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])
352 if __name__ == '__main__' :
353 # command line interface
354 from optparse import OptionParser
356 usage_string = ('%prog <input-file>\n'
357 '2008, W. Trevor King.\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)',
380 parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
381 help='Print gnuplot fit check script to stderr',
383 parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
384 help='Produce pylab fit checks during execution',
386 parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
387 help='Run in tweak-file mode',
389 parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
390 help='Print lots of debugging information',
393 options,args = parser.parse_args()
395 assert len(args) >= 1, "Need an input file"
399 if options.ofilename != None :
400 ofile = file(options.ofilename, 'w')
403 if options.verbose == True :
407 config.TEXT_VERBOSE = options.verbose
408 config.PYLAB_INTERACTIVE = False
409 config.PYLAB_VERBOSE = options.pylab
410 config.GNUPLOT_VERBOSE = options.gnuplot
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
420 Vphoto_vars = vib_load_analyze_tweaked(ifilename,
421 minFreq=options.min_freq,
422 maxFreq=options.max_freq)
423 if options.comma_out :
427 common.write_array(ofile, Vphoto_vars, sep)
429 if common._final_flush_plot != None:
430 common._final_flush_plot()
432 if options.ofilename != None :