calibcant/vibration_analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2013 W. Trevor King <wking@tremily.us>
   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     if config['smooth-window'] and config['smooth-length']:
 310         n = config['smooth-length']
 311         smooth_window = config['smooth-window'](length=n)
 312         smooth_window /= smooth_window.sum()
 313     else:
 314         smooth_window = None
 315     A,B,C,D = fit_psd(
 316         freq_axis, power,
 317         min_frequency=config['minimum-fit-frequency'],
 318         max_frequency=config['maximum-fit-frequency'],
 319         smooth_window=smooth_window,
 320         offset=config['model'] == _OffsetBreitWigner)
 321
 322     _LOG.debug('fit PSD(f) = C / ((A**2-f**2)**2 + (f*B)**2) with '
 323                'A = %g, B = %g, C = %g, D = %g' % (A, B, C, D))
 324
 325     if plot or _package_config['matplotlib']:
 326         _plot(deflection, freq_axis, power, A, B, C, D, config=config)
 327
 328     # Our A is in uV**2, so convert back to Volts**2
 329     fitted_variance = breit_wigner_area(A,B,C) * 1e-12
 330
 331     _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
 332
 333     return min(fitted_variance, naive_variance)
 334
 335 def breit_wigner(f, A, B, C, D=0):
 336     """Breit-Wigner (sortof).
 337
 338     Inputs
 339       f  Frequency
 340       A  Resonant frequency
 341       B  Quality Q = A/B
 342       C  Scaling factor
 343       D  Optional white-noise offset
 344
 345     All parameters must be postive.
 346     """
 347     return abs(C) / ((A**2-f**2)**2 + (B*f)**2) + abs(D)
 348
 349 def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
 350             smooth_window=_FFT_tools.window_hann(10), offset=False):
 351     """Fit the FFTed vibration data to a Breit-Wigner.
 352
 353     Inputs
 354       freq_axis      array of frequencies in Hz
 355       psd_data       array of PSD amplitudes for the frequencies in freq_axis
 356       min_frequency  lower bound of Breit-Wigner fitting region
 357       max_frequency  upper bound of Breit-Wigner fitting region
 358       smooth_window  array for smoothing the PSD before guessing params
 359       offset         add a white-noise offset to the Breit-Wigner model
 360     Output
 361       Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
 362     """
 363     # cut out the relevent frequency range
 364     _LOG.debug('cut the frequency range %g to %g Hz'
 365                % (min_frequency, max_frequency))
 366     imin = 0
 367     while freq_axis[imin] < min_frequency : imin += 1
 368     imax = imin
 369     while freq_axis[imax] < max_frequency : imax += 1
 370     assert imax >= imin + 10 , 'less than 10 points in freq range (%g,%g)' % (
 371         min_frequency, max_frequency)
 372
 373     # generate guesses for Breit-Wigner parameters A, B, C, and D
 374     if hasattr(smooth_window, '__len__'):
 375         smooth_psd_data = _numpy.convolve(psd_data, smooth_window, mode='same')
 376     else:
 377         smooth_psd_data = psd_data
 378     max_psd_index = _numpy.argmax(smooth_psd_data[imin:imax]) + imin
 379     max_psd = psd_data[max_psd_index]
 380     res_freq = freq_axis[max_psd_index]
 381
 382     # Breit-Wigner L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
 383     # is expected power spectrum for
 384     # x'' + B x' + A^2 x'' = F_external(t)/m
 385     # (A = omega_0)
 386     # C = (2 k_B T B) / (pi m)
 387     #
 388     # A = resonant frequency
 389     # peak at  x_res = sqrt(A^2 - B^2/2)  (by differentiating)
 390     #  which could be complex if there isn't a peak (overdamped)
 391     # peak height    = C / (3 x_res^4 + A^4)
 392     #
 393     # Q = A/B
 394     #
 395     # Guessing Q = 1 is pretty safe.
 396
 397     Q_guess = 1
 398
 399     # so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
 400     #    B = x_res / sqrt(Q^2-1/2)
 401     B_guess = res_freq / _numpy.sqrt(Q_guess**2-0.5)
 402     A_guess = Q_guess*B_guess
 403     C_guess = max_psd * (-res_freq**4 + A_guess**4)
 404     if offset:
 405         D_guess = psd_data[-1]
 406         C_guess -= D_guess
 407     else:
 408         D_guess = 0
 409     _LOG.debug(('guessed params: resonant freq %g, max psd %g, Q %g, A %g, '
 410                 'B %g, C %g, D %g') % (
 411             res_freq, max_psd, Q_guess, A_guess, B_guess, C_guess, D_guess))
 412     # Half width w on lower side when L(a-w) = L(a)/2
 413     #  (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
 414     # Let W=(a-w)**2, A=a**2, and B=b**2
 415     #  (A - W)**2 + BW = 2*AB
 416     #  W**2 - 2AW + A**2 + BW = 2AB
 417     #  W**2 + (B-2A)W + (A**2-2AB) = 0
 418     #  W = (2A-B)/2 * [1 +/- sqrt(1 - 4(A**2-2AB)/(B-2A)**2]
 419     #    = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 420     #  (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
 421     #  so w is a disaster ;)
 422     # For some values of A and B (non-underdamped), W is imaginary or negative.
 423     #
 424     # The mass m is given by m = k_B T / (pi A)
 425     # The spring constant k is given by k = m (omega_0)**2
 426     # The quality factor Q is given by Q = omega_0 m / gamma
 427
 428     # Fitting the PSD of V = photoSensitivity*x just rescales the parameters
 429
 430     # fit Breit-Wigner using scipy.optimize.leastsq
 431     def residual(p, y, x):
 432         return breit_wigner(x, *p) - y
 433     if offset:
 434         guess = _numpy.array((A_guess, B_guess, C_guess, D_guess))
 435     else:
 436         guess = _numpy.array((A_guess, B_guess, C_guess))
 437
 438     p,cov,info,mesg,ier = _leastsq(
 439         residual, guess,
 440         args=(psd_data[imin:imax], freq_axis[imin:imax]),
 441         full_output=True, maxfev=10000)
 442     _LOG.debug('fitted params: %s' % p)
 443     _LOG.debug('covariance matrix: %s' % cov)
 444     #_LOG.debug('info: %s' % info)
 445     _LOG.debug('message: %s' % mesg)
 446     if ier == 1:
 447         _LOG.debug('solution converged')
 448     else:
 449         _LOG.debug('solution did not converge')
 450     if offset:
 451         A,B,C,D = p
 452     else:
 453         A,B,C = p
 454         D = 0
 455     A=abs(A) # A and B only show up as squares in f(x)
 456     B=abs(B) # so ensure we get positive values.
 457     C=abs(C) # Only abs(C) is used in breit_wigner().
 458     return (A, B, C, D)
 459
 460 def breit_wigner_area(A, B, C):
 461     # Integrating the the power spectral density per unit time (power) over the
 462     # frequency axis [0, fN] returns the total signal power per unit time
 463     #  int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
 464     #                      = <V(t)**2>, the variance for our AC signal.
 465     # The variance from our fitted Breit-Wigner is the area under the Breit-Wigner
 466     #  <V(t)**2> = (pi*C) / (2*B*A**2)
 467     return (_numpy.pi * C) / (2 * B * A**2)
 468
 469 def breit_wigner_resonant_frequency(A, B):
 470     if (B**2 >= 2*A**2):
 471         return 0  # over- or critically-damped
 472     return _numpy.sqrt(A**2 - B**2/2)
 473
 474 def save(filename=None, group='/', raw=None, config=None,
 475          deflection_channel_config=None, processed=None):
 476     specs = [
 477         _SaveSpec(item=config, relpath='config/vibration',
 478                   config=_VibrationConfig),
 479         _SaveSpec(item=deflection_channel_config,
 480                   relpath='config/deflection',
 481                   config=_InputChannelConfig),
 482         _SaveSpec(item=raw, relpath='raw', units='bits'),
 483         _SaveSpec(item=processed, relpath='processed', units='V^2/Hz'),
 484         ]
 485     _save(filename=filename, group=group, specs=specs)
 486
 487 def load(filename=None, group='/'):
 488     specs = [
 489         _SaveSpec(key=('config', 'vibration'), relpath='config/vibration',
 490                   config=_VibrationConfig),
 491         _SaveSpec(key=('config', 'deflection'), relpath='config/deflection',
 492                   config=_InputChannelConfig),
 493         _SaveSpec(key=('raw',), relpath='raw', array=True, units='bits'),
 494         _SaveSpec(key=('processed',), relpath='processed', units='V^2/Hz'),
 495         ]
 496     return _load(filename=filename, group=group, specs=specs)
 497
 498 def plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
 499          C=None, D=0, config=None, analyze=False):
 500     """Plot 3 subfigures displaying vibration data and analysis.
 501
 502      Time series (Vphoto vs sample index) (show any major drift events),
 503      A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
 504      FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
 505     """
 506     if not _matplotlib:
 507         raise _matplotlib_import_error
 508     figure = _matplotlib_pyplot.figure()
 509
 510     if power is None:
 511         assert deflection != None, (
 512             'must set at least one of `deflection` or `power`.')
 513         time_axes = figure.add_subplot(2, 1, 1)
 514         hist_axes = figure.add_subplot(2, 1, 2)
 515         freq_axes = None
 516     elif deflection is None:
 517         time_axes = hist_axes = None
 518         freq_axes = figure.add_subplot(1, 1, 1)
 519     else:
 520         time_axes = figure.add_subplot(3, 1, 1)
 521         hist_axes = figure.add_subplot(3, 1, 2)
 522         freq_axes = figure.add_subplot(3, 1, 3)
 523
 524     timestamp = _time.strftime('%H%M%S')
 525
 526     if deflection is not None:
 527         time_axes.plot(deflection, 'r.')
 528         time_axes.autoscale(tight=True)
 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 = exp(-0.5*((bins-mean)/std)**2) / (_numpy.sqrt(2*pi) * std)
 541         # Matplotlib's normed histogram uses bin heights of n/(len(x)*dbin)
 542         hist_axes.plot(bins, gauss, 'r-')
 543         hist_axes.autoscale(tight=True)
 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         freq_axes.autoscale(tight=True)
 549         xmin,xmax = freq_axes.get_xbound()
 550         ymin,ymax = freq_axes.get_ybound()
 551
 552         # highlight the region we're fitting
 553         if config:
 554             freq_axes.axvspan(
 555                 config['minimum-fit-frequency'],
 556                 config['maximum-fit-frequency'],
 557                 facecolor='g', alpha=0.1, zorder=-2)
 558
 559         if A is not None:
 560             fitdata = breit_wigner(freq_axis, A, B, C, D)
 561             freq_axes.plot(freq_axis, fitdata, 'b-')
 562             noisefloor = D + 0*freq_axis;
 563             freq_axes.plot(freq_axis, noisefloor)
 564             res_freq = breit_wigner_resonant_frequency(A=A, B=B)
 565             if res_freq > 0:
 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()
 574     return figure
 575 _plot = plot  # alternative name for use inside analyze()