[FEATURE] add rectangular_room_frf and rectangular_room_ir to pytta.rooms

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')

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()