diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..fabcf26 100644 --- a/pytta/__init__.py +++ b/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 diff --git a/pytta/iso354.py b/pytta/iso354.py new file mode 100644 index 0000000..58e5965 --- /dev/null +++ b/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)) diff --git a/pytta/iso9613.py b/pytta/iso9613.py new file mode 100644 index 0000000..d8e5f49 --- /dev/null +++ b/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))) diff --git a/tests/test_iso354.py b/tests/test_iso354.py new file mode 100644 index 0000000..4e66e6b --- /dev/null +++ b/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()