Loudness: EBU R128, K-weighted and gated¶

A programme's loudness, K-weighted and gated.

Requires

the programme's Samples, full scale 1.0, one column a channel (a mono programme may be one-dimensional), at its SampleRate (Hz); the channels' weights, 1.0 each unless given (SURROUND_5_0 for L, R, C, Ls, Rs; the LFE is left out).

Provides

Loudness (LUFS, loudness units relative to full scale), or minus infinity for a programme with no block above the absolute gate: silence has no loudness.

Law

each channel K-weighted (a high shelf and the RLB high-pass); over 400 ms blocks every 100 ms, l_j = -0.691 + 10 log10(sum_i G_i z_ij) with z_ij the block's mean square; the blocks above -70 LKFS kept, then those above that set's loudness less 10 LU; the loudness of the kept blocks' mean power.

Source

ITU-R BS.1770-4 (2015), Annex 1 and 2; EBU R128 (2020) and EBU Tech 3341 (2016), whose test signals are the reference.

Proof

reference, limit, invariant.

No Trait form: a tool, not a step law: filters run sample by sample over a whole signal, to measure it; it runs where a game prepares its data, never in a step, so it has no Trait form.

sound.loudness(samples, rate, weights)

The page runs it over EBU Tech 3341's test signals synthesised here. The meter every sample, voice line and music render is measured and normalised by at build, and the bus panel's meter at runtime.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import dmath as D, library
from lusine.kernels.sound.loudness import SURROUND_5_0, TECH_3341, k_weighting, loudness, tone

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import dmath as D, library from lusine.kernels.sound.loudness import SURROUND_5_0, TECH_3341, k_weighting, loudness, tone

The reference¶

The kernel's reference, in numpy: the body under the declaration above, and what every section below runs. It stays as the kernel's test beside its Trait form.

In [2]:

library.reference(loudness)

library.reference(loudness)

Out[2]:

the reference, loudness.py

def loudness(samples, rate, weights=None):
    """INTEGRATED LOUDNESS: a programme's loudness, K-weighted and gated.

    Requires: the programme's `Samples`, full scale 1.0, one column a
    channel (a mono programme may be one-dimensional), at its `SampleRate`
    (Hz); the channels' weights, 1.0 each unless given (SURROUND_5_0 for
    L, R, C, Ls, Rs; the LFE is left out).
    Provides: `Loudness` (LUFS, loudness units relative to full scale), or
    minus infinity for a programme with no block above the absolute gate:
    silence has no loudness.
    Law: each channel K-weighted (a high shelf and the RLB high-pass); over
    400 ms blocks every 100 ms, l_j = -0.691 + 10 log10(sum_i G_i z_ij)
    with z_ij the block's mean square; the blocks above -70 LKFS kept, then
    those above that set's loudness less 10 LU; the loudness of the kept
    blocks' mean power.
    Source: ITU-R BS.1770-4 (2015), Annex 1 and 2; EBU R128 (2020) and
    EBU Tech 3341 (2016), whose test signals are the reference.
    Proof: reference, limit, invariant."""
    power = block_powers(samples, rate, weights)
    if power.size == 0 or not np.any(power > 0):
        return float("-inf")
    lj = np.where(power > 0, OFFSET + 10.0 * _log10(np.where(power > 0, power, 1.0)), -np.inf)
    kept = lj > ABSOLUTE_GATE
    if not kept.any():
        return float("-inf")
    relative = OFFSET + 10.0 * float(_log10(np.asarray(_mean(power[kept])))) + RELATIVE_GATE
    kept &= lj > relative
    return OFFSET + 10.0 * float(_log10(np.asarray(_mean(power[kept]))))

Reference: the standard's filter, and EBU Tech 3341's programmes¶

BS.1770-4 prints the K-weighting's two biquads at 48 kHz; the kernel derives them at any rate from their analogue prototypes, and at 48 kHz they are the printed ones. The filter's gain at 997 Hz is the 0.691 dB the loudness formula takes back out, so a full-scale sine in both channels reads 0 LUFS. The curve is the filter's magnitude response, evaluated from its coefficients on the unit circle.

In [3]:

printed = [(1.53512485958697, -2.69169618940638, 1.19839281085285, -1.69065929318241, 0.73248077421585),
           (1.0, -2.0, 1.0, -1.99004745483398, 0.99007225036621)]
derived = k_weighting(48000)
gap = max(abs(a - b) for s, p in zip(derived, printed) for a, b in zip(s, p))
print(f'derived against printed: {gap:.1e}')
assert gap < 1e-13

def response(f, rate=48000):
    w = 2 * D.PI * np.asarray(f, dtype=float) / rate
    z1 = D.dcos(w) - 1j * D.dsin(w)            # e^-jw
    z2 = D.dcos(2 * w) - 1j * D.dsin(2 * w)
    h = 1.0
    for b0, b1, b2, a1, a2 in k_weighting(rate):
        h = h * (b0 + b1 * z1 + b2 * z2) / (1 + a1 * z1 + a2 * z2)
    return 20 * D.dlog(np.abs(h)) / D.dlog(10.0)

at_997 = float(response(997.0))
print(f'gain at 997 Hz: {at_997:.4f} dB')
assert abs(at_997 - 0.691) < 0.005
f = np.logspace(1, np.log10(20000), 400)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.semilogx(f, response(f))
ax.axvline(997, color='grey', lw=0.5)
ax.set_xlabel('frequency (Hz)')
ax.set_ylabel('K-weighting (dB)')
plt.show()

printed = [, (1.0, -2.0, 1.0)] derived = k_weighting(48000) gap = max(abs(a - b) for s, p in zip(derived, printed) for a, b in zip(s, p)) print(f'derived against printed: {gap.1e}') assert gap < 1e-13 def response(f, rate=48000): w = 2 * D.PI * np.asarray(f, dtype=float) / rate z1 = D.dcos(w) - 1j * D.dsin(w) # e^-jw z2 = D.dcos(2 * w) - 1j * D.dsin(2 * w) h = 1.0 for b0, b1, b2, a1, a2 in k_weighting(rate): h = h * (b0 + b1 * z1 + b2 * z2) / (1 + a1 * z1 + a2 * z2) return 20 * D.dlog(np.abs(h)) / D.dlog(10.0) at_997 = float(response(997.0)) print(f'gain at 997 Hz: {at_997.4f} dB') assert abs(at_997 - 0.691) < 0.005 f = np.logspace(1, np.log10(20000), 400) fig, ax = plt.subplots(figsize=(7, 3.5)) ax.semilogx(f, response(f)) ax.axvline(997, color='grey', lw=0.5) ax.set_xlabel('frequency (Hz)') ax.set_ylabel('K-weighting (dB)') plt.show

derived against printed: 1.8e-15
gain at 997 Hz: 0.6910 dB

EBU Tech 3341's integrated-loudness cases 1 to 5, synthesised at 48 kHz: 1 kHz stereo tones at the levels and for the lengths the specification gives. Each must read within 0.1 LU of its expected loudness. Cases 3 and 4 hold the gates to account: the -36 and -72 dBFS passages must not pull the -23 dBFS programme down.

In [4]:

for i, (segments, want) in enumerate(TECH_3341, 1):
    got = loudness(np.concatenate([tone(level, s) for level, s in segments]), 48000)
    shape = ' + '.join(f'{s:g} s at {level:g} dBFS' for level, s in segments)
    print(f'case {i}: {got:7.2f} LUFS (expected {want:g}): {shape}')
    assert abs(got - want) <= 0.1

for i, (segments, want) in enumerate(TECH_3341, 1): got = loudness(np.concatenate([tone(level, s) for level, s in segments]), 48000) shape = ' + '.join(f'{s:g} s at {level:g} dBFS' for level, s in segments) print(f'case {i}: {got:7.2f} LUFS (expected {want:g}): {shape}') assert abs(got - want) <= 0.1

case 1:  -22.99 LUFS (expected -23): 20 s at -23 dBFS
case 2:  -32.99 LUFS (expected -33): 20 s at -33 dBFS
case 3:  -23.01 LUFS (expected -23): 10 s at -36 dBFS + 60 s at -23 dBFS + 10 s at -36 dBFS
case 4:  -23.01 LUFS (expected -23): 10 s at -72 dBFS + 10 s at -36 dBFS + 60 s at -23 dBFS + 10 s at -36 dBFS + 10 s at -72 dBFS
case 5:  -22.98 LUFS (expected -23): 20 s at -26 dBFS + 20.1 s at -20 dBFS + 20 s at -26 dBFS

Limit: silence has no loudness¶

Silence, a tone under the absolute gate (-70 LKFS), and a programme shorter than one 400 ms block have no block to measure: minus infinity, never a number a normaliser would chase.

In [5]:

for name, x in (('silence', np.zeros((96000, 2))), ('a tone at -80 dBFS', tone(-80.0, 5.0)),
                ('0.3 s of a -23 dBFS tone', tone(-23.0, 0.3))):
    got = loudness(x, 48000)
    print(f'{name}: {got}')
    assert got == float('-inf')

for name, x in (('silence', np.zeros((96000, 2))), ('a tone at -80 dBFS', tone(-80.0, 5.0)), ('0.3 s of a -23 dBFS tone', tone(-23.0, 0.3))): got = loudness(x, 48000) print(f'{name}: {got}') assert got == float('-inf')

silence: -inf
a tone at -80 dBFS: -inf
0.3 s of a -23 dBFS tone: -inf

Invariant: the meter is linear before the log, and blind to the channels' order¶

Doubling a programme's amplitude adds 20 log10 2 = 6.0206 dB exactly, since the filter and the mean square are linear and a doubling is exact in floating point. Swapping the channels changes nothing. A surround channel weighs 1.41, +1.49 dB against a front one (BS.1770-4 Table 3).

In [6]:

x = tone(-30.0, 5.0)
step = loudness(2 * x, 48000) - loudness(x, 48000)
print(f'doubled: +{step:.12f} dB')
assert abs(step - 20 * float(D.dlog(2.0) / D.dlog(10.0))) < 1e-9
a = np.concatenate([tone(-20.0, 3.0, channels=1), tone(-30.0, 3.0, channels=1)], axis=1)
assert loudness(a, 48000) == loudness(a[:, ::-1], 48000)
mono = tone(-25.0, 3.0, channels=1)[:, 0]
silent = np.zeros_like(mono)
front = loudness(np.stack([mono, silent, silent, silent, silent], axis=1), 48000, SURROUND_5_0)
rear = loudness(np.stack([silent, silent, silent, mono, silent], axis=1), 48000, SURROUND_5_0)
print(f'the same signal in a surround channel: +{rear - front:.4f} dB')
assert abs((rear - front) - 10 * float(D.dlog(1.41) / D.dlog(10.0))) < 1e-9

x = tone(-30.0, 5.0) step = loudness(2 * x, 48000) - loudness(x, 48000) print(f'doubled: +{step:.12f} dB') assert abs(step - 20 * float(D.dlog(2.0) / D.dlog(10.0))) < 1e-9 a = np.concatenate([tone(-20.0, 3.0, channels=1), tone(-30.0, 3.0, channels=1)], axis=1) assert loudness(a, 48000) == loudness(a[:, ::-1], 48000) mono = tone(-25.0, 3.0, channels=1)[:, 0] silent = np.zeros_like(mono) front = loudness(np.stack([mono, silent, silent, silent, silent], axis=1), 48000, SURROUND_5_0) rear = loudness(np.stack([silent, silent, silent, mono, silent], axis=1), 48000, SURROUND_5_0) print(f'the same signal in a surround channel: +{rear - front:.4f} dB') assert abs((rear - front) - 10 * float(D.dlog(1.41) / D.dlog(10.0))) < 1e-9

doubled: +6.020599913280 dB
the same signal in a surround channel: +1.4922 dB