# License along with calibcant. If not, see
# <http://www.gnu.org/licenses/>.
-"""
-Separate the more general vib_analyze() from the other vib_*()
-functions in calibcant. Also provide a command line interface
-for analyzing data acquired through other workflows.
+"""Thermal vibration analysis.
+
+Separate the more general `vib_analyze()` from the other `vib_*()`
+functions in calibcant.
The relevent physical quantities are :
- Vphoto The photodiode vertical deflection voltage (what we measure)
+ Vphoto The photodiode vertical deflection voltage (what we measure)
+
+>>> import os
+>>> from pprint import pprint
+>>> import random
+>>> import tempfile
+>>> import numpy
+>>> from .config import HDF5_VibrationConfig
+>>> from pypiezo.config import pprint_HDF5, HDF5_ChannelConfig
+>>> from pypiezo.base import convert_volts_to_bits
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> vibration_config = HDF5_VibrationConfig(
+... filename=filename, group='/vibration/config/vibration')
+>>> vibration_config['frequency'] = 50e3
+>>> vibration_config.save()
+>>> deflection_channel_config = HDF5_ChannelConfig(filename=None)
+>>> deflection_channel_config['maxdata'] = 200
+>>> deflection_channel_config['conversion-coefficients'] = (0,1)
+>>> deflection_channel_config['conversion-origin'] = 0
+>>> deflection_channel_config['inverse-conversion-coefficients'] = (0,1)
+>>> deflection_channel_config['inverse-conversion-origin'] = 0
+
+We'll be generating a test vibration time series with the following
+parameters. Make sure these are all floats to avoid accidental
+integer division in later steps.
+
+>>> m = 5e-11 # kg
+>>> gamma = 1.6e-6 # N*s/m
+>>> k = 0.05 # N/m
+>>> T = 1/vibration_config['frequency']
+>>> T # doctest: +ELLIPSIS
+2.00...e-05
+>>> N = int(2**15) # count
+>>> F_sigma = 1e3 # N
+
+where `T` is the sampling period, `N` is the number of samples, and
+`F_sigma` is the standard deviation of the white-noise external force.
+Note that the resonant frequency is less than the Nyquist frequency so
+we don't have to worry too much about aliasing.
+
+>>> w0 = numpy.sqrt(k/m)
+>>> f0 = w0/(2*numpy.pi)
+>>> f0 # doctest: +ELLIPSIS
+5032.9...
+>>> f_nyquist = vibration_config['frequency']/2
+>>> f_nyquist # doctest: +ELLIPSIS
+25000.0...
+
+The damping ratio is
+
+>>> damping = gamma / (2*m*w0)
+>>> damping # doctest: +ELLIPSIS
+0.505...
+
+The quality factor is
+
+>>> Q = m*w0 / gamma
+>>> Q # doctest: +ELLIPSIS
+0.988...
+>>> (1 / (2*damping)) / Q # doctest: +ELLIPSIS
+1.000...
+
+We expect the white-noise power spectral density (PSD) to be a flat
+line at
+
+>>> F0 = F_sigma**2 * 2 * T
+
+because the integral from `0` `1/2T` should be `F_sigma**2`.
+
+The expected time series PSD parameters are
+
+>>> A = f0
+>>> B = gamma/(m*2*numpy.pi)
+>>> C = F0/(m**2*(2*numpy.pi)**4)
+
+Simulate a time series with the proper PSD using center-differencing.
+
+ m\ddot{x} + \gamma \dot{x} + kx = F
+
+ m \frac{x_{i+1} - 2x_i + x_{i-1}}{T**2}
+ + \gamma \frac{x_{i+1}-x_{i-1}}{T}
+ + kx_i = F_i
+
+ a x_{i+1} + b x_{i} + c x_{i-1} = F_i
+
+where `T` is the sampling period, `i=t/T` is the measurement index,
+`a=m/T**2+gamma/2T`, `b=-2m/T**2+k`, and `c=m/T**2-gamma/2T`.
+Rearranging and shifting to `j=i-1`
+
+ x_j = \frac{F_{i-1} - bx_{i-1} - cx_{i-2}}{a}
+
+>>> a = m/T**2 + gamma/(2*T)
+>>> b = -2*m/T**2 + k
+>>> c = m/T**2 - gamma/(2*T)
+>>> x = numpy.zeros((N+2,), dtype=numpy.float) # two extra initial points
+>>> F = numpy.zeros((N,), dtype=numpy.float)
+>>> for i in range(2, x.size):
+... Fp = random.gauss(mu=0, sigma=F_sigma)
+... F[i-2] = Fp
+... xp = x[i-1]
+... xpp = x[i-2]
+... x[i] = (Fp - b*xp - c*xpp)/a
+>>> x = x[2:] # drop the initial points
+
+Convert the simulated data to bits.
+
+>>> deflection = x
+>>> deflection_bits = convert_volts_to_bits(deflection_channel_config, x)
+
+Analyze the simulated data.
+
+>>> naive_vibration = vib_analyze_naive(deflection)
+>>> naive_vibration # doctest: +SKIP
+136871517.43486854
+>>> abs(naive_vibration / 136.9e6 - 1) < 0.1
+True
+
+>>> processed_vibration = vib_analyze(
+... deflection_bits, vibration_config, deflection_channel_config)
+>>> processed_vibration # doctest: +SKIP
+136457906.25574699
+
+>>> vib_plot(deflection=deflection_bits, vibration_config=vibration_config)
+>>> vib_save(filename=filename, group='/vibration/',
+... raw_vibration=deflection_bits, vibration_config=vibration_config,
+... deflection_channel_config=deflection_channel_config,
+... processed_vibration=processed_vibration)
+
+>>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+/
+ /vibration
+ /vibration/config
+ /vibration/config/deflection
+ <HDF5 dataset "channel": shape (), type "|S1">
+ 0
+ <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
+ 0, 1
+ <HDF5 dataset "conversion-origin": shape (), type "|S1">
+ 0
+ <HDF5 dataset "device": shape (), type "|S12">
+ /dev/comedi0
+ <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S4">
+ 0, 1
+ <HDF5 dataset "inverse-conversion-origin": shape (), type "|S1">
+ 0
+ <HDF5 dataset "maxdata": shape (), type "|S3">
+ 200
+ <HDF5 dataset "range": shape (), type "|S1">
+ 1
+ <HDF5 dataset "subdevice": shape (), type "|S2">
+ -1
+ /vibration/config/vibration
+ <HDF5 dataset "chunk-size": shape (), type "|S4">
+ 2048
+ <HDF5 dataset "frequency": shape (), type "|S7">
+ 50000.0
+ <HDF5 dataset "maximum-fit-frequency": shape (), type "|S7">
+ 25000.0
+ <HDF5 dataset "minimum-fit-frequency": shape (), type "|S5">
+ 500.0
+ <HDF5 dataset "model": shape (), type "|S12">
+ Breit-Wigner
+ <HDF5 dataset "overlap": shape (), type "|S2">
+ no
+ <HDF5 dataset "sample-time": shape (), type "|S1">
+ 1
+ <HDF5 dataset "window": shape (), type "|S4">
+ Hann
+ <HDF5 dataset "processed": shape (), type "<f8">
+ ...
+ /vibration/raw
+ <HDF5 dataset "deflection": shape (32768,), type "<f8">
+ [...]
+
+>>> (raw_vibration,vibration_config,deflection_channel_config,
+... processed_vibration) = vib_load(
+... filename=filename, group='/vibration/')
+
+>>> processed_vibration # doctest: +SKIP
+136457906.25574699
+>>> abs(processed_vibration / 136.5e6 - 1) < 0.1
+True
+
+>>> os.remove(filename)
"""
-import os, time, numpy
-import GnuplotBiDir # can be used for fitting lorentzian
-import scipy.optimize # can be used for fitting lorentzian
+import os as _os
+import time as _time
-import data_logger
-import FFT_tools
-from splittable_kwargs import splittableKwargsFunction, \
- make_splittable_kwargs_function
+import h5py as _h5py
+import numpy as _numpy
+from scipy.optimize import leastsq as _leastsq
-from . import common
-from . import config
+try:
+ import matplotlib as _matplotlib
+ import matplotlib.pyplot as _matplotlib_pyplot
+ import time as _time # for timestamping lines on plots
+except (ImportError, RuntimeError), e:
+ _matplotlib = None
+ _matplotlib_import_error = e
+import FFT_tools as _FFT_tools
+from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
+from pypiezo.config import HDF5_ChannelConfig as _HDF5_ChannelConfig
+from pypiezo.config import h5_create_group as _h5_create_group
-PLOT_GUESSED_LORENTZIAN=False
+from . import LOG as _LOG
+from . import base_config as _base_config
+from .config import Variance as _Variance
+from .config import BreitWigner as _BreitWigner
+from .config import OffsetBreitWigner as _OffsetBreitWigner
+from .config import HDF5_VibrationConfig as _HDF5_VibrationConfig
-@splittableKwargsFunction()
-def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V,
- zeroVphoto_bits=config.zeroVphoto_bits) :
- """
- Calculate the variance of the raw data, and convert it to Volts**2.
- This method is simple and easy to understand, but it highly succeptible
- to noise, drift, etc.
+
+def vib_analyze_naive(deflection):
+ """Calculate the deflection variance in Volts**2.
+
+ This method is simple and easy to understand, but it highly
+ succeptible to noise, drift, etc.
- Vphoto_in2V is a function converting Vphoto values from bits to Volts.
- This function is assumed linear, see bump_analyze().
- """
- if config.TEXT_VERBOSE :
- print "Calculating the naive (raw) variance"
- Vphoto_std_bits = deflection_bits.std()
- Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zeroVphoto_bits)
- if config.TEXT_VERBOSE :
- print "Average deflection (Volts):", Vphoto_in2V(deflection_bits.mean())
- return Vphoto_std**2
-
-#@splittableKwargsFunction((fit_psd, 'freq_axis', 'psd_data'),
-# (vib_plot, 'deflection_bits', 'freq_axis', 'power',
-# 'A', 'B', 'C', 'minFreq', 'maxFreq'))
-# Some of the child functions aren't yet defined, so postpone
-# make-splittable until later in the module.
-def vib_analyze(deflection_bits, freq,
- chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V, **kwargs) :
+ Inputs
+ deflection : numpy array with deflection timeseries in Volts.
"""
- Calculate the variance in the raw data due to the cantilever's
- thermal oscillation and convert it to Volts**2.
+ std = deflection.std()
+ var = std**2
+ _LOG.debug('naive deflection variance: %g V**2' % var)
+ return var
- Improves on vib_analyze_naive() by first converting Vphoto(t) to
- frequency space, and fitting a Lorentzian in the relevant
+def vib_analyze(deflection, vibration_config, deflection_channel_config,
+ plot=False):
+ """Calculate the deflection variance in Volts**2.
+
+ Improves on `vib_analyze_naive()` by first converting `Vphoto(t)`
+ to frequency space, and fitting a Breit-Wigner in the relevant
frequency range (see cantilever_calib.pdf for derivation).
However, there may be cases where the fit is thrown off by noise
spikes in frequency space. To protect from errors, the fitted
variance is compared to the naive variance (where all noise is
included), and the minimum variance is returned.
- Vphoto_in2V is a function converting Vphoto values from bits to Volts.
- This function is assumed linear, see bump_analyze().
+ Inputs:
+ deflection Vphoto deflection input in bits.
+ vibration_config `.config._VibrationConfig` instance
+ deflection_channel_config
+ deflection `pypiezo.config._ChannelConfig` instance
+ plot boolean overriding matplotlib config setting.
+
+ The conversion to frequency space generates an average power
+ spectrum by breaking the data into windowed chunks and averaging
+ the power spectrums for the chunks together. See
+ `FFT_tools.unitary_avg_power_spectrum()` for details.
"""
- fit_psd_kwargs,vib_plot_kwargs = vib_analyze._splitargs(vib_analyze,kwargs)
- if 'minFreq' in fit_psd_kwargs:
- vib_plot_kwargs['minFreq'] = fit_psd_kwargs['minFreq']
- if 'maxFreq' in fit_psd_kwargs:
- vib_plot_kwargs['maxFreq'] = fit_psd_kwargs['maxFreq']
-
- Vphoto_t = numpy.zeros((len(deflection_bits),),
- dtype=numpy.float)
# convert the data from bits to volts
- if config.TEXT_VERBOSE :
- print "Converting %d bit values to voltages" % len(Vphoto_t)
- for i in range(len(deflection_bits)) :
- Vphoto_t[i] = Vphoto_in2V(deflection_bits[i])
- if config.TEXT_VERBOSE :
- print "Average deflection (Volts):", Vphoto_t.mean()
+ deflection_v = _convert_bits_to_volts(
+ deflection_channel_config, deflection)
+ mean = deflection_v.mean()
+ _LOG.debug('average thermal deflection (Volts): %g' % mean)
+
+ naive_variance = vib_analyze_naive(deflection_v)
+ if vibration_config['model'] == _Variance:
+ return naive_variance
# Compute the average power spectral density per unit time (in uV**2/Hz)
- if config.TEXT_VERBOSE :
- print "Compute the averaged power spectral density in uV**2/Hz"
- (freq_axis, power) = \
- FFT_tools.unitary_avg_power_spectrum((Vphoto_t - Vphoto_t.mean())*1e6,
- freq, chunk_size,overlap, window)
+ _LOG.debug('compute the averaged power spectral density in uV**2/Hz')
+ freq_axis,power = _FFT_tools.unitary_avg_power_spectrum(
+ (deflection_v - mean)*1e6, vibration_config['frequency'],
+ vibration_config['chunk-size'], vibration_config['overlap'],
+ vibration_config['window'])
- A,B,C,D = fit_psd(freq_axis, power, **kwargs)
+ A,B,C,D = fit_psd(
+ freq_axis, power,
+ min_frequency=vibration_config['minimum-fit-frequency'],
+ max_frequency=vibration_config['maximum-fit-frequency'],
+ offset=vibration_config['model'] == _OffsetBreitWigner)
- if config.TEXT_VERBOSE :
- print "Fitted f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
- print "A = %g, \t B = %g, \t C = %g, \t D = %g" % (A, B, C, D)
+ _LOG.debug('fit PSD(f) = C / ((A**2-f**2)**2 + (f*B)**2) with '
+ 'A = %g, B = %g, C = %g, D = %g' % (A, B, C, D))
- vib_plot(deflection_bits, freq_axis, power, A, B, C, D, **kwargs)
+ if plot or _base_config['matplotlib']:
+ vib_plot(deflection, freq_axis, power, A, B, C, D,
+ vibration_config=vibration_config)
# Our A is in uV**2, so convert back to Volts**2
- fitted_variance = lorentzian_area(A,B,C) * 1e-12
+ fitted_variance = breit_wigner_area(A,B,C) * 1e-12
+
+ _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
+
+ if _base_config['matplotlib']:
+ vib_plot(deflection, freq_axis, power, A, B, C, D,
+ vibration_config=vibration_config)
- naive_variance = vib_analyze_naive(deflection_bits,
- Vphoto_in2V=Vphoto_in2V)
- if config.TEXT_VERBOSE:
- print "Fitted variance: %g V**2" % fitted_variance
- print "Naive variance : %g V**2" % naive_variance
-
return min(fitted_variance, naive_variance)
-@splittableKwargsFunction()
-def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=25000) :
- """
- freq_axis : array of frequencies in Hz
- psd_data : array of PSD amplitudes for the frequencies in freq_axis
+def breit_wigner(f, A, B, C, D=0):
+ """Breit-Wigner (sortof).
- Calculate the variance in the raw data due to the cantilever's
- thermal oscillation and convert it to Volts**2.
+ Inputs
+ f Frequency
+ A Resonant frequency
+ B Quality Q = A/B
+ C Scaling factor
+ D Optional white-noise offset
- The conversion to frequency space generates an average power spectrum
- by breaking the data into windowed chunks and averaging the power
- spectrums for the chunks together.
- See FFT_tools.unitary_avg_power_spectrum().
+ All parameters must be postive.
+ """
+ return abs(C) / ((A**2-f**2)**2 + (B*f)**2) + abs(D)
+
+def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
+ offset=False):
+ """Fit the FFTed vibration data to a Breit-Wigner.
+
+ Inputs
+ freq_axis array of frequencies in Hz
+ psd_data array of PSD amplitudes for the frequencies in freq_axis
+ min_frequency lower bound of Breit-Wigner fitting region
+ max_frequency upper bound of Breit-Wigner fitting region
+ offset add a white-noise offset to the Breit-Wigner model
+ Output
+ Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
"""
# cut out the relevent frequency range
- if config.TEXT_VERBOSE :
- print "Cut the frequency range %g to %g Hz" % (minFreq, maxFreq)
+ _LOG.debug('cut the frequency range %g to %g Hz'
+ % (min_frequency, max_frequency))
imin = 0
- while freq_axis[imin] < minFreq : imin += 1
+ while freq_axis[imin] < min_frequency : imin += 1
imax = imin
- while freq_axis[imax] < maxFreq : imax += 1
- assert imax >= imin + 10 , \
- "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq)
+ while freq_axis[imax] < max_frequency : imax += 1
+ assert imax >= imin + 10 , 'less than 10 points in freq range (%g,%g)' % (
+ min_frequency, max_frequency)
- # generate guesses for Lorentzian parameters A,B,C
- max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin
+ # generate guesses for Breit-Wigner parameters A, B, C, and D
+ max_psd_index = _numpy.argmax(psd_data[imin:imax]) + imin
max_psd = psd_data[max_psd_index]
- half_max_index = imin
- while psd_data[half_max_index] < max_psd/2 :
- half_max_index += 1
res_freq = freq_axis[max_psd_index]
- half_width = freq_axis[max_psd_index] - freq_axis[half_max_index]
- # Lorentzian L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
+
+ # Breit-Wigner L(x) = C / ((A**2-x**2)**2 + (B*x)**2)
# is expected power spectrum for
# x'' + B x' + A^2 x'' = F_external(t)/m
# (A = omega_0)
# C = (2 k_B T B) / (pi m)
#
# A = resonant frequency
- # peak at x_res = sqrt(A^2 - B^2/2)
+ # peak at x_res = sqrt(A^2 - B^2/2) (by differentiating)
# which could be complex if there isn't a peak (overdamped)
# peak height = C / (3 x_res^4 + A^4)
+ #
# Q = A/B
#
- # guess Q = 1, and adjust from there
+ # Guessing Q = 1 is pretty safe.
+
Q_guess = 1
+
# so x_res^2 = B^2 Q^2 - B^2/2 = (Q^2-1/2)B^2
# B = x_res / sqrt(Q^2-1/2)
- if config.TEXT_VERBOSE :
- print "params : %g, %g" % (res_freq, max_psd)
- B_guess = res_freq / numpy.sqrt(Q_guess**2-0.5)
+ B_guess = res_freq / _numpy.sqrt(Q_guess**2-0.5)
A_guess = Q_guess*B_guess
C_guess = max_psd * (-res_freq**4 + A_guess**4)
- D_guess = 0 # allow a noise floor offset (DISABLED in fitting below)
- #
+ if offset:
+ D_guess = psd_data[max_psd_index]
+ C_guess -= D_guess
+ else:
+ D_guess = 0
+ _LOG.debug(('guessed params: resonant freq %g, max psd %g, Q %g, A %g, '
+ 'B %g, C %g, D %g') % (
+ res_freq, max_psd, Q_guess, A_guess, B_guess, C_guess, D_guess))
# Half width w on lower side when L(a-w) = L(a)/2
# (a**2 - (a-w)**2)**2 + (b*(a-w))**2 = 2*(b*a)**2
# Let W=(a-w)**2, A=a**2, and B=b**2
# The mass m is given by m = k_B T / (pi A)
# The spring constant k is given by k = m (omega_0)**2
# The quality factor Q is given by Q = omega_0 m / gamma
- if config.TEXT_VERBOSE :
- print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
- if PLOT_GUESSED_LORENTZIAN:
- vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess,
- minFreq, maxFreq, plotVerbose=True)
# Fitting the PSD of V = photoSensitivity*x just rescales the parameters
- if False: # Gnuplot fit worse than scipy.optimize.leastsq fit.
- # fit Lorentzian using Gnuplot's 'fit' command
-
- # save about-to-be-fitted data to a temporary file
- datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time())
- fd = file(datafilename, 'w')
- for i in range(imin, imax) :
- fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i]))
- fd.close()
-
- g = GnuplotBiDir.Gnuplot()
- g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)")
- g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess))
- g("fit f(x) '%s' via A,B,C" % datafilename)
- A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x)
- B=abs(float(g.getvar('B'))) # so ensure we get positive values
- C=float(g.getvar('C'))
-
- os.remove(datafilename)
+ # fit Breit-Wigner using scipy.optimize.leastsq
+ def residual(p, y, x):
+ return breit_wigner(x, *p) - y
+ if offset:
+ guess = _numpy.array((A_guess, B_guess, C_guess, D_guess))
+ else:
+ guess = _numpy.array((A_guess, B_guess, C_guess))
+
+ p,cov,info,mesg,ier = _leastsq(
+ residual, guess,
+ args=(psd_data[imin:imax], freq_axis[imin:imax]),
+ full_output=True, maxfev=10000)
+ _LOG.debug('fitted params: %s' % p)
+ _LOG.debug('covariance matrix: %s' % cov)
+ #_LOG.debug('info: %s' % info)
+ _LOG.debug('message: %s' % mesg)
+ if ier == 1:
+ _LOG.debug('solution converged')
+ else:
+ _LOG.debug('solution did not converge')
+ if offset:
+ A,B,C,D = p
else:
- # fit Lorenzian using scipy.optimize.leastsq
- def residual(p, y, x):
- A = p[0]
- B = p[1]
- C = p[2]
- Y = C / ((A**2-x**2)**2 + (B*x)**2)
- return Y-y
- leastsq = scipy.optimize.leastsq
- p,cov,info,mesg,ier = leastsq(residual, (A_guess, B_guess, C_guess),
- args=(psd_data[imin:imax],
- freq_axis[imin:imax]),
- full_output=True, maxfev=10000)
- if config.TEXT_VERBOSE:
- print "Fitted params:",p
- print "Covariance mx:",cov
- #print "Info:", info # too much information :p
- print "mesg:", mesg
- if ier == 1 :
- print "Solution converged"
- else :
- print "Solution did not converge"
A,B,C = p
- D = 0 # do not fit the noise floor (fit doesn't look very convincing)
- A=abs(A) # A and B only show up as squares in f(x)
- B=abs(B) # so ensure we get positive values
+ D = 0
+ A=abs(A) # A and B only show up as squares in f(x)
+ B=abs(B) # so ensure we get positive values
return (A, B, C, D)
-def lorentzian_area(A, B, C) :
+def breit_wigner_area(A, B, C):
# Integrating the the power spectral density per unit time (power) over the
# frequency axis [0, fN] returns the total signal power per unit time
# int_0^fN power(f)df = 1/T int_0^T |x(t)**2| dt
# = <V(t)**2>, the variance for our AC signal.
- # The variance from our fitted Lorentzian is the area under the Lorentzian
+ # The variance from our fitted Breit-Wigner is the area under the Breit-Wigner
# <V(t)**2> = (pi*C) / (2*B*A**2)
- return (numpy.pi * C) / (2 * B * A**2)
+ return (_numpy.pi * C) / (2 * B * A**2)
+
+def vib_save(filename, group='/', raw_vibration=None, vibration_config=None,
+ deflection_channel_config=None, processed_vibration=None):
+ with _h5py.File(filename, 'a') as f:
+ cwg = _h5_create_group(f, group)
+ if raw_vibration is not None:
+ try:
+ del cwg['raw/deflection']
+ except KeyError:
+ pass
+ cwg['raw/deflection'] = raw_vibration
+ for config,key in [(vibration_config, 'config/vibration'),
+ (deflection_channel_config,
+ 'config/deflection')]:
+ if config is None:
+ continue
+ config_cwg = _h5_create_group(cwg, key)
+ config.save(group=config_cwg)
+ if processed_vibration is not None:
+ try:
+ del cwg['processed']
+ except KeyError:
+ pass
+ cwg['processed'] = processed_vibration
+
+def vib_load(filename, group='/'):
+ assert group.endswith('/')
+ raw_vibration = processed_vibration = None
+ configs = []
+ with _h5py.File(filename, 'a') as f:
+ try:
+ raw_vibration = f[group+'raw/deflection'][...]
+ except KeyError:
+ pass
+ for Config,key in [(_HDF5_VibrationConfig, 'config/vibration'),
+ (_HDF5_ChannelConfig, 'config/deflection')]:
+ config = Config(filename=filename, group=group+key)
+ configs.append(config)
+ try:
+ processed_vibration = f[group+'processed'][...]
+ except KeyError:
+ pass
+ ret = [raw_vibration]
+ ret.extend(configs)
+ ret.append(processed_vibration)
+ for config in configs:
+ config.load()
+ return tuple(ret)
+
+def vib_plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
+ C=None, D=0, vibration_config=None, analyze=False):
+ """Plot 3 subfigures displaying vibration data and analysis.
-def gnuplot_check_fit(datafilename, A, B, C) :
- """
- return a string containing a gnuplot script displaying the fit.
- """
- string = 'f(x) = C / ((A**2-x**2)**2 + (x*B)**2)\n'
- string += 'A= %g ; B= %g ; C= %g\n' % (A, B, C)
- string += 'set logscale y\n'
- string += "plot '%s', f(x)\n" % datafilename
- return string
-
-@splittableKwargsFunction()
-def vib_save(data, log_dir=None) :
- """Save the dictionary data, using data_logger.data_log()
- data should be a dictionary of arrays with the fields
- 'Z piezo output'
- 'Z piezo input'
- 'Deflection input'
- """
- if log_dir :
- freq = data['sample frequency Hz']
- log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
- log_name="vib_%gHz" % freq)
- log.write_dict_of_arrays(data)
-
-def vib_load(datafile=None) :
- """Load the dictionary data, using data_logger.date_load()"
- data should be a dictionary of arrays with the fields
- 'Z piezo output'
- 'Z piezo input'
- 'Deflection input'
- """
- dl = data_logger.data_load()
- data = dl.read_dict_of_arrays(datafile)
- return data
-
-@splittableKwargsFunction()
-def vib_plot(deflection_bits, freq_axis, power, A, B, C, D,
- minFreq=None, maxFreq=None, plotVerbose=False) :
- """
- If plotVerbose or config.PYLAB_VERBOSE == True, plot 3 subfigures:
Time series (Vphoto vs sample index) (show any major drift events),
A histogram of Vphoto, with a gaussian fit (demonstrate randomness), and
FFTed Vphoto data (Vphoto vs Frequency) (show major noise components).
"""
- if plotVerbose or config.PYLAB_VERBOSE :
- common._import_pylab()
- common._pylab.figure(config.BASE_FIGNUM+2)
-
- if deflection_bits != None:
- # plot time series
- common._pylab.subplot(311)
- common._pylab.hold(False)
- common._pylab.plot(deflection_bits, 'r.')
- common._pylab.title("free oscillation")
-
- # plot histogram distribution and gaussian fit
- common._pylab.subplot(312)
- common._pylab.hold(False)
- n, bins, patches = \
- common._pylab.hist(deflection_bits, bins=30,
- normed=1, align='center')
- gaus = numpy.zeros((len(bins),), dtype=numpy.float)
- mean = deflection_bits.mean()
- std = deflection_bits.std()
- pi = numpy.pi
- exp = numpy.exp
- for i in range(len(bins)) :
- gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
- common._pylab.hold(True)
- common._pylab.plot(bins, gaus, 'r-');
- common._pylab.hold(False)
-
- # plot FFTed data
- axes = common._pylab.subplot(313)
- else:
- # use a nice big subplot just for the FFTed data
- axes = common._pylab.subplot(111)
- common._pylab.hold(False)
- common._pylab.semilogy(freq_axis, power, 'r.-')
- common._pylab.hold(True)
- xmin,xmax = axes.get_xbound()
- ymin,ymax = axes.get_ybound()
-
- if minFreq is not None and maxFreq is not None:
- # highlight the region we're fitting in
- common._pylab.axvspan(minFreq, maxFreq, facecolor='g', alpha=0.1,
- zorder = -1)
-
- fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2) + D
- common._pylab.plot(freq_axis, fitdata, 'b-')
- noisefloor = D + 0*freq_axis;
- common._pylab.plot(freq_axis, noisefloor)
- common._pylab.hold(False)
- axes.axis([xmin,xmax,ymin,ymax])
-
- common._flush_plot()
- if (plotVerbose or config.GNUPLOT_VERBOSE) and False: # TODO: cleanup and test
- # write all the ft data now
- fd = file(datafilename, 'w')
- for i in range(len(freq_axis)) :
- fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
- fd.write("\n") # blank line to separate fit data.
- # write the fit data
- for i in range(imin, imax) :
- x = freq_axis[i]
- fd.write("%g\t%g\n" % (freq_axis[i],
- C / ((A**2-x**2)**2 + (x*B)**2) ) )
- fd.close()
- print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
- g("set terminal epslatex color solid")
- g("set output 'calibration.tex'")
- g("set size 2,2") # multipliers on default 5"x3.5"
- #g("set title 'Thermal calibration'")
- g("set logscale y")
- g("set xlabel 'Frequency (Hz)'")
- g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
- g("set xrange [0:20000]")
- g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
- g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
-
- g("set mouse")
- g("pause mouse 'click with mouse'")
- g.getvar("MOUSE_BUTTON")
- del(g)
-
-# Make vib_analyze splittable (was postponed until the child functions
-# were defined).
-make_splittable_kwargs_function(vib_analyze,
- (fit_psd, 'freq_axis', 'psd_data'),
- (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C', 'D',
- 'minFreq', 'maxFreq'))
-
-@splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq'))
-def vib_load_analyze_tweaked(tweak_file, analyze=vib_analyze, **kwargs) :
- """
- Read the vib files listed in tweak_file.
- Return an array of Vphoto variances (due to the cantilever) in Volts**2
- """
- analyze_kwargs, = \
- vib_load_analyze_tweaked._splitargs(vib_load_analyze_tweaked, kwargs)
- dl = data_logger.data_load()
- Vphoto_var = []
- for path in file(tweak_file, 'r') :
- path = path.strip()
- if path[0] == '#': # a comment
- continue
- # read the data
- data = vib_load(path)
- deflection_bits = data['Deflection input']
- freq = float(path.split('_')[-1].split('Hz')[0])
- if config.TEXT_VERBOSE :
- print "Analyzing '%s' at %g Hz" % (path, freq)
- # analyze
- if 'freq' in analyze._kwargs(analyze):
- var = analyze(deflection_bits, freq, **analyze_kwargs)
- else:
- var = analyze(deflection_bits, **analyze_kwargs)
- Vphoto_var.append(var)
- return numpy.array(Vphoto_var, dtype=numpy.float)
-
-# commandline interface functions
-import scipy.io, sys
-
-def read_data(ifile):
- """
- ifile can be a filename string or open (seekable) file object.
- returns an array of data values (1 column)
- """
- #returns (column 1 array, column 2 array)
- #"""
- if ifile == None : ifile = sys.stdin
- unlabeled_data=scipy.io.read_array(ifile)
- return unlabeled_data #(unlabeled_data[:,0], unlabeled_data[:,1])
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <input-file>\n'
- '2007-2009 W. Trevor King.\n'
- '\n'
- 'There are several operation modes, each handling a different <input-file>\n'
- 'type. In each case, the output is the fitted variance (in square Volts).\n'
- '\n'
- 'Single file mode without datalogger mode (the default) :\n'
- '<input-file> should be a 1 column ASCII without a header line of cantilever\n'
- 'deflection (in bits)\n'
- 'e.g: "<deflection bits>\\n..."\n'
- '\n'
- 'Single file mode with datalogger mode :\n'
- 'Same as without datalogger mode, except the input should be a datalogger file\n'
- 'with data stored in bits (e.g. as saved by the unfold module).\n'
- '\n'
- 'Tweak file mode:\n'
- 'Runs the same analysis as in single file mode for each vib file in a tweak\n'
- 'file. Each line in the tweak file specifies a single vib file. Blank lines\n'
- 'and those beginning with a pound sign (#) are ignored. Vib files are valid\n'
- 'datalogger files as per single-file-with-datalogger-mode.\n'
- )
- parser = OptionParser(usage=usage_string, version='%prog 0.1')
- parser.add_option('-f', '--sample-freq', dest='freq', default=100e3,
- help='sample frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
- help='minimum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
- help='maximum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
- help='Output comma-seperated values (default %default)',
- default=False)
- parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
- help='Print gnuplot fit check script to stderr',
- default=False)
- parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
- help='Produce pylab fit checks during execution',
- default=False)
- parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
- help='Produce plots of the pre-fit Lorentzian',
- default=False)
- parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
- help='Run in tweak-file mode',
- default=False)
- parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
- help='Read input files with datalogger.read_dict_of_arrays(). This is useful, for example, to test a single line from a tweakfile.',
- default=False)
- parser.add_option('-s', '--disable-spectrum-fitting', dest='spectrum_fitting',
- action='store_false',
- help='Disable PSD fitting, just use the raw variance',
- default=True)
- parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
- help='Print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
- assert len(args) >= 1, "Need an input file"
-
- ifilename = args[0]
-
- if options.ofilename != None :
- ofile = file(options.ofilename, 'w')
- else :
- ofile = sys.stdout
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_INTERACTIVE = False
- config.PYLAB_VERBOSE = options.pylab
- config.GNUPLOT_VERBOSE = options.gnuplot
- PLOT_GUESSED_LORENTZIAN = options.plot_guess
- if options.spectrum_fitting == True:
- analyze = vib_analyze
- analyze_args = {'freq':options.freq,
- 'minFreq':options.min_freq,
- 'maxFreq':options.max_freq}
+ if not _matplotlib:
+ raise _matplotlib_import_error
+ figure = _matplotlib_pyplot.figure()
+
+ if power is None:
+ assert deflection != None, (
+ 'must set at least one of `deflection` or `power`.')
+ time_axes = figure.add_subplot(2, 1, 1)
+ hist_axes = figure.add_subplot(2, 1, 2)
+ freq_axes = None
+ elif deflection is None:
+ time_axes = hist_axes = None
+ freq_axes = figure.add_subplot(1, 1, 1)
else:
- analyze = vib_analyze_naive
- analyze_args = dict()
- make_splittable_kwargs_function(vib_load_analyze_tweaked,
- (vib_analyze, 'deflection_bits'))
-
- if options.tweakmode == False : # single file mode
- if options.datalogger_mode:
- data = vib_load(ifilename)
- deflection_bits = data['Deflection input']
- else:
- deflection_bits = read_data(ifilename)
- Vphoto_var = analyze(deflection_bits, **analyze_args)
- print >> ofile, Vphoto_var
- else : # tweak mode
- if 'freq' in analyze_args:
- analyze_args.pop('freq') # freq extracted from vib file names
- Vphoto_vars = vib_load_analyze_tweaked(ifilename, analyze=analyze,
- **analyze_args)
- if options.comma_out :
- sep = ','
- else :
- sep = '\n'
- common.write_array(ofile, Vphoto_vars, sep)
-
- if common._final_flush_plot != None:
- common._final_flush_plot()
+ time_axes = figure.add_subplot(3, 1, 1)
+ hist_axes = figure.add_subplot(3, 1, 2)
+ freq_axes = figure.add_subplot(3, 1, 3)
+
+ timestamp = _time.strftime('%H%M%S')
+
+ if deflection is not None:
+ time_axes.plot(deflection, 'r.')
+ time_axes.set_title('free oscillation')
+
+ # plot histogram distribution and gaussian fit
+ hist_axes.hold(True)
+ n,bins,patches = hist_axes.hist(
+ deflection, bins=30, normed=True, align='mid')
+ gauss = _numpy.zeros((len(bins),), dtype=_numpy.float)
+ mean = deflection.mean()
+ std = deflection.std()
+ pi = _numpy.pi
+ exp = _numpy.exp
+ gauss = _numpy.sqrt(2*pi)/std * exp(-0.5*((bins-mean)/std)**2)
+ # Matplotlib's normed histogram uses bin heights of n/(len(x)*dbin)
+ dbin = bins[1]-bins[0]
+ hist_axes.plot(bins, gauss/dbin, 'r-')
+ if power is not None:
+ freq_axes.hold(True)
+ freq_axes.set_yscale('log')
+ freq_axes.plot(freq_axis, power, 'r.-')
+ xmin,xmax = freq_axes.get_xbound()
+ ymin,ymax = freq_axes.get_ybound()
- if options.ofilename != None :
- ofile.close()
-
-# LocalWords: calibcant AFM
+ # highlight the region we're fitting
+ freq_axes.axvspan(
+ vibration_config['minimum-fit-frequency'],
+ vibration_config['maximum-fit-frequency'],
+ facecolor='g', alpha=0.1, zorder=-2)
+
+ if A is not None:
+ fitdata = breit_wigner(freq_axis, A, B, C, D)
+ freq_axes.plot(freq_axis, fitdata, 'b-')
+ noisefloor = D + 0*freq_axis;
+ freq_axes.plot(freq_axis, noisefloor)
+
+ if B**2 < 2*A**2:
+ res_freq = _numpy.sqrt(A**2 - B**2/2)
+ freq_axes.axvline(res_freq, color='b', zorder=-1)
+
+ freq_axes.set_title('power spectral density %s' % timestamp)
+ freq_axes.axis([xmin,xmax,ymin,ymax])
+ freq_axes.set_xlabel('frequency (Hz)')
+ freq_axes.set_ylabel('PSD')
+
+ figure.show()
projects / calibcant.git / blobdiff