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,
99 plotVerbose=plotVerbose)
101 # Our A is in uV**2, so convert back to Volts**2
102 return lorentzian_area(A,B,C) * 1e-12
104 def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
106 freq_axis : array of frequencies in Hz
107 psd_data : array of PSD amplitudes for the frequencies in freq_axis
109 Calculate the variance in the raw data due to the cantilever's
110 thermal oscillation and convert it to Volts**2.
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).
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().
121 # cut out the relevent frequency range
122 if config.TEXT_VERBOSE :
123 print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
125 while freq_axis[imin] < minFreq : imin += 1
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)
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]))
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 :
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
150 # C = (2 k_B T B) / (pi m)
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)
158 # guess Q = 5, and adjust from there
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)
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)
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)
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)
212 def gnuplot_check_fit(datafilename, A, B, C) :
214 return a string containing a gnuplot script displaying the fit.
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
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
230 log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
231 log_name="vib_%gHz" % freq)
232 log.write_dict_of_arrays(data)
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
241 dl = data_logger.data_load()
242 data = dl.read_dict_of_arrays(datafile)
245 def vib_plot(deflection_bits, freq_axis, power, A, B, C,
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).
253 if plotVerbose or config.PYLAB_VERBOSE :
255 common._import_pylab()
256 common._pylab.figure(config.BASE_FIGNUM+2)
259 common._pylab.subplot(311)
260 common._pylab.hold(False)
261 common._pylab.plot(deflection_bits, 'r.')
262 common._pylab.title("free oscillation")
264 # plot histogram distribution and gaussian fit
265 common._pylab.subplot(312)
266 common._pylab.hold(False)
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()
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)
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])
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.
301 for i in range(imin, imax) :
303 fd.write("%g\t%g\n" % (freq_axis[i],
304 C / ((A**2-x**2)**2 + (x*B)**2) ) )
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'")
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)
319 g("pause mouse 'click with mouse'")
320 g.getvar("MOUSE_BUTTON")
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,
329 Read the vib files listed in tweak_file.
330 Return an array of Vphoto variances (due to the cantilever) in Volts**2
332 dl = data_logger.data_load()
334 for path in file(tweak_file, 'r') :
336 path = path.split('\n')[0]
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)
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)
349 # commandline interface functions
352 def read_data(ifile):
354 ifile can be a filename string or open (seekable) file object.
355 returns an array of data values (1 column)
357 #returns (column 1 array, column 2 array)
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])
363 if __name__ == '__main__' :
364 # command line interface
365 from optparse import OptionParser
367 usage_string = ('%prog <input-file>\n'
368 '2008, W. Trevor King.\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)',
391 parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
392 help='Print gnuplot fit check script to stderr',
394 parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
395 help='Produce pylab fit checks during execution',
397 parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
398 help='Run in tweak-file mode',
400 parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
401 help='Print lots of debugging information',
404 options,args = parser.parse_args()
406 assert len(args) >= 1, "Need an input file"
410 if options.ofilename != None :
411 ofile = file(options.ofilename, 'w')
414 if options.verbose == True :
418 config.TEXT_VERBOSE = options.verbose
419 config.PYLAB_INTERACTIVE = False
420 config.PYLAB_VERBOSE = options.pylab
421 config.GNUPLOT_VERBOSE = options.gnuplot
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
431 Vphoto_vars = vib_load_analyze_tweaked(ifilename,
432 minFreq=options.min_freq,
433 maxFreq=options.max_freq)
434 if options.comma_out :
438 common.write_array(ofile, Vphoto_vars, sep)
440 if common._final_flush_plot != None:
441 common._final_flush_plot()
443 if options.ofilename != None :