calibcant/vib_analyze.py

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