[FEATURE] add impact sound insulation (ISO 717-2) to pytta.iso717

pytta/__init__.py

@@ -108,6 +108,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import iso717
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +131,7 @@
            'generate',
            'utils',
            'iso3741',
+           'iso717',
            'rooms',
 
            # Functions

pytta/iso717.py

@@ -0,0 +1,137 @@
+"""
+PyTTa ISO 717 rating:
+--------------------
+Single-number quantities of airborne sound insulation, ISO 717-1:2013, and
+of impact sound insulation, ISO 717-2.
+
+    >>> Rw, C, Ctr = pytta.iso717.airborne(R)
+    >>> Lnw, CI = pytta.iso717.impact(Ln)
+
+"""
+import numpy as np
+
+# ISO 717-1:2013, Table 3 (reference values) and Table 4 (spectra No. 1 and
+# No. 2 for C and Ctr), from 100 Hz to 3 150 Hz in thirds and from 125 Hz
+# to 2 000 Hz in octaves.
+_TABLES = {
+    16: dict(reference=[33, 36, 39, 42, 45, 48, 51, 52, 53, 54, 55, 56, 56,
+                        56, 56, 56],
+             spectra=([-29, -26, -23, -21, -19, -17, -15, -13, -12, -11, -10,
+                       -9, -9, -9, -9, -9],
+                      [-20, -20, -18, -16, -15, -14, -13, -12, -11, -9, -8,
+                       -9, -10, -11, -13, -15]),
+             limit=320, at500=7),
+    5: dict(reference=[36, 45, 52, 55, 56],
+            spectra=([-21, -14, -8, -5, -4], [-14, -10, -7, -4, -6]),
+            limit=100, at500=2),
+}
+
+
+def _half_up(x, step):
+    """Round to a multiple of step, halves up (ISO 80000-1), clearing the
+    binary error of the input first."""
+    return np.floor(np.round(np.asarray(x, dtype=float) / step, 6) + 0.5)
+
+
+def _shift(tenths, reference, limit):
+    """Largest shift s, in dB, of the reference curve (in tenths of a dB)
+    such that the deviations of reference + s above tenths sum to no more
+    than limit."""
+    shift = int(np.floor((tenths - reference).min() / 10))  # no deviation
+    while np.maximum(reference + (shift + 1) * 10 - tenths, 0).sum() \
+            <= limit:
+        shift += 1
+    return shift
+
+
+def airborne(R):
+    """
+    Weighted sound reduction index and spectrum adaptation terms of
+    ISO 717-1:2013, from R, R', Dn or DnT by band.
+
+    Parameters:
+    -----------
+        * R (array-like):
+            16 one-third-octave bands from 100 Hz to 3 150 Hz or 5 octave
+            bands from 125 Hz to 2 000 Hz, in dB.
+
+    Return:
+    -------
+        * (Xw, C, Ctr), integers in dB: Rw (or R'w, Dn,w, DnT,w) and the
+          adaptation terms for pink noise and for urban traffic noise.
+
+    R is rounded to 0.1 dB, halves up (4.4). The reference curve is shifted
+    in steps of 1 dB until the sum of unfavourable deviations is as large as
+    possible but not more than 32.0 dB in thirds or 10.0 dB in octaves, and
+    Xw is its value at 500 Hz. C and Ctr follow Equations (1) and (2).
+    Values are kept in tenths of a decibel as integers, so the limit is
+    compared exactly.
+    """
+    tables = _TABLES.get(len(R))
+    if tables is None:
+        raise ValueError('R must have 16 one-third-octave or 5 octave bands.')
+    tenths = _half_up(R, 0.1)
+    reference = np.array(tables['reference']) * 10
+    shift = _shift(tenths, reference, tables['limit'])
+    Xw = tables['reference'][tables['at500']] + shift
+    C = [int(_half_up(-10 * np.log10(np.sum(10 ** ((np.array(L) - tenths / 10)
+                                                    / 10))), 1)) - Xw
+         for L in tables['spectra']]
+    return Xw, C[0], C[1]
+
+
+# ISO 717-2 reference values for impact sound, from 100 Hz to 3 150 Hz in
+# thirds and from 125 Hz to 2 000 Hz in octaves, as restated by Crocker
+# (ed.), Handbook of Noise and Vibration Control (2007), chapter 106, Table 13
+# and Figure 4. The octave rating is 5 dB below the curve at 500 Hz, and CI
+# sums the bands up to 2 500 Hz (the first 15 thirds, all 5 octaves).
+_IMPACT = {
+    16: dict(reference=[62, 62, 62, 62, 62, 62, 61, 60, 59, 58, 57, 54, 51,
+                        48, 45, 42],
+             limit=320, at500=7, offset=0, upto2500=15),
+    5: dict(reference=[67, 67, 65, 62, 49],
+            limit=100, at500=2, offset=-5, upto2500=5),
+}
+
+
+def impact(Ln):
+    """
+    Weighted normalized impact sound pressure level and spectrum adaptation
+    term of ISO 717-2, from Ln, L'n or L'nT by band.
+
+    Parameters:
+    -----------
+        * Ln (array-like):
+            16 one-third-octave bands from 100 Hz to 3 150 Hz, 5 octave
+            bands from 125 Hz to 2 000 Hz, or 19 one-third-octave bands from
+            50 Hz to 3 150 Hz, in dB.
+
+    Return:
+    -------
+        * (Lnw, CI), integers in dB: Ln,w (or L'n,w, L'nT,w) and the
+          adaptation term for walking noise; (Lnw, CI, CI50) with the
+          19 bands, CI50 being CI,50-2500.
+
+    Ln is rounded to 0.1 dB, halves up. An unfavourable deviation is Ln above
+    the reference curve. The curve is shifted in steps of 1 dB until the sum
+    of unfavourable deviations is as large as possible but not more than
+    32.0 dB in thirds or 10.0 dB in octaves; Ln,w is its value at 500 Hz,
+    less 5 dB in octaves. CI = Ln,sum - 15 - Ln,w, with Ln,sum the energetic
+    sum of the bands from 100 Hz (125 Hz in octaves) to 2 500 Hz;
+    CI,50-2500 sums from 50 Hz. Ln,w is rated on 100 Hz to 3 150 Hz only.
+    """
+    tenths = _half_up(Ln, 0.1)
+    extended = len(tenths) == 19  # 50, 63 and 80 Hz come first
+    low, tenths = (tenths[:3], tenths[3:]) if extended else ([], tenths)
+    tables = _IMPACT.get(len(tenths))
+    if tables is None:
+        raise ValueError('Ln must have 16 or 19 one-third-octave or 5 octave '
+                         'bands.')
+    reference = np.array(tables['reference']) * 10
+    # an impact deviation is the airborne one with both curves negated
+    shift = -_shift(-tenths, -reference, tables['limit'])
+    Lnw = tables['reference'][tables['at500']] + shift + tables['offset']
+    bands = tenths[:tables['upto2500']]
+    ranges = (bands, np.r_[low, bands]) if extended else (bands,)
+    return (Lnw, *[int(_half_up(10 * np.log10(np.sum(10 ** (L / 100))), 1))
+                   - 15 - Lnw for L in ranges])

tests/test_iso717.py

@@ -0,0 +1,50 @@
+"""
+Tests for pytta.iso717.airborne, the single-number rating of airborne
+sound insulation of ISO 717-1:2013. Each test pins one rule.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import pytta
+
+# ISO 717-1:2013, Annex C, Table C.1: sound reduction index of a building
+# element, one-third-octave bands from 100 Hz to 3 150 Hz.
+R_ANNEX_C = [20.4, 16.3, 17.7, 22.6, 22.4, 22.7, 24.8, 26.6, 28.0, 30.5, 31.8,
+             32.5, 33.4, 33.0, 31.0, 25.5]
+
+
+class TestISO717(unittest.TestCase):
+
+    def test_example_of_annex_c(self):
+        # Rw(C;Ctr) = 30 (-2; -3) dB, with 31.8 dB of unfavourable deviation
+        self.assertEqual(pytta.iso717.airborne(R_ANNEX_C), (30, -2, -3))
+
+    def test_deviation_limit_is_inclusive(self):
+        # Octave bands: R 5 dB above the reference curve. Shifting it up by
+        # 7 dB gives 2 dB below R in each of the 5 bands, exactly 10.0 dB,
+        # which is allowed; 8 dB would give 15 dB.
+        reference = [36, 45, 52, 55, 56]
+        Rw, C, Ctr = pytta.iso717.airborne([r + 5 for r in reference])
+        self.assertEqual(Rw, 52 + 7)
+
+    def test_input_is_rounded_half_up_to_one_decimal(self):
+        # The standard rounds R to 0.1 dB before rating, halves up. These
+        # five bands end in 5, and round() of Python, which rounds the
+        # binary value, takes each of them 0.1 dB down: the deviations would
+        # sum 32.3 dB and Rw would drop to 29. Rounded half up, the curve is
+        # that of Annex C.
+        R = list(R_ANNEX_C)
+        for band, value in ((5, 22.65), (8, 27.95), (9, 30.45), (14, 30.95),
+                            (15, 25.45)):
+            R[band] = value
+        self.assertEqual(pytta.iso717.airborne(R), (30, -2, -3))
+
+    def test_rejects_other_band_counts(self):
+        with self.assertRaises(ValueError):
+            pytta.iso717.airborne(R_ANNEX_C[:-1])
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_iso717_2.py

@@ -0,0 +1,100 @@
+"""
+Tests for pytta.iso717.impact, the single-number rating of impact sound
+insulation of ISO 717-2. Each test pins one rule.
+
+The text of ISO 717-2 is not at hand; the worked example is the one of
+Crocker (ed.), Handbook of Noise and Vibration Control (2007), chapter 106,
+Table 13, which follows ISO 717-2.
+"""
+import os
+import unittest
+
+import numpy as np
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import pytta
+
+REFERENCE_THIRDS = [62, 62, 62, 62, 62, 62, 61, 60, 59, 58, 57, 54, 51, 48,
+                    45, 42]
+REFERENCE_OCTAVES = [67, 67, 65, 62, 49]
+
+
+def unfavourable(Ln, reference, Lnw, at500):
+    """Sum of the deviations of Ln above the reference curve shifted to Lnw
+    at 500 Hz."""
+    shifted = np.array(reference) - reference[at500] + Lnw
+    return np.maximum(np.array(Ln) - shifted, 0).sum()
+
+
+class TestISO717Impact(unittest.TestCase):
+
+    def test_example_of_crocker_table_13(self):
+        # Bare floor: Ln,w(CI) = 77 (-11) dB, 29.4 dB of unfavourable
+        # deviation. Floor with covering: Ln,w = 40 dB with 29.2 dB; the
+        # table prints CI = -3 dB, but its own Ln,sum = 57.9 dB gives
+        # 57.9 - 15 - 40 = +2.9 dB, and the floor covering term of the same
+        # table (-11 - Cr = -14, Cr = +3) confirms the sign is a misprint.
+        bare = [65.5, 66.5, 67, 67.4, 67.8, 68.2, 68.6, 69, 69.3, 69.6, 70,
+                70.3, 70.5, 70.7, 71, 71.2]
+        covered = [56.2, 50, 47.4, 43.6, 39.1, 36.4, 33.2, 31.7, 31, 29.6,
+                   24.7, 23.5, 22.6, 24.3, 19.1, 14.3]
+        for Ln, expected, deviation in ((bare, (77, -11), 29.4),
+                                        (covered, (40, 3), 29.2)):
+            result = pytta.iso717.impact(Ln)
+            sum_ = unfavourable(Ln, REFERENCE_THIRDS, result[0], 7)
+            print(f'\n{result} vs {expected}; deviations {sum_:.1f} dB vs '
+                  f'{deviation} dB')
+            self.assertEqual(result, expected)
+            self.assertAlmostEqual(sum_, deviation, places=6)
+
+    def test_octave_bands_rate_5_dB_below_the_curve_at_500_Hz(self):
+        # Ln 5 dB above the octave reference curve. Shifting it up by 3 dB
+        # leaves 2 dB in each of the 5 bands, exactly 10.0 dB, which is
+        # allowed; 2 dB would give 15 dB. The curve then reads 65 + 3 dB at
+        # 500 Hz and Ln,w is 5 dB less.
+        Ln = [r + 5 for r in REFERENCE_OCTAVES]
+        Lnw, CI = pytta.iso717.impact(Ln)
+        print(f'\nLn,w = {Lnw} dB, deviations '
+              f'{unfavourable(Ln, REFERENCE_OCTAVES, Lnw + 5, 2):.1f} dB')
+        self.assertEqual(Lnw, 63)
+
+    def test_octave_rating_within_1_dB_of_the_third_octave_rating(self):
+        # Crocker, chapter 106: the octave reference curve keeps the ratings
+        # from thirds and from octaves within 1 dB. The floors of Table 13,
+        # summed to octaves, rate 77 and 41 dB against 77 and 40 dB in
+        # thirds. This pins the octave curve, 49 dB at 2 kHz included: with
+        # 54 dB there the bare floor would rate 73 dB.
+        floors = ([65.5, 66.5, 67, 67.4, 67.8, 68.2, 68.6, 69, 69.3, 69.6, 70,
+                   70.3, 70.5, 70.7, 71, 71.2],
+                  [56.2, 50, 47.4, 43.6, 39.1, 36.4, 33.2, 31.7, 31, 29.6,
+                   24.7, 23.5, 22.6, 24.3, 19.1, 14.3])
+        for Ln in floors:
+            octaves = [round(10 * np.log10(np.sum(10**(np.array(Ln[i:i + 3])
+                                                      / 10))), 1)
+                       for i in range(0, 15, 3)]
+            thirds, octave = (pytta.iso717.impact(Ln)[0],
+                              pytta.iso717.impact(octaves)[0])
+            print(f'\nLn,w {thirds} dB in thirds, {octave} dB in octaves')
+            self.assertLessEqual(abs(octave - thirds), 1)
+
+    def test_adaptation_terms_sum_up_to_2500_Hz(self):
+        # Flat 59.6 dB from 50 Hz to 3 150 Hz (19 thirds). Shifting the curve
+        # to 66 dB leaves 28.4 dB of deviation (65 dB would give 33.0 dB).
+        # CI sums the 15 bands from 100 Hz to 2 500 Hz:
+        # 59.6 + 10 lg 15 - 15 - 66 = -9.64, so -10; with the 3 150 Hz band
+        # it would be -9.36, so -9. CI,50-2500 adds 50, 63 and 80 Hz:
+        # 59.6 + 10 lg 18 - 15 - 66 = -8.85, so -9.
+        result = pytta.iso717.impact([59.6] * 19)
+        CI = 59.6 + 10 * np.log10(15) - 15 - 66
+        CI50 = 59.6 + 10 * np.log10(18) - 15 - 66
+        print(f'\n{result}; CI {CI:.2f}, CI,50-2500 {CI50:.2f}')
+        self.assertEqual(result, (66, -10, -9))
+
+    def test_rejects_other_band_counts(self):
+        with self.assertRaises(ValueError):
+            pytta.iso717.impact([60.0] * 15)
+
+
+if __name__ == '__main__':
+    unittest.main()