1 # calibcant - tools for thermally calibrating AFM cantilevers
3 # Copyright (C) 2008-2012 W. Trevor King <wking@drexel.edu>
5 # This file is part of calibcant.
7 # calibcant is free software: you can redistribute it and/or modify it under
8 # the terms of the GNU General Public License as published by the Free Software
9 # Foundation, either version 3 of the License, or (at your option) any later
12 # calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14 # A PARTICULAR PURPOSE. See the GNU General Public License for more details.
16 # You should have received a copy of the GNU General Public License along with
17 # calibcant. If not, see <http://www.gnu.org/licenses/>.
19 """Thermal vibration analysis.
21 Separate the more general `vib_analyze()` from the other `vib_*()`
22 functions in calibcant.
24 The relevent physical quantities are :
25 Vphoto The photodiode vertical deflection voltage (what we measure)
28 >>> from pprint import pprint
32 >>> from .config import VibrationConfig
33 >>> from h5config.storage.hdf5 import pprint_HDF5
34 >>> from pypiezo.test import get_piezo_config
35 >>> from pypiezo.base import convert_volts_to_bits
37 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
40 >>> piezo_config = get_piezo_config()
41 >>> vibration_config = VibrationConfig()
42 >>> vibration_config['frequency'] = 50e3
44 We'll be generating a test vibration time series with the following
45 parameters. Make sure these are all floats to avoid accidental
46 integer division in later steps.
49 >>> gamma = 1.6e-6 # N*s/m
51 >>> T = 1/vibration_config['frequency']
52 >>> T # doctest: +ELLIPSIS
54 >>> N = int(2**15) # count
57 where `T` is the sampling period, `N` is the number of samples, and
58 `F_sigma` is the standard deviation of the white-noise external force.
59 Note that the resonant frequency is less than the Nyquist frequency so
60 we don't have to worry too much about aliasing.
62 >>> w0 = numpy.sqrt(k/m)
63 >>> f0 = w0/(2*numpy.pi)
64 >>> f0 # doctest: +ELLIPSIS
66 >>> f_nyquist = vibration_config['frequency']/2
67 >>> f_nyquist # doctest: +ELLIPSIS
72 >>> damping = gamma / (2*m*w0)
73 >>> damping # doctest: +ELLIPSIS
79 >>> Q # doctest: +ELLIPSIS
81 >>> (1 / (2*damping)) / Q # doctest: +ELLIPSIS
84 We expect the white-noise power spectral density (PSD) to be a flat
87 >>> F0 = F_sigma**2 * 2 * T
89 because the integral from `0` `1/2T` should be `F_sigma**2`.
91 The expected time series PSD parameters are
94 >>> B = gamma/(m*2*numpy.pi)
95 >>> C = F0/(m**2*(2*numpy.pi)**4)
97 Simulate a time series with the proper PSD using center-differencing.
99 m\ddot{x} + \gamma \dot{x} + kx = F
101 m \frac{x_{i+1} - 2x_i + x_{i-1}}{T**2}
102 + \gamma \frac{x_{i+1}-x_{i-1}}{T}
105 a x_{i+1} + b x_{i} + c x_{i-1} = F_i
107 where `T` is the sampling period, `i=t/T` is the measurement index,
108 `a=m/T**2+gamma/2T`, `b=-2m/T**2+k`, and `c=m/T**2-gamma/2T`.
109 Rearranging and shifting to `j=i-1`
111 x_j = \frac{F_{i-1} - bx_{i-1} - cx_{i-2}}{a}
113 >>> a = m/T**2 + gamma/(2*T)
114 >>> b = -2*m/T**2 + k
115 >>> c = m/T**2 - gamma/(2*T)
116 >>> x = numpy.zeros((N+2,), dtype=numpy.float) # two extra initial points
117 >>> F = numpy.zeros((N,), dtype=numpy.float)
118 >>> for i in range(2, x.size):
119 ... Fp = random.gauss(mu=0, sigma=F_sigma)
123 ... x[i] = (Fp - b*xp - c*xpp)/a
124 >>> x = x[2:] # drop the initial points
126 Convert the simulated data to bits.
129 >>> deflection_bits = convert_volts_to_bits(
130 ... piezo_config.select_config('inputs', 'deflection'), x)
132 Analyze the simulated data.
134 >>> naive_vibration = vib_analyze_naive(deflection)
135 >>> naive_vibration # doctest: +SKIP
137 >>> abs(naive_vibration / 136.9e6 - 1) < 0.1
140 >>> processed_vibration = vib_analyze(
141 ... deflection_bits, vibration_config,
142 ... piezo_config.select_config('inputs', 'deflection'))
143 >>> processed_vibration # doctest: +SKIP
146 >>> vib_plot(deflection=deflection_bits, vibration_config=vibration_config)
147 >>> vib_save(filename=filename, group='/vibration/',
148 ... raw_vibration=deflection_bits, vibration_config=vibration_config,
149 ... deflection_channel_config=piezo_config.select_config(
150 ... 'inputs', 'deflection'),
151 ... processed_vibration=processed_vibration)
153 >>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
157 /vibration/config/deflection
158 <HDF5 dataset "channel": shape (), type "<i4">
160 <HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
162 <HDF5 dataset "conversion-origin": shape (), type "<i4">
164 <HDF5 dataset "device": shape (), type "|S12">
166 <HDF5 dataset "inverse-conversion-coefficients": shape (2,), type "<i4">
168 <HDF5 dataset "inverse-conversion-origin": shape (), type "<i4">
170 <HDF5 dataset "maxdata": shape (), type "<i4">
172 <HDF5 dataset "name": shape (), type "|S10">
174 <HDF5 dataset "range": shape (), type "<i4">
176 <HDF5 dataset "subdevice": shape (), type "<i4">
178 /vibration/config/vibration
179 <HDF5 dataset "chunk-size": shape (), type "<i4">
181 <HDF5 dataset "frequency": shape (), type "<f8">
183 <HDF5 dataset "maximum-fit-frequency": shape (), type "<f8">
185 <HDF5 dataset "minimum-fit-frequency": shape (), type "<f8">
187 <HDF5 dataset "model": shape (), type "|S12">
189 <HDF5 dataset "overlap": shape (), type "|b1">
191 <HDF5 dataset "sample-time": shape (), type "<i4">
193 <HDF5 dataset "window": shape (), type "|S4">
195 <HDF5 dataset "processed": shape (), type "<f8">
198 <HDF5 dataset "deflection": shape (32768,), type "<f8">
201 >>> (raw_vibration,vibration_config,deflection_channel_config,
202 ... processed_vibration) = vib_load(
203 ... filename=filename, group='/vibration/')
205 >>> processed_vibration # doctest: +SKIP
207 >>> abs(processed_vibration / 136.5e6 - 1) < 0.1
210 >>> os.remove(filename)
217 import numpy as _numpy
218 from scipy.optimize import leastsq as _leastsq
221 import matplotlib as _matplotlib
222 import matplotlib.pyplot as _matplotlib_pyplot
223 import time as _time # for timestamping lines on plots
224 except (ImportError, RuntimeError), e:
226 _matplotlib_import_error = e
228 from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
229 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
230 import FFT_tools as _FFT_tools
231 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
232 from pypiezo.config import ChannelConfig as _ChannelConfig
234 from . import LOG as _LOG
235 from . import package_config as _package_config
236 from .config import Variance as _Variance
237 from .config import BreitWigner as _BreitWigner
238 from .config import OffsetBreitWigner as _OffsetBreitWigner
239 from .config import VibrationConfig as _VibrationConfig
242 def vib_analyze_naive(deflection):
243 """Calculate the deflection variance in Volts**2.
245 This method is simple and easy to understand, but it highly
246 succeptible to noise, drift, etc.
249 deflection : numpy array with deflection timeseries in Volts.
251 std = deflection.std()
253 _LOG.debug('naive deflection variance: %g V**2' % var)
256 def vib_analyze(deflection, vibration_config, deflection_channel_config,
258 """Calculate the deflection variance in Volts**2.
260 Improves on `vib_analyze_naive()` by first converting `Vphoto(t)`
261 to frequency space, and fitting a Breit-Wigner in the relevant
262 frequency range (see cantilever_calib.pdf for derivation).
263 However, there may be cases where the fit is thrown off by noise
264 spikes in frequency space. To protect from errors, the fitted
265 variance is compared to the naive variance (where all noise is
266 included), and the minimum variance is returned.
269 deflection Vphoto deflection input in bits.
270 vibration_config `.config._VibrationConfig` instance
271 deflection_channel_config
272 deflection `pypiezo.config.ChannelConfig` instance
273 plot boolean overriding matplotlib config setting.
275 The conversion to frequency space generates an average power
276 spectrum by breaking the data into windowed chunks and averaging
277 the power spectrums for the chunks together. See
278 `FFT_tools.unitary_avg_power_spectrum()` for details.
280 # convert the data from bits to volts
281 deflection_v = _convert_bits_to_volts(
282 deflection_channel_config, deflection)
283 mean = deflection_v.mean()
284 _LOG.debug('average thermal deflection (Volts): %g' % mean)
286 naive_variance = vib_analyze_naive(deflection_v)
287 if vibration_config['model'] == _Variance:
288 return naive_variance
290 # Compute the average power spectral density per unit time (in uV**2/Hz)
291 _LOG.debug('compute the averaged power spectral density in uV**2/Hz')
292 freq_axis,power = _FFT_tools.unitary_avg_power_spectrum(
293 (deflection_v - mean)*1e6, vibration_config['frequency'],
294 vibration_config['chunk-size'], vibration_config['overlap'],
295 vibration_config['window'])
299 min_frequency=vibration_config['minimum-fit-frequency'],
300 max_frequency=vibration_config['maximum-fit-frequency'],
301 offset=vibration_config['model'] == _OffsetBreitWigner)
303 _LOG.debug('fit PSD(f) = C / ((A**2-f**2)**2 + (f*B)**2) with '
304 'A = %g, B = %g, C = %g, D = %g' % (A, B, C, D))
306 if plot or _package_config['matplotlib']:
307 vib_plot(deflection, freq_axis, power, A, B, C, D,
308 vibration_config=vibration_config)
310 # Our A is in uV**2, so convert back to Volts**2
311 fitted_variance = breit_wigner_area(A,B,C) * 1e-12
313 _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
315 if _package_config['matplotlib']:
316 vib_plot(deflection, freq_axis, power, A, B, C, D,
317 vibration_config=vibration_config)
319 return min(fitted_variance, naive_variance)
321 def breit_wigner(f, A, B, C, D=0):
322 """Breit-Wigner (sortof).
329 D Optional white-noise offset
331 All parameters must be postive.
333 return abs(C) / ((A**2-f**2)**2 + (B*f)**2) + abs(D)
335 def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
337 """Fit the FFTed vibration data to a Breit-Wigner.
340 freq_axis array of frequencies in Hz
341 psd_data array of PSD amplitudes for the frequencies in freq_axis
342 min_frequency lower bound of Breit-Wigner fitting region
343 max_frequency upper bound of Breit-Wigner fitting region
344 offset add a white-noise offset to the Breit-Wigner model
346 Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
348 # cut out the relevent frequency range
349 _LOG.debug('cut the frequency range %g to %g Hz'
350 % (min_frequency, max_frequency))
352 while freq_axis[imin] < min_frequency : imin += 1
354 while freq_axis[imax] < max_frequency : imax += 1
355 assert imax >= imin + 10 , 'less than 10 points in freq range (%g,%g)' % (
356 min_frequency, max_frequency)
358 # generate guesses for Breit-Wigner parameters A, B, C, and D
359 max_psd_index = _numpy.argmax(psd_data[imin:imax]) + imin
360 max_psd = psd_data[max_psd_index]
361 res_freq = freq_axis[max_psd_index]
363 # Breit-Wigner L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
364 # is expected power spectrum for
365 # x'' + B x' + A^2 x'' = F_external(t)/m
367 # C = (2 k_B T B) / (pi m)
369 # A = resonant frequency
370 # peak at x_res = sqrt(A^2 - B^2/2) (by differentiating)
371 # which could be complex if there isn't a peak (overdamped)
372 # peak height = C / (3 x_res^4 + A^4)
376 # Guessing Q = 1 is pretty safe.
380 # so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
381 # B = x_res / sqrt(Q^2-1/2)
382 B_guess = res_freq / _numpy.sqrt(Q_guess**2-0.5)
383 A_guess = Q_guess*B_guess
384 C_guess = max_psd * (-res_freq**4 + A_guess**4)
386 D_guess = psd_data[max_psd_index]
390 _LOG.debug(('guessed params: resonant freq %g, max psd %g, Q %g, A %g, '
391 'B %g, C %g, D %g') % (
392 res_freq, max_psd, Q_guess, A_guess, B_guess, C_guess, D_guess))
393 # Half width w on lower side when L(a-w) = L(a)/2
394 # (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
395 # Let W=(a-w)**2, A=a**2, and B=b**2
396 # (A - W)**2 + BW = 2*AB
397 # W**2 - 2AW + A**2 + BW = 2AB
398 # W**2 + (B-2A)W + (A**2-2AB) = 0
399 # W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
400 # = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
401 # (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
402 # so w is a disaster ;)
403 # For some values of A and B (non-underdamped), W is imaginary or negative.
405 # The mass m is given by m = k_B T / (pi A)
406 # The spring constant k is given by k = m (omega_0)**2
407 # The quality factor Q is given by Q = omega_0 m / gamma
409 # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
411 # fit Breit-Wigner using scipy.optimize.leastsq
412 def residual(p, y, x):
413 return breit_wigner(x, *p) - y
415 guess = _numpy.array((A_guess, B_guess, C_guess, D_guess))
417 guess = _numpy.array((A_guess, B_guess, C_guess))
419 p,cov,info,mesg,ier = _leastsq(
421 args=(psd_data[imin:imax], freq_axis[imin:imax]),
422 full_output=True, maxfev=10000)
423 _LOG.debug('fitted params: %s' % p)
424 _LOG.debug('covariance matrix: %s' % cov)
425 #_LOG.debug('info: %s' % info)
426 _LOG.debug('message: %s' % mesg)
428 _LOG.debug('solution converged')
430 _LOG.debug('solution did not converge')
436 A=abs(A) # A and B only show up as squares in f(x)
437 B=abs(B) # so ensure we get positive values.
438 C=abs(C) # Only abs(C) is used in breit_wigner().
441 def breit_wigner_area(A, B, C):
442 # Integrating the the power spectral density per unit time (power) over the
443 # frequency axis [0, fN] returns the total signal power per unit time
444 # int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
445 # = <V(t)**2>, the variance for our AC signal.
446 # The variance from our fitted Breit-Wigner is the area under the Breit-Wigner
447 # <V(t)**2> = (pi*C) / (2*B*A**2)
448 return (_numpy.pi * C) / (2 * B * A**2)
450 def vib_save(filename, group='/', raw_vibration=None, vibration_config=None,
451 deflection_channel_config=None, processed_vibration=None):
452 with _h5py.File(filename, 'a') as f:
453 cwg = _h5_create_group(f, group)
454 if raw_vibration is not None:
456 del cwg['raw/deflection']
459 cwg['raw/deflection'] = raw_vibration
460 storage = _HDF5_Storage()
461 for config,key in [(vibration_config, 'config/vibration'),
462 (deflection_channel_config,
463 'config/deflection')]:
466 config_cwg = _h5_create_group(cwg, key)
467 storage.save(config=config, group=config_cwg)
468 if processed_vibration is not None:
473 cwg['processed'] = processed_vibration
475 def vib_load(filename, group='/'):
476 assert group.endswith('/')
477 raw_vibration = processed_vibration = None
479 with _h5py.File(filename, 'a') as f:
481 raw_vibration = f[group+'raw/deflection'][...]
484 for Config,key in [(_VibrationConfig, 'config/vibration'),
485 (_ChannelConfig, 'config/deflection')]:
486 config = Config(storage=_HDF5_Storage(
487 filename=filename, group=group+key))
488 configs.append(config)
490 processed_vibration = float(f[group+'processed'][...])
493 ret = [raw_vibration]
495 ret.append(processed_vibration)
496 for config in configs:
500 def vib_plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
501 C=None, D=0, vibration_config=None, analyze=False):
502 """Plot 3 subfigures displaying vibration data and analysis.
504 Time series (Vphoto vs sample index) (show any major drift events),
505 A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
506 FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
509 raise _matplotlib_import_error
510 figure = _matplotlib_pyplot.figure()
513 assert deflection != None, (
514 'must set at least one of `deflection` or `power`.')
515 time_axes = figure.add_subplot(2, 1, 1)
516 hist_axes = figure.add_subplot(2, 1, 2)
518 elif deflection is None:
519 time_axes = hist_axes = None
520 freq_axes = figure.add_subplot(1, 1, 1)
522 time_axes = figure.add_subplot(3, 1, 1)
523 hist_axes = figure.add_subplot(3, 1, 2)
524 freq_axes = figure.add_subplot(3, 1, 3)
526 timestamp = _time.strftime('%H%M%S')
528 if deflection is not None:
529 time_axes.plot(deflection, 'r.')
530 time_axes.set_title('free oscillation')
532 # plot histogram distribution and gaussian fit
534 n,bins,patches = hist_axes.hist(
535 deflection, bins=30, normed=True, align='mid')
536 gauss = _numpy.zeros((len(bins),), dtype=_numpy.float)
537 mean = deflection.mean()
538 std = deflection.std()
541 gauss = _numpy.sqrt(2*pi)/std * exp(-0.5*((bins-mean)/std)**2)
542 # Matplotlib's normed histogram uses bin heights of n/(len(x)*dbin)
543 dbin = bins[1]-bins[0]
544 hist_axes.plot(bins, gauss/dbin, 'r-')
545 if power is not None:
547 freq_axes.set_yscale('log')
548 freq_axes.plot(freq_axis, power, 'r.-')
549 xmin,xmax = freq_axes.get_xbound()
550 ymin,ymax = freq_axes.get_ybound()
552 # highlight the region we're fitting
554 vibration_config['minimum-fit-frequency'],
555 vibration_config['maximum-fit-frequency'],
556 facecolor='g', alpha=0.1, zorder=-2)
559 fitdata = breit_wigner(freq_axis, A, B, C, D)
560 freq_axes.plot(freq_axis, fitdata, 'b-')
561 noisefloor = D + 0*freq_axis;
562 freq_axes.plot(freq_axis, noisefloor)
565 res_freq = _numpy.sqrt(A**2 - B**2/2)
566 freq_axes.axvline(res_freq, color='b', zorder=-1)
568 freq_axes.set_title('power spectral density %s' % timestamp)
569 freq_axes.axis([xmin,xmax,ymin,ymax])
570 freq_axes.set_xlabel('frequency (Hz)')
571 freq_axes.set_ylabel('PSD')
572 if hasattr(figure, 'show'):