diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..c9a72ff 100644 --- a/pytta/__init__.py +++ b/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 diff --git a/pytta/bs1770.py b/pytta/bs1770.py new file mode 100644 index 0000000..b6d82a3 --- /dev/null +++ b/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)) diff --git a/tests/test_bs1770.py b/tests/test_bs1770.py new file mode 100644 index 0000000..ea937e4 --- /dev/null +++ b/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()