[FEATURE] add the room and specimen limits of ISO 354 to pytta.iso354

pytta/__init__.py

@@ -108,6 +108,8 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import iso9613
+from . import iso354
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +132,8 @@
            'generate',
            'utils',
            'iso3741',
+           'iso9613',
+           'iso354',
            'rooms',
 
            # Functions

pytta/iso354.py

@@ -0,0 +1,93 @@
+"""
+PyTTa ISO 354 sound absorption in a reverberation room:
+------------------------------------------------------
+Equivalent sound absorption area and sound absorption coefficient of
+ISO 354:2003 (Equations 5 to 9), with the power attenuation of air from
+ISO 9613-1, the Sabine and Eyring reverberation formulas, and the limits
+of 6.1.4 and 6.2.1.1 on the room and the specimen.
+
+    >>> alpha = pytta.iso354.absorption_coefficient(T1, T2, V, S, t1, t2)
+
+"""
+import numpy as np
+
+from pytta import iso9613
+
+
+def speed_of_sound(temperature):
+    """Speed of sound in m/s, Equation (6): 331 + 0.6 t, for t from 15 to
+    30 degrees Celsius."""
+    return 331 + 0.6 * temperature
+
+
+def power_attenuation(f, temperature=20., humidity=50., pressure=101.325):
+    """Power attenuation coefficient m of air, in 1/m, from the attenuation
+    alpha of ISO 9613-1 (dB/m): m = alpha / (10 lg e)."""
+    return iso9613.attenuation(f, temperature, humidity, pressure) \
+        / (10 * np.log10(np.e))
+
+
+def absorption_area(T, V, c, m=0.):
+    """Equivalent sound absorption area in m2 of a room of volume V (m3)
+    with reverberation time T (s), Equations (5) and (7):
+    A = 55.3 V / (c T) - 4 V m."""
+    return 55.3 * V / (c * np.asarray(T)) - 4 * V * m
+
+
+def absorption_coefficient(T1, T2, V, S, t1=20., t2=None, m1=0., m2=None):
+    """
+    Sound absorption coefficient of a test specimen of area S (m2),
+    Equations (8) and (9), from the reverberation times of the empty room
+    (T1) and with the specimen (T2), with the air temperatures t1 and t2
+    (degrees Celsius) and power attenuations m1 and m2 of each measurement.
+    t2 and m2 default to t1 and m1.
+    """
+    t2 = t1 if t2 is None else t2
+    m2 = m1 if m2 is None else m2
+    A1 = absorption_area(T1, V, speed_of_sound(t1), m1)
+    A2 = absorption_area(T2, V, speed_of_sound(t2), m2)
+    return (A2 - A1) / S
+
+
+def sabine(V, A, c, m=0.):
+    """Reverberation time in s by Sabine, the inverse of Equation (5):
+    T = 55.3 V / (c (A + 4 V m))."""
+    return 55.3 * V / (c * (np.asarray(A) + 4 * V * m))
+
+
+def eyring(V, S, alpha, c, m=0.):
+    """Reverberation time in s by Eyring, for room surface S (m2) with mean
+    absorption coefficient alpha: T = 55.3 V / (c (-S ln(1 - alpha) +
+    4 V m))."""
+    return 55.3 * V / (c * (-S * np.log(1 - np.asarray(alpha)) + 4 * V * m))
+
+
+BANDS = np.array([100, 125, 160, 200, 250, 315, 400, 500, 630, 800, 1000,
+                  1250, 1600, 2000, 2500, 3150, 4000, 5000])
+_MAX_AREA_200 = np.array([6.5] * 10 + [7.0, 7.5, 8.0, 9.5, 10.5, 12.0, 13.0,
+                                       14.0])
+
+
+def max_room_absorption(V):
+    """Largest equivalent sound absorption area in m2 that the empty room of
+    volume V (m3) may have in each third octave band of BANDS: Table 1,
+    for 200 m3, times (V / 200)^(2/3) (6.1.4)."""
+    return _MAX_AREA_200 * (V / 200) ** (2 / 3)
+
+
+def specimen_area(V):
+    """Range of the area in m2 of a plane test specimen: 10 to 12 m2, the
+    upper limit times (V / 200)^(2/3) for V above 200 m3 (6.2.1.1)."""
+    return 10., 12. * max(1., (V / 200) ** (2 / 3))
+
+
+def room_qualifies(A, V):
+    """Whether the absorption area A of the empty room, one value per band of
+    BANDS, meets 6.1.4: not above Table 1 scaled to V, and no band more
+    than 15 % from the mean of its two neighbours. Returns one bool per
+    band (the outer bands are checked against Table 1 only)."""
+    A = np.asarray(A, dtype=float)
+    ok = A <= max_room_absorption(V)
+    mean = (A[:-2] + A[2:]) / 2
+    ok[1:-1] &= np.abs(A[1:-1] - mean) <= 0.15 * mean
+    return ok

pytta/iso9613.py

@@ -0,0 +1,40 @@
+"""
+PyTTa ISO 9613-1 air absorption:
+-------------------------------
+Attenuation of sound by atmospheric absorption of ISO 9613-1:1993, for pure
+tones, from temperature, relative humidity and pressure.
+
+    >>> pytta.iso9613.attenuation(1000, 20, 50)     # dB/m
+
+The equations follow ISO 9613-1, as reproduced in Bies and Hansen,
+Engineering Noise Control (5th ed., Equations 5.69 to 5.72 and 5.75).
+"""
+import numpy as np
+
+_T0 = 293.15        # K, reference temperature
+_T01 = 273.16       # K, triple point of water
+_PR = 101.325       # kPa, reference pressure
+
+
+def attenuation(f, temperature=20., humidity=50., pressure=101.325):
+    """
+    Attenuation coefficient in dB/m at frequency f (Hz) for air at
+    `temperature` (degrees Celsius), relative `humidity` (%) and ambient
+    `pressure` (kPa). Valid, per the standard, for molar concentrations of
+    water vapour of 0.05 % to 5 % and temperatures of -20 to 50 degrees
+    Celsius.
+    """
+    f = np.asarray(f, dtype=float)
+    T = temperature + 273.15
+    pa = pressure / _PR
+    psat = 10 ** (-6.8346 * (_T01 / T) ** 1.261 + 4.6151)
+    h = humidity * psat / pa                  # molar concentration, %
+    frO = pa * (24 + 4.04e4 * h * (0.02 + h) / (0.391 + h))
+    frN = pa * (T / _T0) ** -0.5 * (
+        9 + 280 * h * np.exp(-4.170 * ((T / _T0) ** (-1 / 3) - 1)))
+    return 8.686 * f ** 2 * (
+        1.84e-11 / pa * (T / _T0) ** 0.5
+        + (T / _T0) ** -2.5 * (0.01275 * np.exp(-2239.1 / T)
+                               / (frO + f ** 2 / frO)
+                               + 0.1068 * np.exp(-3352.0 / T)
+                               / (frN + f ** 2 / frN)))

tests/test_iso354.py

@@ -0,0 +1,88 @@
+"""
+Tests for pytta.iso9613 (air absorption, ISO 9613-1) and pytta.iso354
+(sound absorption in a reverberation room, ISO 354, with the Sabine and
+Eyring reverberation formulas). Each test pins one rule.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+
+from pytta import iso354, iso9613
+
+# ISO 9613-1, pure-tone attenuation in dB/km at 20 degrees Celsius, as
+# reproduced in Long, Architectural Acoustics, Table 5.1.
+FREQS = 1000 * 10 ** (np.array([-12, -9, -6, -3, 0, 3, 6, 9]) / 10)
+TABLE = {20: [0.260, 0.712, 1.39, 2.60, 6.53, 21.5, 74.1, 215],
+         50: [0.123, 0.445, 1.32, 2.73, 4.66, 9.86, 29.4, 104],
+         100: [0.0637, 0.247, 0.895, 2.63, 5.42, 9.21, 19.4, 58.1]}
+
+
+class TestAirAbsorption(unittest.TestCase):
+
+    def test_iso_9613_1_against_tabulated_values(self):
+        for rh, expected in TABLE.items():
+            got = 1000 * iso9613.attenuation(FREQS, 20, rh)
+            dev = np.max(np.abs(got / expected - 1))
+            print(f'\n  20 C {rh} %: largest relative deviation {dev:.1e}')
+            # the table has three significant digits
+            self.assertLess(dev, 0.005, rh)
+
+
+class TestISO354(unittest.TestCase):
+
+    def test_speed_of_sound(self):
+        self.assertAlmostEqual(iso354.speed_of_sound(20), 343.0)  # Eq. (6)
+
+    def test_absorption_coefficient_round_trip(self):
+        # A room of 200 m3 with 30 m2 of absorption; adding 10 m2 of a
+        # material with alpha = 0.6 changes T as Sabine says. Equations (5)
+        # to (9) must give back alpha, also with air absorption.
+        V, A1, S, alpha = 200.0, 30.0, 10.0, 0.6
+        for f in (500.0, 4000.0):
+            m = iso354.power_attenuation(f, 20, 50)
+            c = iso354.speed_of_sound(20)
+            T1 = iso354.sabine(V, A1, c, m)
+            T2 = iso354.sabine(V, A1 + alpha * S, c, m)
+            self.assertAlmostEqual(iso354.absorption_area(T1, V, c, m), A1, 9)
+            got = iso354.absorption_coefficient(T1, T2, V, S, 20, 20, m, m)
+            self.assertAlmostEqual(got, alpha, 9)
+
+    def test_eyring_against_sabine(self):
+        # For the same mean absorption, T_Eyring / T_Sabine is
+        # alpha / -ln(1 - alpha): 1 for little absorption, 0.721 at 0.5.
+        V, S, c = 200.0, 220.0, 343.0
+        for alpha in (1e-6, 0.5):
+            ratio = iso354.eyring(V, S, alpha, c) \
+                / iso354.sabine(V, S * alpha, c)
+            self.assertAlmostEqual(ratio, alpha / -np.log(1 - alpha), 12)
+        self.assertAlmostEqual(0.5 / -np.log(0.5), 0.7213, 4)
+
+
+if __name__ == '__main__':
+    unittest.main()
+
+
+class TestISO354Limits(unittest.TestCase):
+
+    def test_table_1_scales_with_volume(self):
+        # Table 1 for 200 m3, times (V / 200)^(2/3).
+        np.testing.assert_allclose(iso354.max_room_absorption(200)[[0, 9, 10, 13, 17]],
+                                   [6.5, 6.5, 7.0, 9.5, 14.0])
+        np.testing.assert_allclose(iso354.max_room_absorption(1600), 4 * iso354.max_room_absorption(200))
+        np.testing.assert_allclose(iso354.specimen_area(1600), (10., 48.))
+        np.testing.assert_allclose(iso354.specimen_area(150), (10., 12.))
+
+    def test_room_qualification(self):
+        # A smooth curve 1 m2 below Table 1 qualifies; a dip of 20 % from
+        # the mean of its neighbours does not, nor a band above Table 1.
+        A = iso354.max_room_absorption(200) - 1
+        self.assertTrue(iso354.room_qualifies(A, 200).all())
+        dip = A.copy()
+        dip[5] *= 0.8
+        self.assertEqual(list(np.flatnonzero(~iso354.room_qualifies(dip, 200))), [5])
+        high = A.copy()
+        high[17] += 2
+        self.assertEqual(list(np.flatnonzero(~iso354.room_qualifies(high, 200))), [17])