calibcant/vibration_analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2012 W. Trevor King <wking@drexel.edu>
   4 #
   5 # This file is part of calibcant.
   6 #
   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
  10 # version.
  11 #
  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.
  15 #
  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/>.
  18
  19 """Thermal vibration analysis.
  20
  21 Separate the more general `analyze()` from the other `vibration`
  22 functions in calibcant.
  23
  24 The relevent physical quantities are :
  25   Vphoto   The photodiode vertical deflection voltage (what we measure)
  26
  27 >>> import os
  28 >>> from pprint import pprint
  29 >>> import random
  30 >>> import tempfile
  31 >>> import numpy
  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
  36
  37 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
  38 >>> os.close(fd)
  39
  40 >>> piezo_config = get_piezo_config()
  41 >>> config = VibrationConfig()
  42 >>> config['frequency'] = 50e3
  43
  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.
  47
  48 >>> m = 5e-11       # kg
  49 >>> gamma = 1.6e-6  # N*s/m
  50 >>> k = 0.05        # N/m
  51 >>> T = 1/config['frequency']
  52 >>> T  # doctest: +ELLIPSIS
  53 2...e-05
  54 >>> N = int(2**15)  # count
  55 >>> F_sigma = 1e3   # N
  56
  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.
  61
  62 >>> w0 = numpy.sqrt(k/m)
  63 >>> f0 = w0/(2*numpy.pi)
  64 >>> f0  # doctest: +ELLIPSIS
  65 5032.9...
  66 >>> f_nyquist = config['frequency']/2
  67 >>> f_nyquist  # doctest: +ELLIPSIS
  68 25000.0...
  69
  70 The damping ratio is
  71
  72 >>> damping = gamma / (2*m*w0)
  73 >>> damping  # doctest: +ELLIPSIS
  74 0.505...
  75
  76 The quality factor is
  77
  78 >>> Q = m*w0 / gamma
  79 >>> Q  # doctest: +ELLIPSIS
  80 0.988...
  81 >>> (1 / (2*damping)) / Q  # doctest: +ELLIPSIS
  82 1.000...
  83
  84 We expect the white-noise power spectral density (PSD) to be a flat
  85 line at
  86
  87 >>> F0 = F_sigma**2 * 2 * T
  88
  89 because the integral from `0` `1/2T` should be `F_sigma**2`.
  90
  91 The expected time series PSD parameters are
  92
  93 >>> A = f0
  94 >>> B = gamma/(m*2*numpy.pi)
  95 >>> C = F0/(m**2*(2*numpy.pi)**4)
  96
  97 Simulate a time series with the proper PSD using center-differencing.
  98
  99   m\ddot{x} + \gamma \dot{x} + kx = F
 100
 101   m \frac{x_{i+1} - 2x_i + x_{i-1}}{T**2}
 102     + \gamma \frac{x_{i+1}-x_{i-1}}{T}
 103     + kx_i = F_i
 104
 105   a x_{i+1} + b x_{i} + c x_{i-1} = F_i
 106
 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`
 110
 111   x_j = \frac{F_{i-1} - bx_{i-1} - cx_{i-2}}{a}
 112
 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)
 120 ...     F[i-2] = Fp
 121 ...     xp = x[i-1]
 122 ...     xpp = x[i-2]
 123 ...     x[i] = (Fp - b*xp - c*xpp)/a
 124 >>> x = x[2:]  # drop the initial points
 125
 126 Convert the simulated data to bits.
 127
 128 >>> deflection = x
 129 >>> deflection_bits = convert_volts_to_bits(
 130 ...     piezo_config.select_config('inputs', 'deflection'), x)
 131
 132 Analyze the simulated data.
 133
 134 >>> naive = analyze_naive(deflection)
 135 >>> naive  # doctest: +SKIP
 136 136871517.43486854
 137 >>> abs(naive / 136.9e6 - 1) < 0.1
 138 True
 139
 140 >>> processed = analyze(
 141 ...     deflection_bits, config,
 142 ...     piezo_config.select_config('inputs', 'deflection'))
 143 >>> processed  # doctest: +SKIP
 144 136457906.25574699
 145
 146 >>> plot(deflection=deflection_bits, config=config)
 147 >>> save(filename=filename, group='/vibration/',
 148 ...     raw=deflection_bits, config=config,
 149 ...     deflection_channel_config=piezo_config.select_config(
 150 ...         'inputs', 'deflection'),
 151 ...     processed=processed)
 152
 153 >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 154 /
 155   /vibration
 156     /vibration/config
 157       /vibration/config/deflection
 158         <HDF5 dataset "analog-reference": shape (), type "|S6">
 159           ground
 160         <HDF5 dataset "channel": shape (), type "<i4">
 161           0
 162         <HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
 163           [0 1]
 164         <HDF5 dataset "conversion-origin": shape (), type "<i4">
 165           0
 166         <HDF5 dataset "device": shape (), type "|S12">
 167           /dev/comedi0
 168         <HDF5 dataset "inverse-conversion-coefficients": shape (2,), type "<i4">
 169           [0 1]
 170         <HDF5 dataset "inverse-conversion-origin": shape (), type "<i4">
 171           0
 172         <HDF5 dataset "maxdata": shape (), type "<i4">
 173           100
 174         <HDF5 dataset "name": shape (), type "|S10">
 175           deflection
 176         <HDF5 dataset "range": shape (), type "<i4">
 177           0
 178         <HDF5 dataset "subdevice": shape (), type "<i4">
 179           -1
 180       /vibration/config/vibration
 181         <HDF5 dataset "chunk-size": shape (), type "<i4">
 182           2048
 183         <HDF5 dataset "frequency": shape (), type "<f8">
 184           50000.0
 185         <HDF5 dataset "maximum-fit-frequency": shape (), type "<f8">
 186           25000.0
 187         <HDF5 dataset "minimum-fit-frequency": shape (), type "<f8">
 188           500.0
 189         <HDF5 dataset "model": shape (), type "|S12">
 190           Breit-Wigner
 191         <HDF5 dataset "overlap": shape (), type "|b1">
 192           False
 193         <HDF5 dataset "sample-time": shape (), type "<i4">
 194           1
 195         <HDF5 dataset "window": shape (), type "|S4">
 196           Hann
 197     /vibration/processed
 198       <HDF5 dataset "data": shape (), type "<f8">
 199         ...
 200       <HDF5 dataset "units": shape (), type "|S6">
 201         V^2/Hz
 202     /vibration/raw
 203       <HDF5 dataset "data": shape (32768,), type "<f8">
 204         [...]
 205       <HDF5 dataset "units": shape (), type "|S4">
 206         bits
 207
 208 >>> data = load(filename=filename, group='/vibration/')
 209
 210 >>> pprint(data)  # doctest: +ELLIPSIS
 211 {'config': {'vibration': <InputChannelConfig ...>},
 212  'processed': ...,
 213  'raw': array([...])}
 214 >>> data['processed']  # doctest: +SKIP
 215 136457906.25574699
 216 >>> abs(data['processed'] / 136.5e6 - 1) < 0.1
 217 True
 218
 219 >>> os.remove(filename)
 220 """
 221
 222 import os as _os
 223 import time as _time
 224
 225 import h5py as _h5py
 226 import numpy as _numpy
 227 from scipy.optimize import leastsq as _leastsq
 228
 229 try:
 230     import matplotlib as _matplotlib
 231     import matplotlib.pyplot as _matplotlib_pyplot
 232     import time as _time  # for timestamping lines on plots
 233 except (ImportError, RuntimeError), e:
 234     _matplotlib = None
 235     _matplotlib_import_error = e
 236
 237 from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
 238 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
 239 import FFT_tools as _FFT_tools
 240 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
 241 from pypiezo.config import InputChannelConfig as _InputChannelConfig
 242
 243 from . import LOG as _LOG
 244 from . import package_config as _package_config
 245 from .config import Variance as _Variance
 246 from .config import BreitWigner as _BreitWigner
 247 from .config import OffsetBreitWigner as _OffsetBreitWigner
 248 from .config import VibrationConfig as _VibrationConfig
 249 from .util import SaveSpec as _SaveSpec
 250 from .util import save as _save
 251 from .util import load as _load
 252
 253
 254 def analyze_naive(deflection):
 255     """Calculate the deflection variance in Volts**2.
 256
 257     This method is simple and easy to understand, but it highly
 258     succeptible to noise, drift, etc.
 259
 260     Inputs
 261       deflection : numpy array with deflection timeseries in Volts.
 262     """
 263     std = deflection.std()
 264     var = std**2
 265     _LOG.debug('naive deflection variance: %g V**2' % var)
 266     return var
 267
 268 def analyze(deflection, config, deflection_channel_config,
 269             plot=False):
 270     """Calculate the deflection variance in Volts**2.
 271
 272     Improves on `analyze_naive()` by first converting `Vphoto(t)`
 273     to frequency space, and fitting a Breit-Wigner in the relevant
 274     frequency range (see cantilever_calib.pdf for derivation).
 275     However, there may be cases where the fit is thrown off by noise
 276     spikes in frequency space.  To protect from errors, the fitted
 277     variance is compared to the naive variance (where all noise is
 278     included), and the minimum variance is returned.
 279
 280     Inputs:
 281       deflection        Vphoto deflection input in bits.
 282       config            `.config.VibrationConfig` instance
 283       deflection_channel_config
 284                         deflection `pypiezo.config.ChannelConfig` instance
 285       plot              boolean overriding matplotlib config setting.
 286
 287     The conversion to frequency space generates an average power
 288     spectrum by breaking the data into windowed chunks and averaging
 289     the power spectrums for the chunks together. See
 290     `FFT_tools.unitary_avg_power_spectrum()` for details.
 291     """
 292     # convert the data from bits to volts
 293     deflection_v = _convert_bits_to_volts(
 294         deflection_channel_config, deflection)
 295     mean = deflection_v.mean()
 296     _LOG.debug('average thermal deflection (Volts): %g' % mean)
 297
 298     naive_variance = analyze_naive(deflection_v)
 299     if config['model'] == _Variance:
 300         return naive_variance
 301
 302     # Compute the average power spectral density per unit time (in uV**2/Hz)
 303     _LOG.debug('compute the averaged power spectral density in uV**2/Hz')
 304     freq_axis,power = _FFT_tools.unitary_avg_power_spectrum(
 305         (deflection_v - mean)*1e6, config['frequency'],
 306         config['chunk-size'], config['overlap'],
 307         config['window'])
 308
 309     A,B,C,D = fit_psd(
 310         freq_axis, power,
 311         min_frequency=config['minimum-fit-frequency'],
 312         max_frequency=config['maximum-fit-frequency'],
 313         offset=config['model'] == _OffsetBreitWigner)
 314
 315     _LOG.debug('fit PSD(f) = C / ((A**2-f**2)**2 + (f*B)**2) with '
 316                'A = %g, B = %g, C = %g, D = %g' % (A, B, C, D))
 317
 318     if plot or _package_config['matplotlib']:
 319         _plot(deflection, freq_axis, power, A, B, C, D, config=config)
 320
 321     # Our A is in uV**2, so convert back to Volts**2
 322     fitted_variance = breit_wigner_area(A,B,C) * 1e-12
 323
 324     _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
 325
 326     return min(fitted_variance, naive_variance)
 327
 328 def breit_wigner(f, A, B, C, D=0):
 329     """Breit-Wigner (sortof).
 330
 331     Inputs
 332       f  Frequency
 333       A  Resonant frequency
 334       B  Quality Q = A/B
 335       C  Scaling factor
 336       D  Optional white-noise offset
 337
 338     All parameters must be postive.
 339     """
 340     return abs(C) / ((A**2-f**2)**2 + (B*f)**2) + abs(D)
 341
 342 def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
 343             offset=False):
 344     """Fit the FFTed vibration data to a Breit-Wigner.
 345
 346     Inputs
 347       freq_axis      array of frequencies in Hz
 348       psd_data       array of PSD amplitudes for the frequencies in freq_axis
 349       min_frequency  lower bound of Breit-Wigner fitting region
 350       max_frequency  upper bound of Breit-Wigner fitting region
 351       offset         add a white-noise offset to the Breit-Wigner model
 352     Output
 353       Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
 354     """
 355     # cut out the relevent frequency range
 356     _LOG.debug('cut the frequency range %g to %g Hz'
 357                % (min_frequency, max_frequency))
 358     imin = 0
 359     while freq_axis[imin] < min_frequency : imin += 1
 360     imax = imin
 361     while freq_axis[imax] < max_frequency : imax += 1
 362     assert imax >= imin + 10 , 'less than 10 points in freq range (%g,%g)' % (
 363         min_frequency, max_frequency)
 364
 365     # generate guesses for Breit-Wigner parameters A, B, C, and D
 366     max_psd_index = _numpy.argmax(psd_data[imin:imax]) + imin
 367     max_psd = psd_data[max_psd_index]
 368     res_freq = freq_axis[max_psd_index]
 369
 370     # Breit-Wigner L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
 371     # is expected power spectrum for
 372     # x'' + B x' + A^2 x'' = F_external(t)/m
 373     # (A = omega_0)
 374     # C = (2 k_B T B) / (pi m)
 375     #
 376     # A = resonant frequency
 377     # peak at  x_res = sqrt(A^2 - B^2/2)  (by differentiating)
 378     #  which could be complex if there isn't a peak (overdamped)
 379     # peak height    = C / (3 x_res^4 + A^4)
 380     #
 381     # Q = A/B
 382     #
 383     # Guessing Q = 1 is pretty safe.
 384
 385     Q_guess = 1
 386
 387     # so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
 388     #    B = x_res / sqrt(Q^2-1/2)
 389     B_guess = res_freq / _numpy.sqrt(Q_guess**2-0.5)
 390     A_guess = Q_guess*B_guess
 391     C_guess = max_psd * (-res_freq**4 + A_guess**4)
 392     if offset:
 393         D_guess = psd_data[max_psd_index]
 394         C_guess -= D_guess
 395     else:
 396         D_guess = 0
 397     _LOG.debug(('guessed params: resonant freq %g, max psd %g, Q %g, A %g, '
 398                 'B %g, C %g, D %g') % (
 399             res_freq, max_psd, Q_guess, A_guess, B_guess, C_guess, D_guess))
 400     # Half width w on lower side when L(a-w) = L(a)/2
 401     #  (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
 402     # Let W=(a-w)**2, A=a**2, and B=b**2
 403     #  (A - W)**2 + BW = 2*AB
 404     #  W**2 - 2AW + A**2 + BW = 2AB
 405     #  W**2 + (B-2A)W + (A**2-2AB) = 0
 406     #  W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
 407     #    = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 408     #  (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 409     #  so w is a disaster ;)
 410     # For some values of A and B (non-underdamped), W is imaginary or negative.
 411     #
 412     # The mass m is given by m = k_B T / (pi A)
 413     # The spring constant k is given by k = m (omega_0)**2
 414     # The quality factor Q is given by Q = omega_0 m / gamma
 415
 416     # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
 417
 418     # fit Breit-Wigner using scipy.optimize.leastsq
 419     def residual(p, y, x):
 420         return breit_wigner(x, *p) - y
 421     if offset:
 422         guess = _numpy.array((A_guess, B_guess, C_guess, D_guess))
 423     else:
 424         guess = _numpy.array((A_guess, B_guess, C_guess))
 425
 426     p,cov,info,mesg,ier = _leastsq(
 427         residual, guess,
 428         args=(psd_data[imin:imax], freq_axis[imin:imax]),
 429         full_output=True, maxfev=10000)
 430     _LOG.debug('fitted params: %s' % p)
 431     _LOG.debug('covariance matrix: %s' % cov)
 432     #_LOG.debug('info: %s' % info)
 433     _LOG.debug('message: %s' % mesg)
 434     if ier == 1:
 435         _LOG.debug('solution converged')
 436     else:
 437         _LOG.debug('solution did not converge')
 438     if offset:
 439         A,B,C,D = p
 440     else:
 441         A,B,C = p
 442         D = 0
 443     A=abs(A) # A and B only show up as squares in f(x)
 444     B=abs(B) # so ensure we get positive values.
 445     C=abs(C) # Only abs(C) is used in breit_wigner().
 446     return (A, B, C, D)
 447
 448 def breit_wigner_area(A, B, C):
 449     # Integrating the the power spectral density per unit time (power) over the
 450     # frequency axis [0, fN] returns the total signal power per unit time
 451     #  int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
 452     #                      = <V(t)**2>, the variance for our AC signal.
 453     # The variance from our fitted Breit-Wigner is the area under the Breit-Wigner
 454     #  <V(t)**2> = (pi*C) / (2*B*A**2)
 455     return (_numpy.pi * C) / (2 * B * A**2)
 456
 457 def save(filename, group='/', raw=None, config=None,
 458          deflection_channel_config=None, processed=None):
 459     specs = [
 460         _SaveSpec(item=config, relpath='config/vibration',
 461                   config=_VibrationConfig),
 462         _SaveSpec(item=deflection_channel_config,
 463                   relpath='config/deflection',
 464                   config=_InputChannelConfig),
 465         _SaveSpec(item=raw, relpath='raw', units='bits'),
 466         _SaveSpec(item=processed, relpath='processed', units='V^2/Hz'),
 467         ]
 468     _save(filename=filename, group=group, specs=specs)
 469
 470 def load(filename=None, group='/'):
 471     specs = [
 472         _SaveSpec(key=('config', 'vibration'), relpath='config/vibration',
 473                   config=_VibrationConfig),
 474         _SaveSpec(key=('config', 'deflection'), relpath='config/deflection',
 475                   config=_InputChannelConfig),
 476         _SaveSpec(key=('raw',), relpath='raw', array=True, units='bits'),
 477         _SaveSpec(key=('processed',), relpath='processed', units='V^2/Hz'),
 478         ]
 479     return _load(filename=filename, group=group, specs=specs)
 480
 481 def plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
 482              C=None, D=0, config=None, analyze=False):
 483     """Plot 3 subfigures displaying vibration data and analysis.
 484
 485      Time series (Vphoto vs sample index) (show any major drift events),
 486      A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
 487      FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
 488     """
 489     if not _matplotlib:
 490         raise _matplotlib_import_error
 491     figure = _matplotlib_pyplot.figure()
 492
 493     if power is None:
 494         assert deflection != None, (
 495             'must set at least one of `deflection` or `power`.')
 496         time_axes = figure.add_subplot(2, 1, 1)
 497         hist_axes = figure.add_subplot(2, 1, 2)
 498         freq_axes = None
 499     elif deflection is None:
 500         time_axes = hist_axes = None
 501         freq_axes = figure.add_subplot(1, 1, 1)
 502     else:
 503         time_axes = figure.add_subplot(3, 1, 1)
 504         hist_axes = figure.add_subplot(3, 1, 2)
 505         freq_axes = figure.add_subplot(3, 1, 3)
 506
 507     timestamp = _time.strftime('%H%M%S')
 508
 509     if deflection is not None:
 510         time_axes.plot(deflection, 'r.')
 511         time_axes.autoscale(tight=True)
 512         time_axes.set_title('free oscillation')
 513
 514         # plot histogram distribution and gaussian fit
 515         hist_axes.hold(True)
 516         n,bins,patches = hist_axes.hist(
 517             deflection, bins=30, normed=True, align='mid')
 518         gauss = _numpy.zeros((len(bins),), dtype=_numpy.float)
 519         mean = deflection.mean()
 520         std = deflection.std()
 521         pi = _numpy.pi
 522         exp = _numpy.exp
 523         gauss = exp(-0.5*((bins-mean)/std)**2) / (_numpy.sqrt(2*pi) * std)
 524         # Matplotlib's normed histogram uses bin heights of n/(len(x)*dbin)
 525         hist_axes.plot(bins, gauss, 'r-')
 526         hist_axes.autoscale(tight=True)
 527     if power is not None:
 528         freq_axes.hold(True)
 529         freq_axes.set_yscale('log')
 530         freq_axes.plot(freq_axis, power, 'r.-')
 531         freq_axes.autoscale(tight=True)
 532         xmin,xmax = freq_axes.get_xbound()
 533         ymin,ymax = freq_axes.get_ybound()
 534
 535         # highlight the region we're fitting
 536         freq_axes.axvspan(
 537             config['minimum-fit-frequency'],
 538             config['maximum-fit-frequency'],
 539             facecolor='g', alpha=0.1, zorder=-2)
 540
 541         if A is not None:
 542             fitdata = breit_wigner(freq_axis, A, B, C, D)
 543             freq_axes.plot(freq_axis, fitdata, 'b-')
 544             noisefloor = D + 0*freq_axis;
 545             freq_axes.plot(freq_axis, noisefloor)
 546
 547             if B**2 < 2*A**2:
 548                 res_freq = _numpy.sqrt(A**2 - B**2/2)
 549                 freq_axes.axvline(res_freq, color='b', zorder=-1)
 550
 551         freq_axes.set_title('power spectral density %s' % timestamp)
 552         freq_axes.axis([xmin,xmax,ymin,ymax])
 553         freq_axes.set_xlabel('frequency (Hz)')
 554         freq_axes.set_ylabel('PSD')
 555     if hasattr(figure, 'show'):
 556         figure.show()
 557     return figure
 558 _plot = plot  # alternative name for use inside analyze()