[FEATURE] add pytta.iso9613 and pytta.iso354 for air and room absorption

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,61 @@
+"""
+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, and the Sabine and Eyring reverberation formulas.
+
+    >>> 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))

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,65 @@
+"""
+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()