[FEATURE] add pytta.vibration with human vibration weightings

pytta/__init__.py

@@ -108,6 +108,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import vibration
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +131,7 @@
            'generate',
            'utils',
            'iso3741',
+           'vibration',
            'rooms',
 
            # Functions

pytta/vibration.py

@@ -0,0 +1,122 @@
+"""
+PyTTa human vibration:
+---------------------
+Frequency weightings of whole-body vibration (ISO 2631-1, ISO 2631-2,
+ISO 2631-4) and hand-arm vibration (ISO 5349-1), as specified for
+instruments in ISO 8041, the evaluation values built on them (weighted
+rms, maximum transient vibration value, vibration dose value), and
+integration and differentiation of signals (acceleration, velocity,
+displacement).
+
+    >>> aw = pytta.vibration.rms(acceleration, 'Wk')      # m/s2
+
+"""
+import numpy as np
+import scipy.signal as ss
+
+from pytta.classes import SignalObj
+
+_INF = np.inf
+# ISO 8041 parameters: band limiting f1, f2 (Q = 1/sqrt(2)), a-v
+# transition f3, f4, Q4, upward step f5, Q5, f6, Q6, gain K.
+_PARAMETERS = {
+    'Wb': (0.4, 100, 16, 16, 0.55, 2.5, 0.9, 4, 0.95, 1.024),
+    'Wc': (0.4, 100, 8, 8, 0.63, _INF, 1, _INF, 1, 1),
+    'Wd': (0.4, 100, 2, 2, 0.63, _INF, 1, _INF, 1, 1),
+    'We': (0.4, 100, 1, 1, 0.63, _INF, 1, _INF, 1, 1),
+    'Wf': (0.08, 0.63, _INF, 0.25, 0.86, 0.0625, 0.80, 0.10, 0.80, 1),
+    'Wh': (10 ** 0.8, 10 ** 3.1, 100 / (2 * np.pi), 100 / (2 * np.pi), 0.64,
+           _INF, 1, _INF, 1, 1),
+    'Wj': (0.4, 100, _INF, _INF, 1, 3.75, 0.91, 5.32, 0.91, 1),
+    'Wk': (0.4, 100, 12.5, 12.5, 0.63, 2.37, 0.91, 3.35, 0.91, 1),
+    'Wm': (10 ** -0.1, 100, 1 / (0.028 * 2 * np.pi), 1 / (0.028 * 2 * np.pi),
+           0.5, _INF, 1, _INF, 1, 1),
+}
+
+
+def weighting_filter(kind, fs):
+    """
+    Second-order sections of the frequency weighting `kind` ('Wb', 'Wc',
+    'Wd', 'We', 'Wf', 'Wh', 'Wj', 'Wk' or 'Wm') at the sampling rate fs:
+    the product of the band-limiting high-pass and low-pass (second-order
+    Butterworth), the acceleration-velocity transition and the upward step
+    of ISO 8041, through the bilinear transform.
+    """
+    if kind not in _PARAMETERS:
+        raise ValueError(f'kind must be one of {sorted(_PARAMETERS)}.')
+    f1, f2, f3, f4, Q4, f5, Q5, f6, Q6, K = _PARAMETERS[kind]
+    w = lambda f: 2 * np.pi * f
+    q = 1 / np.sqrt(2)
+    # each factor as zeros, poles and gain, so that the low-frequency
+    # poles keep their precision (a single polynomial does not)
+    zeros, poles = [0.0, 0.0], []
+    poles += list(np.roots([1, w(f1) / q, w(f1) ** 2]))          # Hh
+    poles += list(np.roots([1, w(f2) / q, w(f2) ** 2]))          # Hl
+    gain = K * w(f2) ** 2
+    if np.isfinite(f4):                                          # Ht
+        poles += list(np.roots([1, w(f4) / Q4, w(f4) ** 2]))
+        gain *= w(f4) ** 2
+        if np.isfinite(f3):
+            zeros.append(-w(f3))
+            gain /= w(f3)
+    if np.isfinite(f5):                                          # Hs
+        zeros += list(np.roots([1, w(f5) / Q5, w(f5) ** 2]))
+        poles += list(np.roots([1, w(f6) / Q6, w(f6) ** 2]))
+    return ss.zpk2sos(*ss.bilinear_zpk(zeros, poles, gain, fs))
+
+
+def _weighted(signalObj, kind):
+    x = np.asarray(signalObj.timeSignal, dtype=float)
+    if kind is None:
+        return x
+    return ss.sosfilt(weighting_filter(kind, signalObj.samplingRate), x,
+                      axis=0)
+
+
+def rms(signalObj, kind=None):
+    """Frequency-weighted rms value of each channel (no weighting with
+    kind None)."""
+    return np.sqrt(np.mean(_weighted(signalObj, kind) ** 2, axis=0))
+
+
+def mtvv(signalObj, kind=None, tau=1.0):
+    """Maximum transient vibration value: the maximum of the running rms
+    with exponential averaging of time constant tau (1 s)."""
+    a = _weighted(signalObj, kind)
+    alpha = np.exp(-1 / (tau * signalObj.samplingRate))
+    running = ss.lfilter([1 - alpha], [1, -alpha], a ** 2, axis=0)
+    return np.sqrt(np.max(running, axis=0))
+
+
+def vdv(signalObj, kind=None):
+    """Fourth power vibration dose value, (integral of aw^4 dt)^(1/4)."""
+    a = _weighted(signalObj, kind)
+    return (np.sum(a ** 4, axis=0) / signalObj.samplingRate) ** 0.25
+
+
+def _spectral(signalObj, power, fmin):
+    x = np.asarray(signalObj.timeSignal, dtype=float)
+    n = x.shape[0]
+    f = np.fft.rfftfreq(n, 1 / signalObj.samplingRate)
+    factor = np.zeros(f.size, dtype=complex)
+    keep = f > max(fmin, 0)
+    factor[keep] = (2j * np.pi * f[keep]) ** power
+    if power >= 0:
+        factor[0] = 0 if power > 0 else 1
+    y = np.fft.irfft(np.fft.rfft(x, axis=0) * factor[:, None], n, axis=0)
+    return SignalObj(y, 'time', samplingRate=signalObj.samplingRate,
+                     signalType=signalObj.signalType)
+
+
+def integrate(signalObj, times=1, fmin=0.):
+    """Integrate `times` times in the frequency domain (divide by j 2 pi f),
+    e.g. acceleration to velocity. DC and the bins up to fmin are set to 0,
+    so the result has no drift; the signal is taken as periodic."""
+    return _spectral(signalObj, -times, fmin)
+
+
+def differentiate(signalObj, times=1):
+    """Differentiate `times` times in the frequency domain (multiply by
+    j 2 pi f), e.g. displacement to velocity. The signal is taken as
+    periodic."""
+    return _spectral(signalObj, times, 0.)

tests/test_vibration.py

@@ -0,0 +1,68 @@
+"""
+Tests for pytta.vibration: frequency weightings of human vibration
+(ISO 2631-1, ISO 5349-1, parameters of ISO 8041), the weighted rms, MTVV
+and VDV, and integration and differentiation. Each test pins one rule.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+import scipy.signal as ss
+
+import pytta
+from pytta import vibration
+
+# Weightings in dB at nominal frequencies, computed from the ISO 8041
+# parameters (tabulated to 0.01 dB in the Virtins note "Frequency
+# weightings for sound and vibration perceived by humans").
+TABLE = {'Wk': {0.1: -30.11, 1: -6.33, 4: -0.29, 8: 0.31, 16: -2.28,
+                63: -14.61},
+         'Wd': {0.1: -24.09, 1: 0.10, 2: -1.01, 8: -11.93, 31.5: -23.98},
+         'Wh': {8: -1.18, 16: -0.96, 100: -15.91, 1000: -37.42}}
+
+
+class TestVibration(unittest.TestCase):
+
+    def test_weightings(self):
+        for kind, values in TABLE.items():
+            fs = 48000 if kind == 'Wh' else 2000
+            sos = vibration.weighting_filter(kind, fs)
+            for f, expected in values.items():
+                _, h = ss.sosfreqz(sos, worN=np.array([float(f)]), fs=fs)
+                self.assertAlmostEqual(20 * np.log10(abs(h[0])), expected,
+                                       delta=0.05, msg=(kind, f))
+        with self.assertRaises(ValueError):
+            vibration.weighting_filter('Wz', 2000)
+
+    def test_rms_mtvv_and_vdv_of_a_sine(self):
+        # Unweighted sine of amplitude A over T seconds: rms A / sqrt(2),
+        # MTVV the same once settled, VDV A (3 T / 8)^(1/4).
+        fs, A, T = 2000, 2.0, 20.0
+        t = np.arange(int(T * fs)) / fs
+        sig = pytta.SignalObj(A * np.sin(2 * np.pi * 10 * t), 'time', fs)
+        self.assertAlmostEqual(vibration.rms(sig)[0], A / np.sqrt(2), 4)
+        self.assertAlmostEqual(vibration.mtvv(sig)[0], A / np.sqrt(2),
+                               delta=0.01)
+        self.assertAlmostEqual(vibration.vdv(sig)[0],
+                               A * (3 * T / 8) ** 0.25, 4)
+
+    def test_integration_and_differentiation_of_a_sine(self):
+        # A whole number of cycles of sin(w t): the integral is
+        # -cos(w t) / w and the derivative w cos(w t).
+        fs = 1000
+        t = np.arange(fs) / fs
+        w = 2 * np.pi * 50
+        sig = pytta.SignalObj(np.sin(w * t), 'time', fs)
+        integral = vibration.integrate(sig).timeSignal[:, 0]
+        derivative = vibration.differentiate(sig).timeSignal[:, 0]
+        self.assertTrue(np.allclose(integral, -np.cos(w * t) / w, atol=1e-6))
+        self.assertTrue(np.allclose(derivative / w, np.cos(w * t), atol=1e-5))
+        twice = vibration.integrate(sig, times=2).timeSignal[:, 0]
+        self.assertTrue(np.allclose(twice, -np.sin(w * t) / w ** 2,
+                                    atol=1e-8))
+
+
+if __name__ == '__main__':
+    unittest.main()