[FEATURE] add pytta.iec60268 with the speech transmission index

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

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

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