From 586550537b9723e1a9aba6cc3e3fee2fad80c8b3 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 01:16:29 +0100 Subject: [PATCH] [FEATURE] add rectangular_room_frf and rectangular_room_ir to pytta.rooms rectangular_room_frf gives the transfer function from a point source of volume velocity to a receiver in a rectangular room with rigid walls and uniform modal damping, by modal superposition (Fahy, Foundations of Engineering Acoustics, 9.8), as ita_roomacoustics_analytic_FRF of the ITA-Toolbox. The uniform mode is computed in closed form, so the response is finite at 0 Hz. rectangular_room_ir gives its impulse response, band-limited with a raised-cosine fade below the Nyquist frequency. Tests check the frequency and half-power bandwidth of an isolated mode (2.2 / T), reciprocity, the nodal plane of a mode, and that the impulse response decays with the reverberation time of the model with no floor. Co-Authored-By: Claude Fable 5.1 --- pytta/rooms.py | 90 ++++++++++++++++++++++++++++++++++++++++ tests/test_room_modes.py | 77 ++++++++++++++++++++++++++++++++++ 2 files changed, 167 insertions(+) create mode 100644 tests/test_room_modes.py diff --git a/pytta/rooms.py b/pytta/rooms.py index 8ab8949..9e64dcc 100644 --- a/pytta/rooms.py +++ b/pytta/rooms.py @@ -22,6 +22,8 @@ >>> pytta.rooms.strength_factor(...) >>> pytta.rooms.G_Lpe >>> pytta.rooms.G_Lps + >>> pytta.rooms.rectangular_room_frf(dimensions, source, receiver, freqs, T) + >>> pytta.rooms.rectangular_room_ir(dimensions, source, receiver, fs, nSamples, T) Authors: João Vitor Gutkoski Paes, joao.paes@eac.ufsm.br @@ -417,3 +419,91 @@ def analyse(obj, *params, return result + +def _room_modes(dimensions, f_max, c): + """Mode numbers (n, 3) and frequencies of a rigid rectangular room up to f_max.""" + L = np.asarray(dimensions, dtype=float) + kmax = 2 * f_max / c + nx = np.arange(int(np.floor(kmax * L[0])) + 1) + ny = np.arange(int(np.floor(kmax * L[1])) + 1) + nz = np.arange(int(np.floor(kmax * L[2])) + 1) + n = np.stack(np.meshgrid(nx, ny, nz, indexing='ij'), -1).reshape(-1, 3) + f = c / 2 * np.sqrt(np.sum((n / L) ** 2, axis=1)) + keep = f <= f_max + return n[keep], f[keep] + + +def rectangular_room_frf(dimensions, source, receiver, frequencies, T=2., + f_max=None, c=343.2, rho=1.204): + """ + Transfer function p / Q (Pa per m^3/s) from a point source of volume + velocity Q to a receiver in a rectangular room with rigid walls and + uniform modal damping, by modal superposition (Fahy, Foundations of + Engineering Acoustics, 9.8; Kuttruff, Room Acoustics, chapter 3), as + ita_roomacoustics_analytic_FRF: + + H = j w rho c^2 sum_N psi_N(r) psi_N(r0) / (L_N (w_N^2 - w^2 + 2 j d w)) + + with psi_N = cos(nx pi x / Lx) cos(ny pi y / Ly) cos(nz pi z / Lz), + L_N = V prod(1 for n = 0, 1/2 otherwise), w_N = pi c sqrt(sum (n/L)^2) + and d = 3 ln(10) / T for every mode. Time dependence exp(j w t), the + convention of numpy.fft, so the impulse response is causal. + + Parameters + ---------- + dimensions, source, receiver : array_like + (Lx, Ly, Lz) of the room and the positions (m), with 0 <= x <= Lx. + frequencies : array_like + Frequencies (Hz). + T : float + Reverberation time (s) that sets the damping. Default 2. + f_max : float, optional + Highest mode frequency; by default 1.2 times the highest frequency. + """ + L = np.asarray(dimensions, dtype=float) + r0, r = np.asarray(source, dtype=float), np.asarray(receiver, dtype=float) + if np.any(r0 < 0) or np.any(r0 > L) or np.any(r < 0) or np.any(r > L): + raise ValueError('source and receiver must lie inside the room.') + f = np.asarray(frequencies, dtype=float) + f_max = 1.2 * np.max(f) if f_max is None else f_max + n, fn = _room_modes(L, f_max, c) + psi = np.prod(np.cos(np.pi * n * r0 / L), axis=1) * np.prod(np.cos(np.pi * n * r / L), axis=1) + lam = np.prod(L) * 0.5 ** np.sum(n > 0, axis=1) + delta = 3 * np.log(10) / T + # the (0, 0, 0) mode in closed form, rho c^2 / (V (2 d + j w)), finite at + # f = 0 (the compression of the whole volume); the other modes vanish there + w = 2 * np.pi * f + H = rho * c ** 2 / (np.prod(L) * (2 * delta + 1j * w)) + rest = np.any(n > 0, axis=1) + wn2, a = (2 * np.pi * fn[rest]) ** 2, (psi / lam)[rest] + step = max(1, 2_000_000 // max(1, a.size)) # bounds the memory per block + for k in range(0, f.size, step): + wk = w[k:k + step, None] + H[k:k + step] += 1j * wk[:, 0] * rho * c ** 2 * np.sum( + a / (wn2 - wk ** 2 + 2j * delta * wk), axis=1) + return H + + +def rectangular_room_ir(dimensions, source, receiver, samplingRate, numSamples, + T=2., f_max=None, c=343.2, rho=1.204): + """ + Impulse response of rectangular_room_frf at the FFT bins of numSamples + samples, as a SignalObj, so that its convolution with the volume + velocity gives the pressure. + + The modal sum is band-limited: modes up to f_max (0.4 samplingRate by + default) and a raised-cosine fade of the response from f_max to + samplingRate / 2. Cutting the response where it is large, at the + Nyquist frequency, would ring through the whole circular response. The + band-limited direct sound still rings a little before it arrives; the + part before t = 0 wraps to the last samples of the circular response. + """ + from pytta.classes import SignalObj + f = np.fft.rfftfreq(numSamples, 1 / samplingRate) + f_max = 0.4 * samplingRate if f_max is None else f_max + H = rectangular_room_frf(dimensions, source, receiver, f, T, f_max, c, rho) + nyq = samplingRate / 2 + fade = np.clip((f - f_max) / max(nyq - f_max, 1e-12), 0, 1) + H = H * 0.5 * (1 + np.cos(np.pi * fade)) + return SignalObj(np.fft.irfft(H, numSamples), 'time', samplingRate, + signalType='energy') diff --git a/tests/test_room_modes.py b/tests/test_room_modes.py new file mode 100644 index 0000000..be91ba7 --- /dev/null +++ b/tests/test_room_modes.py @@ -0,0 +1,77 @@ +""" +Tests for the modal model of a rectangular room in pytta.rooms +(rectangular_room_frf, rectangular_room_ir), after +ita_roomacoustics_analytic_FRF of the ITA-Toolbox. Each test pins one rule +of the modal superposition (Fahy, Foundations of Engineering Acoustics, +9.8). +""" +import unittest + +import numpy as np + +from pytta import rooms + +L = (5., 4., 3.) +C = 343.2 + + +class TestRoomModes(unittest.TestCase): + + def test_isolated_mode_frequency_and_bandwidth(self): + # The (1, 0, 0) mode sits at c / (2 Lx) = 34.32 Hz, alone up to the + # (0, 1, 0) mode at 42.9 Hz; its half-power bandwidth is + # d / pi = 3 ln(10) / (pi T) = 2.2 / T. + T = 2. + f = np.arange(25, 39, 0.001) + H = np.abs(rooms.rectangular_room_frf(L, (0, 2, 1.5), (5, 2, 1.5), f, T, f_max=200)) + peak = f[np.argmax(H)] + band = f[H >= H.max() / np.sqrt(2)] + width = band[-1] - band[0] + print(f'\n peak {peak:.3f} Hz (mode {C / 10:.3f} Hz), half-power width {width:.4f} Hz ' + f'(3 ln 10 / (pi T) = {3 * np.log(10) / (np.pi * T):.4f} Hz)') + self.assertAlmostEqual(peak, C / 10, delta=0.05) + self.assertAlmostEqual(width, 3 * np.log(10) / (np.pi * T), delta=0.05) + + def test_reciprocity_and_nodal_plane(self): + # Swapping source and receiver leaves H unchanged; a receiver on + # the plane x = Lx / 2 does not excite the (1, 0, 0) mode. + f = np.linspace(10, 120, 500) + a = rooms.rectangular_room_frf(L, (0.7, 1.1, 0.4), (3.9, 2.7, 2.2), f) + b = rooms.rectangular_room_frf(L, (3.9, 2.7, 2.2), (0.7, 1.1, 0.4), f) + fm = np.array([C / 10]) + at_node = abs(rooms.rectangular_room_frf(L, (0, 2, 1.5), (2.5, 2, 1.5), fm, f_max=60)[0]) + off_node = abs(rooms.rectangular_room_frf(L, (0, 2, 1.5), (5, 2, 1.5), fm, f_max=60)[0]) + print(f'\n reciprocity {np.abs(a - b).max():.1e}; at the mode, receiver on the nodal plane ' + f'{20 * np.log10(at_node / off_node):.1f} dB below the wall') + self.assertLess(np.abs(a - b).max() / np.abs(a).max(), 1e-12) + self.assertLess(at_node / off_node, 0.1) + + def test_impulse_response_decays_with_T(self): + # The modal damping d = 3 ln(10) / T makes the energy fall 60 dB in + # T. The decay curve of the energy summed over six receivers (one + # receiver beats between modes) gives T back from -5 to -35 dB, and + # it keeps falling with no floor: from 2.0 s to 2.5 s it lies more + # than 100 dB below its first 0.25 s. + label = 'decay' + T, fs, n = 0.8, 2000, 6000 + rng = np.random.default_rng(2) + e = 0 + for k in range(6): + r = rng.uniform((0.3, 0.3, 0.3), (4.7, 3.7, 2.7)) + ir = rooms.rectangular_room_ir(L, (0.7, 1.1, 0.4), r, fs, n, T) + e = e + ir.timeSignal[:, 0].astype(float) ** 2 + edc = 10 * np.log10(np.cumsum(e[::-1])[::-1] / e.sum()) + t = np.arange(n) / fs + T30 = 2 * (t[np.argmax(edc < -35)] - t[np.argmax(edc < -5)]) + tail = 10 * np.log10(e[4000:5000].mean() / e[:500].mean()) + print(f'\n T from the decay {T30:.4f} s, T of the model {T} s; 2.0 to 2.5 s at {tail:.1f} dB') + self.assertAlmostEqual(T30 / T, 1, delta=0.05) + self.assertLess(tail, -100) + + def test_rejects_positions_outside(self): + with self.assertRaises(ValueError): + rooms.rectangular_room_frf(L, (6, 1, 1), (1, 1, 1), [100.]) + + +if __name__ == '__main__': + unittest.main()