-#!/usr/bin/python
-#
# calibcant - tools for thermally calibrating AFM cantilevers
#
-# Copyright (C) 2007,2008, William Trevor King
+# Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
#
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License as
-# published by the Free Software Foundation; either version 3 of the
-# License, or (at your option) any later version.
+# This file is part of calibcant.
#
-# This program is distributed in the hope that it will be useful, but
-# WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-# See the GNU General Public License for more details.
+# calibcant is free software: you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation, either
+# version 3 of the License, or (at your option) any later version.
#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
-# 02111-1307, USA.
+# calibcant is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
#
-# The author may be contacted at <wking@drexel.edu> on the Internet, or
-# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
-# Philadelphia PA 19104, USA.
+# You should have received a copy of the GNU Lesser General Public
+# 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 VibrationConfig
+>>> from h5config.storage.hdf5 import pprint_HDF5
+>>> from pypiezo.test import get_piezo_config
+>>> from pypiezo.base import convert_volts_to_bits
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> piezo_config = get_piezo_config()
+>>> vibration_config = VibrationConfig()
+>>> vibration_config['frequency'] = 50e3
+
+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...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(
+... piezo_config.select_config('inputs', 'deflection'), 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,
+... piezo_config.select_config('inputs', 'deflection'))
+>>> 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=piezo_config.select_config(
+... 'inputs', 'deflection'),
+... processed_vibration=processed_vibration)
+
+>>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
+/
+ /vibration
+ /vibration/config
+ /vibration/config/deflection
+ <HDF5 dataset "channel": shape (), type "<i4">
+ 0
+ <HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
+ [0 1]
+ <HDF5 dataset "conversion-origin": shape (), type "<i4">
+ 0
+ <HDF5 dataset "device": shape (), type "|S12">
+ /dev/comedi0
+ <HDF5 dataset "inverse-conversion-coefficients": shape (2,), type "<i4">
+ [0 1]
+ <HDF5 dataset "inverse-conversion-origin": shape (), type "<i4">
+ 0
+ <HDF5 dataset "maxdata": shape (), type "<i4">
+ 100
+ <HDF5 dataset "name": shape (), type "|S10">
+ deflection
+ <HDF5 dataset "range": shape (), type "<i4">
+ 1
+ <HDF5 dataset "subdevice": shape (), type "<i4">
+ -1
+ /vibration/config/vibration
+ <HDF5 dataset "chunk-size": shape (), type "<i4">
+ 2048
+ <HDF5 dataset "frequency": shape (), type "<f8">
+ 50000.0
+ <HDF5 dataset "maximum-fit-frequency": shape (), type "<f8">
+ 25000.0
+ <HDF5 dataset "minimum-fit-frequency": shape (), type "<f8">
+ 500.0
+ <HDF5 dataset "model": shape (), type "|S12">
+ Breit-Wigner
+ <HDF5 dataset "overlap": shape (), type "|b1">
+ False
+ <HDF5 dataset "sample-time": shape (), type "<i4">
+ 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 # used for fitting lorentzian
-import common # common module for the calibcant package
-import config # config module for the calibcant package
-import data_logger
-import FFT_tools
+import os as _os
+import time as _time
-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.
-
- Vphoto_in2V is a function converting Vphoto values from bits to Volts.
- This function is assumed linear, see bump_analyze().
- zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
- """
- Vphoto_std_bits = deflection_bits.std()
- Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in)
- return Vphoto_std**2
-
-def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000,
- chunk_size=4096, overlap=True, window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V, plotVerbose=False) :
- """
- Calculate the variance in the raw data due to the cantilever's
- thermal oscillation and convert it to Volts**2.
+import h5py as _h5py
+import numpy as _numpy
+from scipy.optimize import leastsq as _leastsq
+
+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
+
+from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
+from h5config.storage.hdf5 import h5_create_group as _h5_create_group
+import FFT_tools as _FFT_tools
+from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
+from pypiezo.config import ChannelConfig as _ChannelConfig
- Improves on vib_analyze_naive() by first converting Vphoto(t) to
- frequency space, and fitting a Lorentzian in the relevant frequency
- range (see cantilever_calib.pdf for derivation).
+from . import LOG as _LOG
+from . import package_config as _package_config
+from .config import Variance as _Variance
+from .config import BreitWigner as _BreitWigner
+from .config import OffsetBreitWigner as _OffsetBreitWigner
+from .config import VibrationConfig as _VibrationConfig
- 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.avg_power_spectrum
+
+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().
- zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts.
+ Inputs
+ deflection : numpy array with deflection timeseries in Volts.
+ """
+ std = deflection.std()
+ var = std**2
+ _LOG.debug('naive deflection variance: %g V**2' % var)
+ return var
+
+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.
+
+ 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.
"""
- 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])
-
+ 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 = fit_psd(freq_axis, power, minFreq, maxFreq)
+ 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" % (A, B, C)
+ _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, plotVerbose=plotVerbose)
+ if plot or _package_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
- return lorentzian_area(A,B,C) * 1e-12
+ fitted_variance = breit_wigner_area(A,B,C) * 1e-12
-def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) :
- """
- freq_axis : array of frequencies in Hz
- psd_data : array of PSD amplitudes for the frequencies in freq_axis
+ _LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
- Calculate the variance in the raw data due to the cantilever's
- thermal oscillation and convert it to Volts**2.
+ if _package_config['matplotlib']:
+ vib_plot(deflection, freq_axis, power, A, B, C, D,
+ vibration_config=vibration_config)
- Improves on vib_analyze_naive() by working on frequency space data
- and fitting a Lorentzian in the relevant frequency range (see
- cantilever_calib.pdf for derivation).
+ return min(fitted_variance, naive_variance)
- 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().
+def breit_wigner(f, A, B, C, D=0):
+ """Breit-Wigner (sortof).
+
+ Inputs
+ f Frequency
+ A Resonant frequency
+ B Quality Q = A/B
+ C Scaling factor
+ D Optional white-noise offset
+
+ 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)
-
- # 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()
-
- # generate guesses for Lorentzian parameters A,B,C
- max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin
+ 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 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 = 5, and adjust from there
- Q_guess = 5
+ # 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)
- #
+ 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
# (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2]
# so w is a disaster ;)
# For some values of A and B (non-underdamped), W is imaginary or negative.
- if config.TEXT_VERBOSE :
- print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess)
-
- # fit Lorentzian using Gnuplot's 'fit' command
- g = GnuplotBiDir.Gnuplot()
- # The Lorentzian is the predicted one-sided power spectral density
- # For an oscillator whose position x obeys
- # m x'' + gamma x' + kx = F_thermal(t)
- # A is the ?driving noise?
- # B is gamma, the damping term
- # C is omega_0, the resonant frequency
+ #
+ # 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
+
# Fitting the PSD of V = photoSensitivity*x just rescales the parameters
- g("f(x) = C / ((A**2-x**2)**2 + (B*x)**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
- 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)
- return (A, B, C)
-
-def lorentzian_area(A, B, C) :
+
+ # 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:
+ A,B,C = p
+ 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.
+ C=abs(C) # Only abs(C) is used in breit_wigner().
+ return (A, B, C, D)
+
+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
+ storage = _HDF5_Storage()
+ 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)
+ storage.save(config=config, 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 [(_VibrationConfig, 'config/vibration'),
+ (_ChannelConfig, 'config/deflection')]:
+ config = Config(storage=_HDF5_Storage(
+ filename=filename, group=group+key))
+ configs.append(config)
+ try:
+ processed_vibration = float(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
-
-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 :
- 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
-
-def vib_plot(deflection_bits, freq_axis, power, A, B, C,
- plotVerbose=True) :
- """
- 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 :
- print "plotting"
- common._import_pylab()
- common._pylab.figure(config.BASE_FIGNUM+2)
- common._pylab.hold(False)
+ 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:
+ 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')
- # plot time series
- common._pylab.subplot(311)
- common._pylab.plot(deflection_bits, 'r.')
- common._pylab.title("free oscillation")
+ if deflection is not None:
+ time_axes.plot(deflection, 'r.')
+ time_axes.set_title('free oscillation')
# plot histogram distribution and gaussian fit
- common._pylab.subplot(312)
- n, bins, patches = \
- common._pylab.hist(deflection_bits, bins=30,
- normed=1, align='center')
- gaus = numpy.zeros((len(bins),))
- 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
- common._pylab.subplot(313)
- common._pylab.semilogy(freq_axis, power, 'r.-')
- common._flush_plot()
- if (plotVerbose or config.GNUPLOT_VERBOSE):
- # 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)
-
-def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
- chunk_size=2048, overlap=True,
- window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V,
- plotVerbose=False) :
- """
- Read the vib files listed in tweak_file.
- Return an array of Vphoto variances (due to the cantilever) in Volts**2
- """
- dl = data_logger.data_load()
- Vphoto_var = []
- for path in file(tweak_file, 'r') :
- if path[-1] == '\n':
- path = path.split('\n')[0]
- # 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
- Vphoto_var.append(vib_analyze(deflection_bits, freq, minFreq, maxFreq,
- chunk_size, overlap, window,
- Vphoto_in2V, plotVerbose))
- 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'
- '2008, W. Trevor King.\n'
- '\n'
- 'Deflection power spectral density (PSD) data (X^2/Hz)\n'
- 'and returns the variance in X units (X^2)'
- '<input-file> should be whitespace-delimited, 2 column ASCII\n'
- 'without a header line.\n'
- 'e.g: "<frequency_Hz>\\t<deflection_psd X^2/Hz>\\n"\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('-t', '--tweak-mode', dest='tweakmode', action='store_true',
- help='Run in tweak-file mode',
- default=False)
- 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
- if options.verbose == True :
- vfile = sys.stderr
- else :
- vfile = None
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_INTERACTIVE = False
- config.PYLAB_VERBOSE = options.pylab
- config.GNUPLOT_VERBOSE = options.gnuplot
-
- if options.tweakmode == False :
- data = read_data(ifilename)
- deflection_bits = data['Deflection input']
- Vphoto_var = vib_analyze(deflection_bits, freq=options.freq,
- minFreq=options.min_freq,
- maxFreq=options.max_freq)
- print >> ofile, Vphoto_var
- else : # tweak mode
- Vphoto_vars = vib_load_analyze_tweaked(ifilename,
- minFreq=options.min_freq,
- maxFreq=options.max_freq)
- 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()
-
- if options.ofilename != None :
- ofile.close()
+ 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()
+ # 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