From: W. Trevor King <wking@tremily.us>
Date: Sun, 5 May 2013 00:13:21 +0000 (-0400)
Subject: vibration_analyze: Add smooth_window to PSD fitting
X-Git-Tag: 0.9~3
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=3418c8c094616215515bf0a6f349ba06e450b16a;p=calibcant.git

vibration_analyze: Add smooth_window to PSD fitting

Sometimes the broad resonant peak from the cantilever has sharper
unidentified peaks sitting on its shoulders.  In the rare case that
these peaks are taller than the cantilever peak, they can throw off
the resonant frequency portion of the parameter-guessing heuristic.
By smoothing the PSD data somewhat, we can damp down the sharp peaks
and bring the cantilever peak back to the top.

Because the amount of smoothing you'll need depends on your noise
sources and the expected width of the cantilever peak, make the
smoothing configurable.  A 10-sample-wide Hann filter (rescaled to
have unit weight) seems like a pretty conservative default.
---

diff --git a/calibcant/config.py b/calibcant/config.py
index 87ec917..11686f7 100644
--- a/calibcant/config.py
+++ b/calibcant/config.py
@@ -157,6 +157,17 @@ class VibrationConfig (_config.Config):
             name='maximum-fit-frequency',
             help='Upper bound of Lorentzian fitting region.',
             default=25e3),
+        _config.ChoiceSetting(
+            name='smooth-window',
+            help='Smoothing window type (for guessing PSD parameters).',
+            default=_window_hann,
+            choices=[
+                ('Hann', _window_hann),
+                ]),
+        _config.FloatSetting(
+            name='smooth-length',
+            help='Size of the smoothing window (in points).',
+            default=10),
         ]
 
 
diff --git a/calibcant/vibration_analyze.py b/calibcant/vibration_analyze.py
index f5b83de..a0d383e 100644
--- a/calibcant/vibration_analyze.py
+++ b/calibcant/vibration_analyze.py
@@ -306,10 +306,17 @@ def analyze(deflection, config, deflection_channel_config,
         config['chunk-size'], config['overlap'],
         config['window'])
 
+    if config['smooth-window'] and config['smooth-length']:
+        n = config['smooth-length']
+        smooth_window = config['smooth-window'](length=n)
+        smooth_window /= smooth_window.sum()
+    else:
+        smooth_window = None
     A,B,C,D = fit_psd(
         freq_axis, power,
         min_frequency=config['minimum-fit-frequency'],
         max_frequency=config['maximum-fit-frequency'],
+        smooth_window=smooth_window,
         offset=config['model'] == _OffsetBreitWigner)
 
     _LOG.debug('fit PSD(f) = C / ((A**2-f**2)**2 + (f*B)**2) with '
@@ -340,7 +347,7 @@ def breit_wigner(f, A, B, C, D=0):
     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):
+            smooth_window=_FFT_tools.window_hann(10), offset=False):
     """Fit the FFTed vibration data to a Breit-Wigner.
 
     Inputs
@@ -348,6 +355,7 @@ def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
       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
+      smooth_window  array for smoothing the PSD before guessing params
       offset         add a white-noise offset to the Breit-Wigner model
     Output
       Breit-Wigner model fit parameters `A`, `B`, `C`, and `D`.
@@ -363,7 +371,11 @@ def fit_psd(freq_axis, psd_data, min_frequency=500, max_frequency=25000,
         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
+    if hasattr(smooth_window, '__len__'):
+        smooth_psd_data = _numpy.convolve(psd_data, smooth_window, mode='same')
+    else:
+        smooth_psd_data = psd_data
+    max_psd_index = _numpy.argmax(smooth_psd_data[imin:imax]) + imin
     max_psd = psd_data[max_psd_index]
     res_freq = freq_axis[max_psd_index]