From 8f742bbb59366fc337c85ef20aa74ceb5650cc22 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 20:42:53 +0100 Subject: [PATCH] [FEATURE] add pytta.iec60268 with the speech transmission index pytta.iec60268.sti gives the STI of IEC 60268-16:2011 from an impulse response by the indirect method of clause 6: the modulation transfer function from the Schroeder equation in the seven octave bands and 14 modulation frequencies, the noise term when background levels are given, the auditory masking of Table A.1 and the reception threshold of Table A.2 when speech levels are given, and the male or female weights of Table A.3. It returns the STI, the MTI per band and the MTF. Tests pin a perfect channel at STI 1, the masking at 100 dB, noise at 0 dB SNR, and the MTF of an exponential decay, within the 0.01 STI that the standard allows for the octave filters. Co-Authored-By: Claude Fable 5.1 --- pytta/__init__.py | 2 ++ pytta/iec60268.py | 89 +++++++++++++++++++++++++++++++++++++++++++++ tests/test_sti.py | 92 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 183 insertions(+) create mode 100644 pytta/iec60268.py create mode 100644 tests/test_sti.py diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..62c6ef8 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -108,6 +108,7 @@ from . import utils from . import rooms from . import iso3741 +from . import iec60268 from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +131,7 @@ 'generate', 'utils', 'iso3741', + 'iec60268', 'rooms', # Functions diff --git a/pytta/iec60268.py b/pytta/iec60268.py new file mode 100644 index 0000000..411d11f --- /dev/null +++ b/pytta/iec60268.py @@ -0,0 +1,89 @@ +""" +PyTTa IEC 60268-16 speech transmission index: +-------------------------------------------- +Speech transmission index (STI) of IEC 60268-16:2011 from an impulse +response, by the indirect method of clause 6. + + >>> sti, mti, mtf = pytta.iec60268.sti(ir, levels, noise) + +""" +import numpy as np + +from pytta.classes import OctFilter + +BANDS = (125, 250, 500, 1000, 2000, 4000, 8000) +MODULATION = np.array([0.63, 0.8, 1, 1.25, 1.6, 2, 2.5, 3.15, 4, 5, 6.3, 8, + 10, 12.5]) # A.2.2 +_ART = np.array([46, 27, 12, 6.5, 7.5, 8, 12]) # Table A.2 +_WEIGHTS = { # Table A.3 + False: (np.array([0.085, 0.127, 0.230, 0.233, 0.309, 0.224, 0.173]), + np.array([0.085, 0.078, 0.065, 0.011, 0.047, 0.095])), + True: (np.array([0, 0.117, 0.223, 0.216, 0.328, 0.250, 0.194]), + np.array([0, 0.099, 0.066, 0.062, 0.025, 0.076])), +} + + +def _masking_db(L): + """Auditory masking in dB from the level of the band below (Table A.1).""" + return np.select([L < 63, L < 67, L < 100], + [0.5 * L - 65, 1.8 * L - 146.9, 0.5 * L - 59.8], -10) + + +def sti(ir, levels=None, noise=None, female=False): + """ + Speech transmission index from an impulse response (IEC 60268-16:2011, + clause 6 and Annex A). + + Parameters (default), (type): + ----------------------------- + * ir (), (SignalObj): impulse response, one channel, at least 1.6 s + long, with a sampling rate above 16 kHz; + * levels (None), (7 values): speech level at the listener in the + octave bands from 125 Hz to 8 kHz, in dB SPL. With None the + auditory masking and the reception threshold (A.3) are left out, + as for signals not measured in dB SPL; + * noise (None), (7 values): background noise in the same bands, in + dB SPL; needs levels; + * female (False), (bool): female weighting of Table A.3. + + Return: + ------- + * STI, limited to 1; + * MTI: modulation transfer index of each octave band; + * MTF: modulation transfer ratio, 14 modulation frequencies (rows) + by 7 bands (columns), before the masking correction. + + The MTF comes from the Schroeder equation of 6.1 on each octave band, + times 1 / (1 + 10**(-SNR / 10)) when noise is given. Masking (Table A.1, + from the total level of the band below) and the reception threshold + correct it as in A.5.3, and A.5.4 to A.5.6 give the STI. + """ + for name, value in (('levels', levels), ('noise', noise)): + if value is not None and len(value) != 7: + raise ValueError(f'{name} must have 7 octave band values.') + if noise is not None and levels is None: + raise ValueError('noise needs the speech levels.') + fs = ir.samplingRate + bands = OctFilter(order=4, nthOct=1, samplingRate=fs, minFreq=125, + maxFreq=8000, refFreq=1000, base=10) + h2 = bands.filter(ir)[0].timeSignal.astype(float) ** 2 + t = np.arange(h2.shape[0]) / fs + mtf = np.abs(np.exp(-2j * np.pi * MODULATION[:, None] * t) @ h2) \ + / h2.sum(axis=0) + m = mtf.copy() + if levels is not None: + signal = 10 ** (np.asarray(levels, dtype=float) / 10) + total = signal.copy() + if noise is not None: + noiseI = 10 ** (np.asarray(noise, dtype=float) / 10) + m = m / (1 + noiseI / signal) + total = total + noiseI + masking = np.r_[0, total[:-1] * 10 ** (_masking_db( + 10 * np.log10(total[:-1])) / 10)] + m = m * signal / (signal + masking + 10 ** (_ART / 10)) + m = np.clip(m, 1e-12, 1 - 1e-12) + snr = np.clip(10 * np.log10(m / (1 - m)), -15, 15) + mti = np.mean((snr + 15) / 30, axis=0) + alpha, beta = _WEIGHTS[female] + value = np.sum(alpha * mti) - np.sum(beta * np.sqrt(mti[:-1] * mti[1:])) + return min(float(value), 1.0), mti, mtf diff --git a/tests/test_sti.py b/tests/test_sti.py new file mode 100644 index 0000000..0881a54 --- /dev/null +++ b/tests/test_sti.py @@ -0,0 +1,92 @@ +""" +Tests for pytta.iec60268.sti, the speech transmission index of +IEC 60268-16:2011 by the indirect method (clause 6). Expected values are +worked out here from the equations of the standard; the octave filters may +add a systematic error up to 0.01 STI (A.5.1.2), which is the tolerance. +Each test pins one rule. +""" +import os +import unittest + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import numpy as np + +import pytta + +FS = 48000 +TOL = 0.01 # A.5.1.2 +ALPHA = np.array([0.085, 0.127, 0.230, 0.233, 0.309, 0.224, 0.173]) +BETA = np.array([0.085, 0.078, 0.065, 0.011, 0.047, 0.095]) +FM = np.array([0.63, 0.8, 1, 1.25, 1.6, 2, 2.5, 3.15, 4, 5, 6.3, 8, 10, + 12.5]) + + +def perfect_channel(): + # no local variables: SignalObj inspects the locals of its caller + return pytta.SignalObj(np.eye(1, 2 * FS, 100)[0], 'time', FS) + + +def sti_from_m(m): + """STI from one modulation transfer ratio per band (A.5.4 to A.5.6).""" + snr = np.clip(10 * np.log10(m / (1 - m)), -15, 15) + mti = (snr + 15) / 30 + return np.sum(ALPHA * mti) - np.sum(BETA * np.sqrt(mti[:-1] * mti[1:])) + + +class TestSTI(unittest.TestCase): + + def test_perfect_channel_gives_one(self): + # The weights and redundancies of Table A.3 add up to 1. + for female in (False, True): + sti, mti, mtf = pytta.iec60268.sti(perfect_channel(), + female=female) + self.assertGreater(sti, 1 - TOL) + self.assertLessEqual(sti, 1) + + def test_masking_and_reception_threshold(self): + # 100 dB in every band: masking of -10 dB from the band below + # (Table A.1), so m = 1 / 1.1 from 250 Hz up; the threshold of + # Table A.2 is 54 dB down or more. + sti, mti, mtf = pytta.iec60268.sti(perfect_channel(), + levels=[100] * 7) + m = np.r_[1 / (1 + 10 ** ((46 - 100) / 10)), np.full(6, 1 / 1.1)] + expected = sti_from_m(np.minimum(m, 1 - 1e-12)) + print(f'\n STI {sti:.4f} (expected {expected:.4f})') + self.assertAlmostEqual(sti, expected, delta=TOL) + + def test_noise_at_0_db_snr(self): + # Noise as loud as the signal halves m (6.1), and the masking + # comes from the total level, 103 dB: m = 0.5 / 1.2 from 250 Hz up. + sti, mti, mtf = pytta.iec60268.sti(perfect_channel(), + levels=[100] * 7, noise=[100] * 7) + m = np.r_[0.5, np.full(6, 0.5 / 1.2)] + expected = sti_from_m(m) + print(f'\n STI {sti:.4f} (expected {expected:.4f})') + self.assertAlmostEqual(sti, expected, delta=TOL) + + def test_mtf_of_an_exponential_decay(self): + # Diffuse decay with T = 1 s: m(F) = 1 / sqrt(1 + (2 pi F T / 13.8)^2) + # (Houtgast and Steeneken). A noise decay only follows it on + # average: its squared envelope fluctuates by about + # 1 / sqrt(bandwidth x decay time), which in the 125 Hz band makes + # m wrong by 0.2 or more. The bands from 2 kHz up stay within 0.07 + # (largest of 8 noise seeds). + rng = np.random.default_rng(0) + t = np.arange(2 * FS) / FS + h = rng.standard_normal(t.size) * 10 ** (-3 * t) + _, _, mtf = pytta.iec60268.sti(pytta.SignalObj(h, 'time', FS)) + theory = 1 / np.sqrt(1 + (2 * np.pi * FM * 1.0 / 13.8) ** 2) + worst = np.max(np.abs(mtf[:, 4:] - theory[:, None])) + print(f'\n largest |m - theory| = {worst:.3f}') + self.assertLess(worst, 0.07) + + def test_rejects_levels_of_other_lengths(self): + with self.assertRaises(ValueError): + pytta.iec60268.sti(perfect_channel(), levels=[60] * 6) + with self.assertRaises(ValueError): + pytta.iec60268.sti(perfect_channel(), noise=[40] * 7) + + +if __name__ == '__main__': + unittest.main()