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     C=abs(C) # Only abs(C) is used in breit_wigner().
 441     return (A, B, C, D)
 442
 443 def breit_wigner_area(A, B, C):
 444     # Integrating the the power spectral density per unit time (power) over the
 445     # frequency axis [0, fN] returns the total signal power per unit time
 446     #  int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
 447     #                      = <V(t)**2>, the variance for our AC signal.
 448     # The variance from our fitted Breit-Wigner is the area under the Breit-Wigner
 449     #  <V(t)**2> = (pi*C) / (2*B*A**2)
 450     return (_numpy.pi * C) / (2 * B * A**2)
 451
 452 def vib_save(filename, group='/', raw_vibration=None, vibration_config=None,
 453              deflection_channel_config=None, processed_vibration=None):
 454     with _h5py.File(filename, 'a') as f:
 455         cwg = _h5_create_group(f, group)
 456         if raw_vibration is not None:
 457             try:
 458                 del cwg['raw/deflection']
 459             except KeyError:
 460                 pass
 461             cwg['raw/deflection'] = raw_vibration
 462         for config,key in [(vibration_config, 'config/vibration'),
 463                            (deflection_channel_config,
 464                             'config/deflection')]:
 465             if config is None:
 466                 continue
 467             config_cwg = _h5_create_group(cwg, key)
 468             config.save(group=config_cwg)
 469         if processed_vibration is not None:
 470             try:
 471                 del cwg['processed']
 472             except KeyError:
 473                 pass
 474             cwg['processed'] = processed_vibration
 475
 476 def vib_load(filename, group='/'):
 477     assert group.endswith('/')
 478     raw_vibration = processed_vibration = None
 479     configs = []
 480     with _h5py.File(filename, 'a') as f:
 481         try:
 482             raw_vibration = f[group+'raw/deflection'][...]
 483         except KeyError:
 484             pass
 485         for Config,key in [(_HDF5_VibrationConfig, 'config/vibration'),
 486                            (_HDF5_ChannelConfig, 'config/deflection')]:
 487             config = Config(filename=filename, group=group+key)
 488             configs.append(config)
 489         try:
 490             processed_vibration = 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
 573     figure.show()