calibcant/vib_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 `vib_analyze()` from the other `vib_*()`
  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 >>> vibration_config = VibrationConfig()
  42 >>> vibration_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/vibration_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 = vibration_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_vibration = vib_analyze_naive(deflection)
 135 >>> naive_vibration  # doctest: +SKIP
 136 136871517.43486854
 137 >>> abs(naive_vibration / 136.9e6 - 1) < 0.1
 138 True
 139
 140 >>> processed_vibration = vib_analyze(
 141 ...     deflection_bits, vibration_config,
 142 ...     piezo_config.select_config('inputs', 'deflection'))
 143 >>> processed_vibration  # doctest: +SKIP
 144 136457906.25574699
 145
 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)
 152
 153 >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 154 /
 155   /vibration
 156     /vibration/config
 157       /vibration/config/deflection
 158         <HDF5 dataset "channel": shape (), type "<i4">
 159           0
 160         <HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
 161           [0 1]
 162         <HDF5 dataset "conversion-origin": shape (), type "<i4">
 163           0
 164         <HDF5 dataset "device": shape (), type "|S12">
 165           /dev/comedi0
 166         <HDF5 dataset "inverse-conversion-coefficients": shape (2,), type "<i4">
 167           [0 1]
 168         <HDF5 dataset "inverse-conversion-origin": shape (), type "<i4">
 169           0
 170         <HDF5 dataset "maxdata": shape (), type "<i4">
 171           100
 172         <HDF5 dataset "name": shape (), type "|S10">
 173           deflection
 174         <HDF5 dataset "range": shape (), type "<i4">
 175           1
 176         <HDF5 dataset "subdevice": shape (), type "<i4">
 177           -1
 178       /vibration/config/vibration
 179         <HDF5 dataset "chunk-size": shape (), type "<i4">
 180           2048
 181         <HDF5 dataset "frequency": shape (), type "<f8">
 182           50000.0
 183         <HDF5 dataset "maximum-fit-frequency": shape (), type "<f8">
 184           25000.0
 185         <HDF5 dataset "minimum-fit-frequency": shape (), type "<f8">
 186           500.0
 187         <HDF5 dataset "model": shape (), type "|S12">
 188           Breit-Wigner
 189         <HDF5 dataset "overlap": shape (), type "|b1">
 190           False
 191         <HDF5 dataset "sample-time": shape (), type "<i4">
 192           1
 193         <HDF5 dataset "window": shape (), type "|S4">
 194           Hann
 195     <HDF5 dataset "processed": shape (), type "<f8">
 196       ...
 197     /vibration/raw
 198       <HDF5 dataset "deflection": shape (32768,), type "<f8">
 199         [...]
 200
 201 >>> (raw_vibration,vibration_config,deflection_channel_config,
 202 ...  processed_vibration) = vib_load(
 203 ...     filename=filename, group='/vibration/')
 204
 205 >>> processed_vibration  # doctest: +SKIP
 206 136457906.25574699
 207 >>> abs(processed_vibration / 136.5e6 - 1) < 0.1
 208 True
 209
 210 >>> os.remove(filename)
 211 """
 212
 213 import os as _os
 214 import time as _time
 215
 216 import h5py as _h5py
 217 import numpy as _numpy
 218 from scipy.optimize import leastsq as _leastsq
 219
 220 try:
 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:
 225     _matplotlib = None
 226     _matplotlib_import_error = e
 227
 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
 233
 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
 240
 241
 242 def vib_analyze_naive(deflection):
 243     """Calculate the deflection variance in Volts**2.
 244
 245     This method is simple and easy to understand, but it highly
 246     succeptible to noise, drift, etc.
 247
 248     Inputs
 249       deflection : numpy array with deflection timeseries in Volts.
 250     """
 251     std = deflection.std()
 252     var = std**2
 253     _LOG.debug('naive deflection variance: %g V**2' % var)
 254     return var
 255
 256 def vib_analyze(deflection, vibration_config, deflection_channel_config,
 257                 plot=False):
 258     """Calculate the deflection variance in Volts**2.
 259
 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.
 267
 268     Inputs:
 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.
 274
 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.
 279     """
 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)
 285
 286     naive_variance = vib_analyze_naive(deflection_v)
 287     if vibration_config['model'] == _Variance:
 288         return naive_variance
 289
 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'])
 296
 297     A,B,C,D = fit_psd(
 298         freq_axis, power,
 299         min_frequency=vibration_config['minimum-fit-frequency'],
 300         max_frequency=vibration_config['maximum-fit-frequency'],
 301         offset=vibration_config['model'] == _OffsetBreitWigner)
 302
 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))
 305
 306     if plot or _package_config['matplotlib']:
 307         vib_plot(deflection, freq_axis, power, A, B, C, D,
 308                  vibration_config=vibration_config)
 309
 310     # Our A is in uV**2, so convert back to Volts**2
 311     fitted_variance = breit_wigner_area(A,B,C) * 1e-12
 312
 313     _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
 314
 315     if _package_config['matplotlib']:
 316         vib_plot(deflection, freq_axis, power, A, B, C, D,
 317                  vibration_config=vibration_config)
 318
 319     return min(fitted_variance, naive_variance)
 320
 321 def breit_wigner(f, A, B, C, D=0):
 322     """Breit-Wigner (sortof).
 323
 324     Inputs
 325       f  Frequency
 326       A  Resonant frequency
 327       B  Quality Q = A/B
 328       C  Scaling factor
 329       D  Optional white-noise offset
 330
 331     All parameters must be postive.
 332     """
 333     return abs(C) / ((A**2-f**2)**2 + (B*f)**2) + abs(D)
 334
 335 def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
 336             offset=False):
 337     """Fit the FFTed vibration data to a Breit-Wigner.
 338
 339     Inputs
 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
 345     Output
 346       Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
 347     """
 348     # cut out the relevent frequency range
 349     _LOG.debug('cut the frequency range %g to %g Hz'
 350                % (min_frequency, max_frequency))
 351     imin = 0
 352     while freq_axis[imin] < min_frequency : imin += 1
 353     imax = imin
 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)
 357
 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]
 362
 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
 366     # (A = omega_0)
 367     # C = (2 k_B T B) / (pi m)
 368     #
 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)
 373     #
 374     # Q = A/B
 375     #
 376     # Guessing Q = 1 is pretty safe.
 377
 378     Q_guess = 1
 379
 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)
 385     if offset:
 386         D_guess = psd_data[max_psd_index]
 387         C_guess -= D_guess
 388     else:
 389         D_guess = 0
 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.
 404     #
 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
 408
 409     # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
 410
 411     # fit Breit-Wigner using scipy.optimize.leastsq
 412     def residual(p, y, x):
 413         return breit_wigner(x, *p) - y
 414     if offset:
 415         guess = _numpy.array((A_guess, B_guess, C_guess, D_guess))
 416     else:
 417         guess = _numpy.array((A_guess, B_guess, C_guess))
 418
 419     p,cov,info,mesg,ier = _leastsq(
 420         residual, guess,
 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)
 427     if ier == 1:
 428         _LOG.debug('solution converged')
 429     else:
 430         _LOG.debug('solution did not converge')
 431     if offset:
 432         A,B,C,D = p
 433     else:
 434         A,B,C = p
 435         D = 0
 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().
 439     return (A, B, C, D)
 440
 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)
 449
 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:
 455             try:
 456                 del cwg['raw/deflection']
 457             except KeyError:
 458                 pass
 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')]:
 464             if config is None:
 465                 continue
 466             config_cwg = _h5_create_group(cwg, key)
 467             storage.save(config=config, 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 [(_VibrationConfig, 'config/vibration'),
 485                            (_ChannelConfig, 'config/deflection')]:
 486             config = Config(storage=_HDF5_Storage(
 487                     filename=filename, group=group+key))
 488             configs.append(config)
 489         try:
 490             processed_vibration = float(f[group+'processed'][...])
 491         except KeyError:
 492             pass
 493     ret = [raw_vibration]
 494     ret.extend(configs)
 495     ret.append(processed_vibration)
 496     for config in configs:
 497         config.load()
 498     return tuple(ret)
 499
 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.
 503
 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).
 507     """
 508     if not _matplotlib:
 509         raise _matplotlib_import_error
 510     figure = _matplotlib_pyplot.figure()
 511
 512     if power is None:
 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)
 517         freq_axes = None
 518     elif deflection is None:
 519         time_axes = hist_axes = None
 520         freq_axes = figure.add_subplot(1, 1, 1)
 521     else:
 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)
 525
 526     timestamp = _time.strftime('%H%M%S')
 527
 528     if deflection is not None:
 529         time_axes.plot(deflection, 'r.')
 530         time_axes.set_title('free oscillation')
 531
 532         # plot histogram distribution and gaussian fit
 533         hist_axes.hold(True)
 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()
 539         pi = _numpy.pi
 540         exp = _numpy.exp
 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:
 546         freq_axes.hold(True)
 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()
 551
 552         # highlight the region we're fitting
 553         freq_axes.axvspan(
 554             vibration_config['minimum-fit-frequency'],
 555             vibration_config['maximum-fit-frequency'],
 556             facecolor='g', alpha=0.1, zorder=-2)
 557
 558         if A is not None:
 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)
 563
 564             if B**2 < 2*A**2:
 565                 res_freq = _numpy.sqrt(A**2 - B**2/2)
 566                 freq_axes.axvline(res_freq, color='b', zorder=-1)
 567
 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'):
 573         figure.show()