[FEATURE] add pytta.iec61672 with the sound level meter measures

pytta/__init__.py

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

pytta/iec61672.py

@@ -0,0 +1,113 @@
+"""
+PyTTa IEC 61672-1 sound level meter:
+-----------------------------------
+Frequency weightings A, C and Z, time weightings F and S, and the levels
+of a sound level meter of IEC 61672-1:2013, plus the percentile levels of a
+time-weighted level.
+
+    >>> pytta.iec61672.leq(signalObj, weighting='A')    # LAeq, dB
+
+Levels are in dB re the dBRef of each channel (20 uPa for a channel in Pa)
+and are returned with one value per channel.
+"""
+import numpy as np
+import scipy.signal as ss
+
+_TIME_CONSTANTS = {'F': 0.125, 'S': 1.0}            # 5.8.1
+
+
+def _poles():
+    """Pole frequencies f1 to f4 in Hz, from Annex E (E.2.4, E.3.4)."""
+    fr, fL, fH, D = 1000.0, 10 ** 1.5, 10 ** 3.9, 1 / np.sqrt(2)
+    b = (fr ** 2 + fL ** 2 * fH ** 2 / fr ** 2 - D * (fL ** 2 + fH ** 2)) \
+        / (1 - D)
+    c = fL ** 2 * fH ** 2
+    f1 = np.sqrt((-b - np.sqrt(b ** 2 - 4 * c)) / 2)
+    f4 = np.sqrt((-b + np.sqrt(b ** 2 - 4 * c)) / 2)
+    fA = 10 ** 2.45
+    return f1, (3 - 5 ** 0.5) / 2 * fA, (3 + 5 ** 0.5) / 2 * fA, f4
+
+
+def weighting_filter(kind, fs):
+    """
+    Second-order sections of the frequency weighting 'A', 'C' or 'Z' at
+    the sampling rate fs: the poles and zeros of Annex E through the
+    bilinear transform, with 0 dB at 1 kHz. 'Z' returns None (no filter).
+    """
+    if kind not in ('A', 'C', 'Z'):
+        raise ValueError("weighting must be 'A', 'C' or 'Z'.")
+    if kind == 'Z':
+        return None
+    f1, f2, f3, f4 = _poles()
+    poles = [f1, f1, f4, f4] + ([f2, f3] if kind == 'A' else [])
+    zeros = [0.0] * (4 if kind == 'A' else 2)
+    sos = ss.zpk2sos(*ss.bilinear_zpk(zeros, [-2 * np.pi * p for p in poles],
+                                      1.0, fs))
+    _, h = ss.sosfreqz(sos, worN=np.array([1000.0]), fs=fs)
+    sos[0, :3] /= abs(h[0])
+    return sos
+
+
+def _weighted(signalObj, weighting):
+    """Frequency-weighted samples (samples, channels) and the dB reference
+    of each channel."""
+    x = np.asarray(signalObj.timeSignal, dtype=float)
+    sos = weighting_filter(weighting, signalObj.samplingRate)
+    if sos is not None:
+        x = ss.sosfilt(sos, x, axis=0)
+    refs = np.array([signalObj.channels[ch].dBRef
+                     for ch in signalObj.channels.mapping], dtype=float)
+    return x, refs
+
+
+def time_weighted_level(signalObj, weighting='A', time='F'):
+    """
+    Exponentially time-weighted sound level of 5.8 (F: 0.125 s, S: 1 s),
+    at every sample: (time in s, level in dB of shape (samples, channels)).
+    """
+    if time not in _TIME_CONSTANTS:
+        raise ValueError("time must be 'F' or 'S'.")
+    x, refs = _weighted(signalObj, weighting)
+    alpha = np.exp(-1 / (_TIME_CONSTANTS[time] * signalObj.samplingRate))
+    meanSquare = ss.lfilter([1 - alpha], [1, -alpha], x ** 2, axis=0)
+    with np.errstate(divide='ignore'):
+        level = 10 * np.log10(meanSquare / refs ** 2)
+    return np.arange(x.shape[0]) / signalObj.samplingRate, level
+
+
+def leq(signalObj, weighting='A'):
+    """Time-averaged (equivalent continuous) sound level over the signal."""
+    x, refs = _weighted(signalObj, weighting)
+    return 10 * np.log10(np.mean(x ** 2, axis=0) / refs ** 2)
+
+
+def sel(signalObj, weighting='A'):
+    """Sound exposure level, the time integral of the squared weighted
+    signal re the reference squared times 1 s."""
+    x, refs = _weighted(signalObj, weighting)
+    return 10 * np.log10(np.sum(x ** 2, axis=0) / signalObj.samplingRate
+                         / refs ** 2)
+
+
+def lmax(signalObj, weighting='A', time='F'):
+    """Maximum time-weighted sound level."""
+    return np.max(time_weighted_level(signalObj, weighting, time)[1], axis=0)
+
+
+def lpeak(signalObj, weighting='C'):
+    """Peak sound level, from the largest magnitude of the weighted
+    signal (C by default, as LCpeak)."""
+    x, refs = _weighted(signalObj, weighting)
+    return 20 * np.log10(np.max(np.abs(x), axis=0) / refs)
+
+
+def percentile_levels(signalObj, percents=(10, 50, 90), weighting='A',
+                      time='F'):
+    """
+    Percentile levels L_N: the time-weighted level exceeded during N % of
+    the signal, for each N in percents. Returns an array of shape
+    (len(percents), channels).
+    """
+    level = time_weighted_level(signalObj, weighting, time)[1]
+    return np.percentile(level, 100 - np.asarray(percents, dtype=float),
+                         axis=0)

tests/test_iec61672.py

@@ -0,0 +1,112 @@
+"""
+Tests for pytta.iec61672, the sound level meter measures of IEC 61672-1:2013
+(frequency weightings A, C and Z, time weightings F and S, Leq, SEL, Lmax,
+Lpeak) and the percentile levels of a time-weighted level. Expected values
+and acceptance limits are those of the standard. 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 iec61672
+
+FS = 48000
+P94 = np.sqrt(2)                       # sine amplitude of 94 dB SPL, in Pa
+
+# Table 3: nominal frequency, A, C, class 1 acceptance limits (+, -)
+TABLE3 = [(10, -70.4, -14.3, 3.0, np.inf), (16, -56.7, -8.5, 2.0, 4.0),
+          (31.5, -39.4, -3.0, 1.5, 1.5), (63, -26.2, -0.8, 1.0, 1.0),
+          (125, -16.1, -0.2, 1.0, 1.0), (250, -8.6, 0.0, 1.0, 1.0),
+          (500, -3.2, 0.0, 1.0, 1.0), (1000, 0.0, 0.0, 0.7, 0.7),
+          (2000, 1.2, -0.2, 1.0, 1.0), (4000, 1.0, -0.8, 1.0, 1.0),
+          (8000, -1.1, -3.0, 1.5, 2.5), (10000, -2.5, -4.4, 2.0, 3.0),
+          (12500, -4.3, -6.2, 2.0, 5.0), (16000, -6.6, -8.5, 2.5, 16.0)]
+
+
+def in_pascal(signalObj):
+    for ch in signalObj.channels.mapping:
+        signalObj.channels[ch].unit = 'Pa'
+    return signalObj
+
+
+def sine(seconds, amplitude, freq=1000.0):
+    t = np.arange(int(round(seconds * FS))) / FS
+    return amplitude * np.sin(2 * np.pi * freq * t)
+
+
+class TestIEC61672(unittest.TestCase):
+
+    def test_frequency_weightings_within_class_1(self):
+        # At the exact frequencies of Table 3 (base-ten, 1000 * 10^(n/10))
+        # the response deviates from the design goals by less than the
+        # class 1 acceptance limits, at 48 kHz.
+        for kind, column in (('A', 1), ('C', 2)):
+            sos = iec61672.weighting_filter(kind, FS)
+            for row in TABLE3:
+                f = 1000 * 10 ** (round(10 * np.log10(row[0] / 1000)) / 10)
+                _, h = ss.sosfreqz(sos, worN=np.array([f]), fs=FS)
+                dev = 20 * np.log10(abs(h[0])) - row[column]
+                self.assertLessEqual(dev, row[3] + 0.05, (kind, row[0]))
+                self.assertGreaterEqual(dev, -row[4] - 0.05, (kind, row[0]))
+
+    def test_toneburst_response_of_time_weighting_f(self):
+        # Table 4: LAFmax - LA of single 4 kHz tonebursts, and LAE - LA,
+        # which is 10 lg(Tb / 1 s).
+        table4 = [(1000, 0.0, 0.5, 0.5), (500, -0.1, 0.5, 0.5),
+                  (200, -1.0, 0.5, 0.5), (100, -2.6, 1.0, 1.0),
+                  (50, -4.8, 1.0, 1.0), (20, -8.3, 1.0, 1.0),
+                  (10, -11.1, 1.0, 1.0), (5, -14.1, 1.0, 1.0),
+                  (2, -18.0, 1.0, 1.5)]
+        steady = pytta.SignalObj(sine(2, P94, 4000), 'time', FS)
+        LA = iec61672.leq(in_pascal(steady))[0]
+        for ms, ref, up, down in table4:
+            burst = np.zeros(2 * FS)
+            burst[FS // 2:FS // 2 + int(ms * FS / 1000)] = \
+                sine(ms / 1000, P94, 4000)
+            sig = in_pascal(pytta.SignalObj(burst, 'time', FS))
+            response = iec61672.lmax(sig, time='F')[0] - LA
+            self.assertLessEqual(response, ref + up, ms)
+            self.assertGreaterEqual(response, ref - down, ms)
+            self.assertAlmostEqual(iec61672.sel(sig)[0] - LA,
+                                   10 * np.log10(ms / 1000), delta=0.05)
+
+    def test_levels_of_a_steady_94_db_tone(self):
+        # 1 Pa rms at 1 kHz: Leq and Lmax 94 dB with A (0 dB at 1 kHz), SEL
+        # of 10 s 104 dB, LCpeak of the amplitude sqrt(2) Pa 97.0 dB.
+        sig = in_pascal(pytta.SignalObj(sine(10, P94), 'time', FS))
+        self.assertAlmostEqual(iec61672.leq(sig)[0], 94.0, delta=0.05)
+        self.assertAlmostEqual(iec61672.lmax(sig, time='S')[0], 94.0,
+                               delta=0.05)
+        self.assertAlmostEqual(iec61672.sel(sig)[0], 104.0, delta=0.05)
+        # The peak is that of the steady tone: a 50 ms fade-in keeps the
+        # onset transient of the C filter out of it.
+        x = sine(2, P94)
+        x[:FS // 20] *= np.sin(np.linspace(0, np.pi / 2, FS // 20)) ** 2
+        faded = in_pascal(pytta.SignalObj(x, 'time', FS))
+        self.assertAlmostEqual(iec61672.lpeak(faded)[0],
+                               94 + 20 * np.log10(np.sqrt(2)), delta=0.05)
+
+    def test_percentile_levels(self):
+        # 30 s at 60 dB then 30 s at 80 dB: L10 is 80 dB and L90 60 dB.
+        x = np.concatenate((sine(30, P94 / 10 ** (34 / 20)),
+                            sine(30, P94 / 10 ** (14 / 20))))
+        sig = in_pascal(pytta.SignalObj(x, 'time', FS))
+        L10, L90 = iec61672.percentile_levels(sig, (10, 90))[:, 0]
+        self.assertAlmostEqual(L10, 80.0, delta=0.1)
+        self.assertAlmostEqual(L90, 60.0, delta=0.1)
+
+    def test_rejects_unknown_weightings(self):
+        sig = pytta.SignalObj(np.zeros(FS), 'time', FS)
+        with self.assertRaises(ValueError):
+            iec61672.leq(sig, weighting='B')
+        with self.assertRaises(ValueError):
+            iec61672.lmax(sig, time='I')
+
+
+if __name__ == '__main__':
+    unittest.main()