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()