"""Thermal vibration analysis.
-Separate the more general `vib_analyze()` from the other `vib_*()`
+Separate the more general `analyze()` from the other `vibration`
functions in calibcant.
The relevent physical quantities are :
>>> os.close(fd)
>>> piezo_config = get_piezo_config()
->>> vibration_config = VibrationConfig()
->>> vibration_config['frequency'] = 50e3
+>>> config = VibrationConfig()
+>>> 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
>>> m = 5e-11 # kg
>>> gamma = 1.6e-6 # N*s/m
>>> k = 0.05 # N/m
->>> T = 1/vibration_config['frequency']
+>>> T = 1/config['frequency']
>>> T # doctest: +ELLIPSIS
2...e-05
>>> N = int(2**15) # count
>>> f0 = w0/(2*numpy.pi)
>>> f0 # doctest: +ELLIPSIS
5032.9...
->>> f_nyquist = vibration_config['frequency']/2
+>>> f_nyquist = config['frequency']/2
>>> f_nyquist # doctest: +ELLIPSIS
25000.0...
Analyze the simulated data.
->>> naive_vibration = vib_analyze_naive(deflection)
->>> naive_vibration # doctest: +SKIP
+>>> naive = analyze_naive(deflection)
+>>> naive # doctest: +SKIP
136871517.43486854
->>> abs(naive_vibration / 136.9e6 - 1) < 0.1
+>>> abs(naive / 136.9e6 - 1) < 0.1
True
->>> processed_vibration = vib_analyze(
-... deflection_bits, vibration_config,
+>>> processed = analyze(
+... deflection_bits, config,
... piezo_config.select_config('inputs', 'deflection'))
->>> processed_vibration # doctest: +SKIP
+>>> processed # 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,
+>>> plot(deflection=deflection_bits, config=config)
+>>> save(filename=filename, group='/vibration/',
+... raw=deflection_bits, config=config,
... deflection_channel_config=piezo_config.select_config(
... 'inputs', 'deflection'),
-... processed_vibration=processed_vibration)
+... processed=processed)
>>> pprint_HDF5(filename) # doctest: +ELLIPSIS, +REPORT_UDIFF
/
/vibration
/vibration/config
/vibration/config/deflection
+ <HDF5 dataset "analog-reference": shape (), type "|S6">
+ ground
<HDF5 dataset "channel": shape (), type "<i4">
0
<HDF5 dataset "conversion-coefficients": shape (2,), type "<i4">
<HDF5 dataset "name": shape (), type "|S10">
deflection
<HDF5 dataset "range": shape (), type "<i4">
- 1
+ 0
<HDF5 dataset "subdevice": shape (), type "<i4">
-1
/vibration/config/vibration
1
<HDF5 dataset "window": shape (), type "|S4">
Hann
- <HDF5 dataset "processed": shape (), type "<f8">
- ...
+ /vibration/processed
+ <HDF5 dataset "data": shape (), type "<f8">
+ ...
+ <HDF5 dataset "units": shape (), type "|S6">
+ V^2/Hz
/vibration/raw
- <HDF5 dataset "deflection": shape (32768,), type "<f8">
+ <HDF5 dataset "data": shape (32768,), type "<f8">
[...]
+ <HDF5 dataset "units": shape (), type "|S4">
+ bits
->>> (raw_vibration,vibration_config,deflection_channel_config,
-... processed_vibration) = vib_load(
-... filename=filename, group='/vibration/')
+>>> data = load(filename=filename, group='/vibration/')
->>> processed_vibration # doctest: +SKIP
+>>> pprint(data) # doctest: +ELLIPSIS
+{'config': {'vibration': <InputChannelConfig ...>},
+ 'processed': ...,
+ 'raw': array([...])}
+>>> data['processed'] # doctest: +SKIP
136457906.25574699
->>> abs(processed_vibration / 136.5e6 - 1) < 0.1
+>>> abs(data['processed'] / 136.5e6 - 1) < 0.1
True
>>> os.remove(filename)
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
+from pypiezo.config import InputChannelConfig as _InputChannelConfig
from . import LOG as _LOG
from . import package_config as _package_config
from .config import BreitWigner as _BreitWigner
from .config import OffsetBreitWigner as _OffsetBreitWigner
from .config import VibrationConfig as _VibrationConfig
+from .util import SaveSpec as _SaveSpec
+from .util import save as _save
+from .util import load as _load
-def vib_analyze_naive(deflection):
+def analyze_naive(deflection):
"""Calculate the deflection variance in Volts**2.
This method is simple and easy to understand, but it highly
_LOG.debug('naive deflection variance: %g V**2' % var)
return var
-def vib_analyze(deflection, vibration_config, deflection_channel_config,
- plot=False):
+def analyze(deflection, config, deflection_channel_config,
+ plot=False):
"""Calculate the deflection variance in Volts**2.
- Improves on `vib_analyze_naive()` by first converting `Vphoto(t)`
+ Improves on `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
Inputs:
deflection Vphoto deflection input in bits.
- vibration_config `.config._VibrationConfig` instance
+ config `.config.VibrationConfig` instance
deflection_channel_config
deflection `pypiezo.config.ChannelConfig` instance
plot boolean overriding matplotlib config setting.
mean = deflection_v.mean()
_LOG.debug('average thermal deflection (Volts): %g' % mean)
- naive_variance = vib_analyze_naive(deflection_v)
- if vibration_config['model'] == _Variance:
+ naive_variance = analyze_naive(deflection_v)
+ if config['model'] == _Variance:
return naive_variance
# Compute the average power spectral density per unit time (in uV**2/Hz)
_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'])
+ (deflection_v - mean)*1e6, config['frequency'],
+ config['chunk-size'], config['overlap'],
+ config['window'])
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)
+ min_frequency=config['minimum-fit-frequency'],
+ max_frequency=config['maximum-fit-frequency'],
+ offset=config['model'] == _OffsetBreitWigner)
_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))
if plot or _package_config['matplotlib']:
- vib_plot(deflection, freq_axis, power, A, B, C, D,
- vibration_config=vibration_config)
+ _plot(deflection, freq_axis, power, A, B, C, D, config=config)
# Our A is in uV**2, so convert back to Volts**2
fitted_variance = breit_wigner_area(A,B,C) * 1e-12
_LOG.debug('fitted deflection variance: %g V**2' % fitted_variance)
if _package_config['matplotlib']:
- vib_plot(deflection, freq_axis, power, A, B, C, D,
- vibration_config=vibration_config)
+ plot(deflection, freq_axis, power, A, B, C, D,
+ config=config)
return min(fitted_variance, naive_variance)
# <V(t)**2> = (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):
+def save(filename, group='/', raw=None, config=None,
+ deflection_channel_config=None, processed=None):
+ specs = [
+ _SaveSpec(item=config, relpath='config/vibration',
+ config=_VibrationConfig),
+ _SaveSpec(item=deflection_channel_config,
+ relpath='config/deflection',
+ config=_InputChannelConfig),
+ _SaveSpec(item=raw, relpath='raw', units='bits'),
+ _SaveSpec(item=processed, relpath='processed', units='V^2/Hz'),
+ ]
+ _save(filename=filename, group=group, specs=specs)
+
+def load(filename=None, group='/'):
+ specs = [
+ _SaveSpec(key=('config', 'vibration'), relpath='config/vibration',
+ config=_VibrationConfig),
+ _SaveSpec(key=('config', 'deflection'), relpath='config/deflection',
+ config=_InputChannelConfig),
+ _SaveSpec(key=('raw',), relpath='raw', array=True, units='bits'),
+ _SaveSpec(key=('processed',), relpath='processed', units='V^2/Hz'),
+ ]
+ return _load(filename=filename, group=group, specs=specs)
+
+def plot(deflection=None, freq_axis=None, power=None, A=None, B=None,
+ C=None, D=0, config=None, analyze=False):
"""Plot 3 subfigures displaying vibration data and analysis.
Time series (Vphoto vs sample index) (show any major drift events),
# highlight the region we're fitting
freq_axes.axvspan(
- vibration_config['minimum-fit-frequency'],
- vibration_config['maximum-fit-frequency'],
+ config['minimum-fit-frequency'],
+ config['maximum-fit-frequency'],
facecolor='g', alpha=0.1, zorder=-2)
if A is not None:
freq_axes.set_ylabel('PSD')
if hasattr(figure, 'show'):
figure.show()
+_plot = plot # alternative name for use inside analyze()
projects / calibcant.git / blobdiff