From 8cab8a57de7442192795b87f8731175b07534307 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 00:59:49 +0100 Subject: [PATCH] [FEATURE] add air_constants, convert_room_quantity and din18041 to pytta.rooms air_constants gives the speed of sound, density and impedance of humid air as ita_constants of the ITA-Toolbox. convert_room_quantity converts between reverberation time, absorption coefficient, reflection factor and surface impedance, as ita_convert_RT_alpha_R_Z, with the same formula (Eyring or Sabine) in both directions. din18041 gives the recommended reverberation time and its tolerance band of DIN 18041:2004, as ita_din18041_reverberation_times. Tests check the constants and the DIN 18041 band against the ITA-Toolbox on golden data from tests/golden/make_room_tools_ita.m (identical), and that the conversions are consistent in both directions. Co-Authored-By: Claude Fable 5.1 --- pytta/rooms.py | 140 +++++++++++++++++++++++++++++ tests/golden/make_room_tools_ita.m | 17 ++++ tests/golden/room_tools_ita.mat | Bin 0 -> 1921 bytes tests/test_room_tools.py | 78 ++++++++++++++++ 4 files changed, 235 insertions(+) create mode 100644 tests/golden/make_room_tools_ita.m create mode 100644 tests/golden/room_tools_ita.mat create mode 100644 tests/test_room_tools.py diff --git a/pytta/rooms.py b/pytta/rooms.py index 8ab8949..d24dcd1 100644 --- a/pytta/rooms.py +++ b/pytta/rooms.py @@ -22,6 +22,9 @@ >>> pytta.rooms.strength_factor(...) >>> pytta.rooms.G_Lpe >>> pytta.rooms.G_Lps + >>> pytta.rooms.air_constants(temperature, humidity, pressure) + >>> pytta.rooms.convert_room_quantity(value, 'T', 'alpha', volume, surface) + >>> pytta.rooms.din18041(volume, purpose) Authors: João Vitor Gutkoski Paes, joao.paes@eac.ufsm.br @@ -417,3 +420,140 @@ def analyse(obj, *params, return result + +def air_constants(temperature=20., humidity=0.5, pressure=101325.): + """ + Speed of sound, density and characteristic impedance of humid air, as + ita_constants of the ITA-Toolbox: the gas constant of air with the + molar concentration of water vapour of ISO 9613-1, c = sqrt(1.4 R T) + and rho = p / (R T). + + Parameters + ---------- + temperature : float + Celsius. Default 20. + humidity : float + Relative humidity, 0 to 1. Default 0.5. + pressure : float + Static pressure in Pa. Default 101325. + + Returns + ------- + dict + 'c' (m/s), 'rho' (kg/m^3) and 'z0' (Pa s/m). The power attenuation + of air, m, is pytta.iso354.power_attenuation. + """ + if not 0 <= humidity <= 1: + raise ValueError('humidity is a relative humidity from 0 to 1.') + T = temperature + 273.15 + p_sat = 101325. * 10 ** (-6.8346 * (273.15 / T) ** 1.261 + 4.6151) + h = 100 * humidity * p_sat / pressure # molar concentration, % + R_dry = 8.31 / 0.0289644 # J/(kg K) + R = R_dry / (1 - h / 100 * (1 - R_dry / 461.)) + c = np.sqrt(1.4 * R * T) + rho = pressure / (R * T) + return {'c': float(c), 'rho': float(rho), 'z0': float(rho * c)} + + +_QUANTITIES = ('T', 'alpha', 'R', 'Z') + + +def convert_room_quantity(value, source, target, volume=None, surface=None, + formula='eyring', c=343.2, z0=413.3): + """ + Conversion between reverberation time T (s), absorption coefficient + alpha, normal-incidence reflection factor R and surface impedance Z + (Pa s/m), as ita_convert_RT_alpha_R_Z. + + T and alpha are related by Eyring, T = 24 ln(10) V / (-c S ln(1 - + alpha)), or Sabine, T = 24 ln(10) V / (c S alpha), in both directions + (ita_convert_RT_alpha_R_Z goes from T to alpha by Eyring and back by + Sabine with air absorption); alpha = 1 - |R|^2, R = sqrt(1 - alpha) + (real, the phase is not known), and Z = z0 (1 + R) / (1 - R). + + Parameters + ---------- + value : array_like + source, target : str + 'T', 'alpha', 'R' or 'Z'. + volume, surface : float + Room volume (m^3) and surface (m^2), for conversions with T. + formula : str + 'eyring' (default) or 'sabine'. + c, z0 : float + Speed of sound and characteristic impedance of air (see + air_constants). Defaults for 20 C. + """ + if source not in _QUANTITIES or target not in _QUANTITIES: + raise ValueError(f'quantities are {_QUANTITIES}.') + if formula not in ('eyring', 'sabine'): + raise ValueError("formula must be 'eyring' or 'sabine'.") + if 'T' in (source, target) and source != target and (volume is None or surface is None): + raise ValueError('volume and surface are needed for T.') + v = np.asarray(value, dtype=complex if source in ('R', 'Z') else float) + k = 24 * np.log(10) * volume / (c * surface) if volume is not None and surface is not None else None + # to alpha (or R for Z) + if source == 'T': + x = k / v + alpha = 1 - np.exp(-x) if formula == 'eyring' else x + elif source == 'alpha': + alpha = v + elif source == 'R': + alpha = 1 - np.abs(v) ** 2 + if target == 'Z': + return z0 * (1 + v) / (1 - v) + else: + R = (v / z0 - 1) / (v / z0 + 1) + if target == 'R': + return R + alpha = 1 - np.abs(R) ** 2 + if target == 'alpha': + return alpha + if target == 'T': + return k / (-np.log(1 - alpha) if formula == 'eyring' else alpha) + R = np.sqrt(1 - alpha) + return R if target == 'R' else z0 * (1 + R) / (1 - R) + + +def din18041(volume, purpose='speech', frequencies=None): + """ + Recommended reverberation time and its tolerance band by frequency, as + ita_din18041_reverberation_times (DIN 18041:2004, rooms of groups + music, speech and education): T_soll = a lg(V) + b, with (a, b) = + (0.45, 0.07), (0.37, -0.14) and (0.32, -0.17), times the tolerance + curves of the standard. + + Parameters + ---------- + volume : float + Room volume, m^3. + purpose : str + 'music', 'speech' or 'education'. + frequencies : array_like, optional + Frequencies (Hz); by default 1025 points from 0 to 22.05 kHz, as + ita_din18041_reverberation_times. + + Returns + ------- + dict + 'T_soll' (s), 'frequency' (Hz), 'T_min' and 'T_max' (s). + """ + coef = {'music': (0.45, 0.07), 'speech': (0.37, -0.14), 'education': (0.32, -0.17)} + if purpose not in coef: + raise ValueError(f'purpose must be one of {tuple(coef)}.') + f = np.linspace(0, 22050, 1025) if frequencies is None else np.asarray(frequencies, dtype=float) + a, b = coef[purpose] + T = a * np.log10(volume) + b + lg = np.log10(np.maximum(f, 1e-300)) + lo = np.full(f.size, 0.8) + if purpose == 'music': + lo[f <= 250] = (np.log10(250) - lg[f <= 250]) * 0.7 + 0.8 + lo[f >= 2000] = (np.log10(2000) - lg[f >= 2000]) * 0.5 + 0.8 + hi = lo * 1.2 / 0.8 + hi[f >= 2000] = 1.2 + else: + lo[f <= 250] = (lg[f <= 250] - np.log10(250)) * 0.5 + 0.8 + lo[f >= 2000] = (np.log10(2000) - lg[f >= 2000]) * 0.5 + 0.8 + hi = np.full(f.size, 1.2) + return {'T_soll': float(T), 'frequency': f, + 'T_min': np.maximum(lo * T, 1e-6), 'T_max': np.maximum(hi * T, 1e-6)} diff --git a/tests/golden/make_room_tools_ita.m b/tests/golden/make_room_tools_ita.m new file mode 100644 index 0000000..04309ed --- /dev/null +++ b/tests/golden/make_room_tools_ita.m @@ -0,0 +1,17 @@ +% Golden data for tests/test_room_tools.py: ita_constants and +% ita_din18041_reverberation_times of the ITA-Toolbox. Run from tests/golden +% with the ITA-Toolbox on the path. +conds = [20 0.5 101325; 0 0.0 101325; 35 0.9 95000]; +c = zeros(3,1); rho = c; z0 = c; +for k = 1:3 + [c(k), rho(k), z0(k)] = ita_constants({'c','rho_0','z_0'}, 'T', conds(k,1), 'phi', conds(k,2), 'p', conds(k,3)); + c(k) = double(c(k)); rho(k) = double(rho(k)); z0(k) = double(z0(k)); +end +r = ita_din18041_reverberation_times('v', 1000, 'purpose', 'speech'); +f = r.freqVector; din_speech = r.freqData; +r = ita_din18041_reverberation_times('v', 1000, 'purpose', 'music'); +din_music = r.freqData; +idx = 2:16:numel(f); % a subset of the 1025 bins +f = f(idx); din_speech = din_speech(idx,:); din_music = din_music(idx,:); +save('room_tools_ita.mat', 'conds', 'c', 'rho', 'z0', 'f', 'din_speech', 'din_music'); +fprintf('c %s rho %s\n', mat2str(c', 6), mat2str(rho', 6)); diff --git a/tests/golden/room_tools_ita.mat b/tests/golden/room_tools_ita.mat new file mode 100644 index 0000000000000000000000000000000000000000..5a0aff143279f7ff8ca7191704163a1ecf729703 GIT binary patch literal 1921 zcma)-c}!D@7RC!GB#2@_0TpQHRZtn0vhN@=PXt*chc`J8Or9c|8` z47D*RCmVOI;H%*TO_XbRKx{A&O*H#5v$ZiY&_vlr69QrhK`0=?40S0k0_6DDcGonu z)+tlE5F=hAB_VU)1cOSIv^x#ouykDDA3K~V!NrZOw~OURV1&yELF zZ^v4*SwvNY z9{By^a%Zz=;Gosrtfz$b7DbS}*{@+vz>UFk=5Ojr12IrrHer+(S3Ci*g9oPQTltg$ z9V(^mQ3%LvJ+7lk-5d6QF9=@~-YN%0{{;Z!Fsp_t@|)EKkaY_5dSI+h$hE3zoBdHy z=pjm{EaoJ@$?vS_^SP{`Szl;ODyEC)OZF;W*gL4qBge2+ z4^+yHIX5VK!XC}`CsJVlG;Q8U!>XFbN{adI^*soLb1qNaFQ&3Z)!|w$yKFKjTbW#) zQ9xq);}g0YEcVDCF21Pnee!Ol<)^k6nYa30f>RW(p;uvoK<(ImYI;@c721-=0c8>bjt!<+{LURB4iRXCEx-PEF^m^#Aj*42lQg&nuXD=1bFls zyD|?6jGs`BzcTEbIOI=t;?rkbwbRHN;aQg+(k5;XqBA5B%K6LZzq-lnpWGz=%}ui@ z!^3kjKbbin78`EM0L0-OQ$@aFa|9&Z>N${7-!xp4Cg-QaXXh$qXaS@M=0~{Sjxwze zu{^L~IWFlM9L1{3D(Bp&pSztWGIZe^^CI2hByh~?XD#EW)?f`Is_*es#Fi?nxFwEp z4cJ9PmNC}GexlX){*l}cWQ@8+J+}IywCoT$v~9_KVtFA$mT#Uxtj*Ctj)f!~Vd!0^ z;F)V~&1yj(C{SLoLyB9) zodQswm|>C^#IA2N7=(-;@*XpKvo-~tVG(1+wFBA*rC$Uv4Yyp!Il!nOVT?)4w}P7) zpEud@PX40;>ot-5b0M|g|F;5DbP~%D%B|n4^!_-?&+Dt-nEpMfpIVjvgMi7v3X0Qp zMq2!mKTtcV@9L(s)Qc$<$EB@S3|cItTlGS?YD0puUuZC?GA(UB`pfVi~dWp<|zBHWLiQ~2s2 z8?mqdNUO)?bZV5cEy{$y&)>o3JsdTOsf1RI8Lf_X6o7#t6P_SD`F-D8kaS$oEUsjOZOa`v+{n Bym alpha -> T returns T by both formulas; alpha = 1 - R^2 and + # Z = z0 (1 + R) / (1 - R), so R = 0 is alpha = 1 and Z = z0; Eyring + # tends to Sabine for small alpha. + V, S, T = 200., 210., np.array([0.5, 1.0, 2.0]) + for formula in ('eyring', 'sabine'): + a = rooms.convert_room_quantity(T, 'T', 'alpha', V, S, formula) + back = rooms.convert_room_quantity(a, 'alpha', 'T', V, S, formula) + np.testing.assert_allclose(back, T, rtol=1e-12) + R = rooms.convert_room_quantity(0.75, 'alpha', 'R') + Z = rooms.convert_room_quantity(R, 'R', 'Z') + print(f'\n alpha 0.75 -> R {R:.4f} -> Z {abs(Z):.2f} -> alpha ' + f'{rooms.convert_room_quantity(Z, "Z", "alpha"):.6f}') + self.assertAlmostEqual(float(R), 0.5) + self.assertAlmostEqual(float(abs(Z)), 3 * 413.3, places=9) + self.assertAlmostEqual(float(rooms.convert_room_quantity(413.3, 'Z', 'alpha')), 1) + small_e = rooms.convert_room_quantity(1e-3, 'alpha', 'T', V, S, 'eyring') + small_s = rooms.convert_room_quantity(1e-3, 'alpha', 'T', V, S, 'sabine') + self.assertAlmostEqual(small_e / small_s, 1, delta=1e-3) + with self.assertRaises(ValueError): + rooms.convert_room_quantity(1., 'T', 'alpha') + + def test_din18041_matches_ita(self): + # T_soll = 0.37 lg V - 0.14 for speech; the tolerance band of + # ita_din18041_reverberation_times for speech and music. + d = golden() + out = rooms.din18041(1000, 'speech', d['f']) + music = rooms.din18041(1000, 'music', d['f']) + e = max(np.abs(out['T_min'] - d['din_speech'][:, 0]).max(), + np.abs(out['T_max'] - d['din_speech'][:, 1]).max(), + np.abs(music['T_min'] - d['din_music'][:, 0]).max(), + np.abs(music['T_max'] - d['din_music'][:, 1]).max()) + print(f"\n T_soll {out['T_soll']:.4f} s for 1000 m^3; band against ITA {e:.1e}") + self.assertAlmostEqual(out['T_soll'], 0.37 * 3 - 0.14) + self.assertLess(e, 1e-12) + + +if __name__ == '__main__': + unittest.main()