[FEATURE] add pytta.bs1770 with loudness, loudness range and true peak

pytta/__init__.py

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

pytta/bs1770.py

@@ -0,0 +1,135 @@
+"""
+PyTTa ITU-R BS.1770 loudness and true peak:
+------------------------------------------
+Loudness of ITU-R BS.1770-5 (K-weighting, gated integrated loudness), the
+momentary and short-term loudness of EBU Tech 3341, the loudness range of
+EBU Tech 3342, and the true-peak level of BS.1770-5 Annex 2.
+
+    >>> pytta.bs1770.loudness(signalObj)       # LUFS (LKFS)
+
+Signals are taken in digital full scale (a sample of 1.0 is 0 dBFS), with
+one channel per loudspeaker of the programme.
+"""
+import numpy as np
+import scipy.signal as ss
+
+# BS.1770-5, Annex 2: order 48, 4-phase FIR interpolating filter
+_PHASES = np.array([
+    [0.0017089843750, -0.0291748046875, -0.0189208984375, -0.0083007812500],
+    [0.0109863281250, 0.0292968750000, 0.0330810546875, 0.0148925781250],
+    [-0.0196533203125, -0.0517578125000, -0.0582275390625, -0.0266113281250],
+    [0.0332031250000, 0.0891113281250, 0.1015625000000, 0.0476074218750],
+    [-0.0594482421875, -0.1665039062500, -0.2003173828125, -0.1022949218750],
+    [0.1373291015625, 0.4650878906250, 0.7797851562500, 0.9721679687500],
+    [0.9721679687500, 0.7797851562500, 0.4650878906250, 0.1373291015625],
+    [-0.1022949218750, -0.2003173828125, -0.1665039062500, -0.0594482421875],
+    [0.0476074218750, 0.1015625000000, 0.0891113281250, 0.0332031250000],
+    [-0.0266113281250, -0.0582275390625, -0.0517578125000, -0.0196533203125],
+    [0.0148925781250, 0.0330810546875, 0.0292968750000, 0.0109863281250],
+    [-0.0083007812500, -0.0189208984375, -0.0291748046875, 0.0017089843750],
+])
+
+
+def k_weighting(fs):
+    """
+    The two biquads of the K-weighting (stage 1 shelf, stage 2 high-pass),
+    as ((b1, a1), (b2, a2)) at the sampling rate fs. The analog prototypes
+    (shelf at 1681.97 Hz, +4.0 dB; high-pass at 38.135 Hz, as also used by
+    libebur128) give Tables 1 and 2 of BS.1770-5 at 48 kHz, and the same
+    response at other rates, as the Recommendation asks.
+    """
+    K = np.tan(np.pi * 1681.974450955533 / fs)
+    Q = 0.7071752369554196
+    Vh = 10 ** (3.999843853973347 / 20)
+    Vb = Vh ** 0.4996667741545416
+    a0 = 1 + K / Q + K ** 2
+    b1 = np.array([Vh + Vb * K / Q + K ** 2, 2 * (K ** 2 - Vh),
+                   Vh - Vb * K / Q + K ** 2]) / a0
+    a1 = np.array([a0, 2 * (K ** 2 - 1), 1 - K / Q + K ** 2]) / a0
+    K = np.tan(np.pi * 38.13547087602444 / fs)
+    Q = 0.5003270373238773
+    a2 = np.array([1 + K / Q + K ** 2, 2 * (K ** 2 - 1),
+                   1 - K / Q + K ** 2]) / (1 + K / Q + K ** 2)
+    return (b1, a1), (np.array([1.0, -2.0, 1.0]), a2)
+
+
+def _block_loudness(signalObj, weights, seconds, hop=0.1):
+    """Loudness of each complete block of `seconds`, every `hop` seconds,
+    and the time of the end of each block."""
+    fs = signalObj.samplingRate
+    x = np.asarray(signalObj.timeSignal, dtype=float)
+    for b, a in k_weighting(fs):
+        x = ss.lfilter(b, a, x, axis=0)
+    weights = np.ones(x.shape[1]) if weights is None else np.asarray(weights)
+    if weights.size != x.shape[1]:
+        raise ValueError('weights must have one value per channel.')
+    length, step = int(round(seconds * fs)), int(round(hop * fs))
+    energy = np.vstack((np.zeros(x.shape[1]), np.cumsum(x ** 2, axis=0)))
+    ends = np.arange(length, x.shape[0] + 1, step)
+    z = (energy[ends] - energy[ends - length]) / length
+    with np.errstate(divide='ignore'):
+        return ends / fs, -0.691 + 10 * np.log10(z @ weights)
+
+
+def _gated(z, threshold):
+    """Mean of the block powers above the threshold, in LUFS."""
+    return -0.691 + 10 * np.log10(np.mean(z[threshold]))
+
+
+def loudness(signalObj, weights=None):
+    """
+    Integrated, gated loudness in LUFS (LKFS) of BS.1770-5: blocks of
+    400 ms with 75 % overlap, an absolute gate at -70 LUFS and a relative
+    gate 10 LU below the absolute-gated loudness.
+
+    * weights (None): weight of each channel (Table 3); 1 for all by
+      default; 1.41 for the surround channels of a 5.1 programme, 0 for LFE.
+    """
+    _, lj = _block_loudness(signalObj, weights, 0.4)
+    power = 10 ** ((lj + 0.691) / 10)
+    absolute = lj > -70
+    relative = _gated(power, absolute) - 10
+    return float(_gated(power, absolute & (lj > relative)))
+
+
+def momentary(signalObj, weights=None):
+    """Momentary loudness of EBU Tech 3341 (400 ms window, no gating),
+    every 100 ms: (time of the end of each window in s, LUFS)."""
+    return _block_loudness(signalObj, weights, 0.4)
+
+
+def short_term(signalObj, weights=None):
+    """Short-term loudness of EBU Tech 3341 (3 s window, no gating), every
+    100 ms: (time of the end of each window in s, LUFS)."""
+    return _block_loudness(signalObj, weights, 3.0)
+
+
+def loudness_range(signalObj, weights=None):
+    """
+    Loudness range in LU of EBU Tech 3342: the difference between the 95th
+    and 10th percentiles of the short-term loudness, after an absolute gate
+    at -70 LUFS and a relative gate 20 LU below the mean of the absolute-
+    gated values. Percentiles as in the MATLAB code of Tech 3342.
+    """
+    _, S = short_term(signalObj, weights)
+    S = S[S >= -70]
+    S = np.sort(S[S >= 10 * np.log10(np.mean(10 ** (S / 10))) - 20])
+    n = S.size
+
+    def pick(percent):              # MATLAB round: halves away from zero
+        return S[int(np.floor((n - 1) * percent / 100 + 1 + 0.5)) - 1]
+    return float(pick(95) - pick(10))
+
+
+def true_peak(signalObj):
+    """
+    True-peak level in dBTP of BS.1770-5 Annex 2: the largest magnitude of
+    the signal oversampled 4 times with the interpolating filter of the
+    Annex, over all channels. Computed in floating point, so the 12.04 dB
+    attenuation of the Annex is not needed.
+    """
+    x = np.asarray(signalObj.timeSignal, dtype=float)
+    peak = np.max(np.abs(x))
+    for phase in _PHASES.T:
+        peak = max(peak, np.max(np.abs(ss.lfilter(phase, 1, x, axis=0))))
+    return float(20 * np.log10(peak))

tests/test_bs1770.py

@@ -0,0 +1,110 @@
+"""
+Tests for pytta.bs1770: loudness and true peak of ITU-R BS.1770-5, with
+the time scales of EBU Tech 3341 and the loudness range of EBU Tech 3342.
+The signals are the minimum requirement test signals of those documents,
+synthesized here, and each test pins their expected response and
+tolerance.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+
+import pytta
+from pytta import bs1770
+
+FS = 48000
+
+
+def tone(segments, channels=2, freq=1000.0, phase=0.0):
+    """1 kHz sine of (seconds, dBFS peak level) segments, on all channels."""
+    parts = []
+    for seconds, level in segments:
+        t = np.arange(int(round(seconds * FS))) / FS
+        amp = 0.0 if level is None else 10 ** (level / 20)
+        parts.append(amp * np.sin(2 * np.pi * freq * t + phase))
+    return np.tile(np.concatenate(parts)[:, None], (1, channels))
+
+
+class TestBS1770(unittest.TestCase):
+
+    def test_k_weighting_at_48_khz_is_that_of_tables_1_and_2(self):
+        (b1, a1), (b2, a2) = bs1770.k_weighting(48000)
+        self.assertTrue(np.allclose(b1, [1.53512485958697, -2.69169618940638,
+                                         1.19839281085285], atol=1e-10))
+        self.assertTrue(np.allclose(a1, [1, -1.69065929318241,
+                                         0.73248077421585], atol=1e-10))
+        self.assertTrue(np.allclose(b2, [1.0, -2.0, 1.0]))
+        self.assertTrue(np.allclose(a2, [1, -1.99004745483398,
+                                         0.99007225036621], atol=1e-10))
+
+    def test_steady_tones_tech_3341_cases_1_2_and_6(self):
+        # M, S, I = -23.0 and -33.0 LUFS for stereo tones; the 5.0 case
+        # pins the weights of Table 3 (1.41 on the surrounds).
+        for level in (-23.0, -33.0):
+            sig = pytta.SignalObj(tone([(20, level)]), 'time', FS)
+            I = bs1770.loudness(sig)
+            _, M = bs1770.momentary(sig)
+            _, S = bs1770.short_term(sig)
+            print(f'\n  {level} dBFS: I {I:.2f}, M {np.median(M):.2f}, '
+                  f'S {np.median(S):.2f}')
+            for value in (I, np.median(M), np.median(S)):
+                self.assertAlmostEqual(value, level, delta=0.1)
+        levels = [-28.0, -28.0, -24.0, -30.0, -30.0]
+        sig = pytta.SignalObj(np.column_stack(
+            [tone([(20, L)], 1)[:, 0] for L in levels]), 'time', FS)
+        I = bs1770.loudness(sig, weights=[1, 1, 1, 1.41, 1.41])
+        self.assertAlmostEqual(I, -23.0, delta=0.1)
+
+    def test_gating_tech_3341_cases_3_4_and_5(self):
+        cases = ([(10, -36), (60, -23), (10, -36)],
+                 [(10, -72), (10, -36), (60, -23), (10, -36), (10, -72)],
+                 [(20, -26), (20.1, -20), (20, -26)])
+        for segments in cases:
+            sig = pytta.SignalObj(tone(segments), 'time', FS)
+            I = bs1770.loudness(sig)
+            print(f'\n  {segments}: I {I:.2f}')
+            self.assertAlmostEqual(I, -23.0, delta=0.1)
+
+    def test_window_lengths_tech_3341_cases_9_and_12(self):
+        # Alternating tones whose energy averages to -23 LUFS over exactly
+        # 3 s (S) and 400 ms (M): constant after the first window.
+        sig = pytta.SignalObj(tone([(1.34, -20), (1.66, -30)] * 5), 'time', FS)
+        t, S = bs1770.short_term(sig)
+        self.assertTrue(np.allclose(S[t >= 3.0], -23.0, atol=0.1))
+        sig = pytta.SignalObj(tone([(0.18, -20), (0.22, -30)] * 25), 'time', FS)
+        t, M = bs1770.momentary(sig)
+        self.assertTrue(np.allclose(M[t >= 1.0], -23.0, atol=0.1))
+
+    def test_loudness_range_tech_3342_cases_1_to_4(self):
+        cases = (([(20, -20), (20, -30)], 10), ([(20, -20), (20, -15)], 5),
+                 ([(20, -40), (20, -20)], 20),
+                 ([(20, -50), (20, -35), (20, -20), (20, -35), (20, -50)], 15))
+        for segments, expected in cases:
+            sig = pytta.SignalObj(tone(segments + [(1.5, None)]), 'time', FS)
+            lra = bs1770.loudness_range(sig)
+            print(f'\n  LRA {lra:.2f} (expected {expected} +-1)')
+            self.assertAlmostEqual(lra, expected, delta=1)
+
+    def test_true_peak_tech_3341_cases_15_to_19(self):
+        # Sines between samples: the true peak of an amplitude 0.5 is
+        # -6.0 dBTP and of 1.41 is +3.0, within +0.2 / -0.4 dB.
+        cases = ((FS / 4, 0.50, 0.0, -6.0), (FS / 4, 0.50, 45.0, -6.0),
+                 (FS / 6, 0.50, 60.0, -6.0), (FS / 8, 0.50, 67.5, -6.0),
+                 (FS / 4, 1.41, 45.0, 3.0))
+        t = np.arange(FS) / FS
+        fade = np.minimum(1, np.minimum(t, t[::-1]) / 0.01)
+        for freq, amp, phase, expected in cases:
+            x = amp * fade * np.sin(2 * np.pi * freq * t + np.radians(phase))
+            sig = pytta.SignalObj(np.column_stack((x, x)), 'time', FS)
+            tp = bs1770.true_peak(sig)
+            print(f'\n  {freq:.0f} Hz {phase} deg: {tp:.2f} dBTP '
+                  f'(expected {expected})')
+            self.assertGreaterEqual(tp, expected - 0.4)
+            self.assertLessEqual(tp, expected + 0.2)
+
+
+if __name__ == '__main__':
+    unittest.main()