[FEATURE] add image_source_ir to pytta.rooms

pytta/rooms.py

@@ -22,6 +22,7 @@
     >>> pytta.rooms.strength_factor(...)
     >>> pytta.rooms.G_Lpe
     >>> pytta.rooms.G_Lps
+    >>> pytta.rooms.image_source_ir(dimensions, source, receiver, fs, duration)
 
 Authors:
     João Vitor Gutkoski Paes, joao.paes@eac.ufsm.br
@@ -417,3 +418,76 @@ def analyse(obj, *params,
 
     return result
 
+
+def image_source_ir(dimensions, source, receiver, samplingRate, duration,
+                    reflection=1., c=343.2, half_width=8):
+    """
+    Impulse response of a rectangular room by the image source method of
+    Allen and Berkley (1979), as ita_roomacoustics_imagesources: the free
+    field Green's function delta(t - d/c) / (4 pi d) of every image source
+    within duration, weighted by the reflection factors of the walls it
+    goes through.
+
+    Parameters
+    ----------
+    dimensions, source, receiver : array_like
+        (Lx, Ly, Lz) and the positions (m), with 0 <= x <= Lx.
+    samplingRate : int
+    duration : float
+        Length of the response (s); image sources farther than c duration
+        are left out.
+    reflection : float or array_like
+        Pressure reflection factor of every wall, or six of them for the
+        walls x = 0, x = Lx, y = 0, y = Ly, z = 0, z = Lz (frequency and
+        angle independent; ita_roomacoustics_imagesources derives them
+        from a wall impedance Z as (Z - rho c) / (Z + rho c)). Default 1.
+    half_width : int
+        Half length in samples of the Hann-windowed sinc that places each
+        arrival between samples. Default 8.
+
+    Returns
+    -------
+    SignalObj
+        The response; convolved with rho dQ/dt of the source (Q the volume
+        velocity) it gives the pressure.
+    """
+    from pytta.classes import SignalObj
+    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.')
+    beta = np.broadcast_to(np.asarray(reflection, dtype=float), (6,))
+    n = int(round(duration * samplingRate))
+    out = np.zeros(n + 2 * half_width + 1)
+    reach = c * duration
+    lmax = [int(np.ceil(reach / (2 * L[i]))) + 1 for i in range(3)]
+    k = np.arange(-half_width, half_width + 1)
+    # per axis: image coordinates, distances and reflection weights of
+    # x = 2 l Lx + (1 - 2u) x0 (u = 0, 1): |l - u| reflections at x = 0 and
+    # |l| at x = Lx
+    axes = []
+    for i in range(3):
+        ls = np.arange(-lmax[i], lmax[i] + 1)
+        coord, weight = [], []
+        for u in (0, 1):
+            coord.append(2 * ls * L[i] + (1 - 2 * u) * r0[i] - r[i])
+            weight.append(beta[2 * i] ** np.abs(ls - u) * beta[2 * i + 1] ** np.abs(ls))
+        axes.append((np.concatenate(coord), np.concatenate(weight)))
+    (dx, wx), (dy, wy), (dz, wz) = axes
+    d2xy = dx[:, None] ** 2 + dy[None, :] ** 2
+    wxy = wx[:, None] * wy[None, :]
+    for z, w in zip(dz, wz):
+        d = np.sqrt(d2xy + z ** 2).ravel()
+        keep = d <= reach
+        if not np.any(keep):
+            continue
+        d, a = d[keep], (wxy.ravel() * w)[keep] / (4 * np.pi * d[keep])
+        t = d / c * samplingRate
+        i0 = np.floor(t).astype(int)
+        frac = t - i0
+        x = k[None, :] - frac[:, None]
+        kernel = np.sinc(x) * 0.5 * (1 + np.cos(np.pi * x / (half_width + 1)))
+        idx = i0[:, None] + k[None, :] + half_width
+        np.add.at(out, idx.ravel(), (a[:, None] * kernel).ravel())
+    return SignalObj(out[half_width:half_width + n], 'time', samplingRate,
+                     signalType='energy')

tests/test_image_sources.py

@@ -0,0 +1,78 @@
+"""
+Tests for pytta.rooms.image_source_ir, the image source method of Allen and
+Berkley (1979) for a rectangular room, after ita_roomacoustics_imagesources
+of the ITA-Toolbox. Each test pins one rule.
+"""
+import unittest
+
+import numpy as np
+
+from pytta import rooms
+
+C, FS = 343.2, 16000
+
+
+def spectrum(ir):
+    h = ir.timeSignal[:, 0].astype(float)
+    return np.fft.rfftfreq(h.size, 1 / FS), np.fft.rfft(h)
+
+
+class TestImageSources(unittest.TestCase):
+
+    def test_anechoic_and_one_wall(self):
+        # With no reflection the response is the free-field Green's
+        # function exp(-j w d / c) / (4 pi d); with only the wall x = 0
+        # reflecting (factor 1) the mirror image at (-x0, y0, z0) adds its
+        # own spherical wave. Compared below FS / 4, where the windowed sinc
+        # that places each arrival is flat.
+        label = 'image'
+        Lr, s, r = (6., 5., 4.), (1.2, 2.0, 1.5), (4.1, 3.2, 2.6)
+        d0 = np.linalg.norm(np.subtract(s, r))
+        d1 = np.linalg.norm(np.subtract((-s[0], s[1], s[2]), r))
+        free = rooms.image_source_ir(Lr, s, r, FS, 0.05, reflection=0)
+        wall = rooms.image_source_ir(Lr, s, r, FS, 0.05, reflection=[1, 0, 0, 0, 0, 0])
+        f, Hf = spectrum(free)
+        _, Hw = spectrum(wall)
+        band = (f > 50) & (f < FS / 4)
+        w = 2 * np.pi * f[band]
+        g0 = np.exp(-1j * w * d0 / C) / (4 * np.pi * d0)
+        g1 = np.exp(-1j * w * d1 / C) / (4 * np.pi * d1)
+        e_free = np.abs(Hf[band] - g0).max() / np.abs(g0).max()
+        e_wall = np.abs(Hw[band] - g0 - g1).max() / np.abs(g0).max()
+        print(f'\n  anechoic {e_free:.1e}, one wall {e_wall:.1e} relative to 1 / (4 pi d)')
+        self.assertLess(e_free, 2e-3)
+        self.assertLess(e_wall, 2e-3)
+
+    def test_each_wall_weights_its_images(self):
+        # Walls x = 0 and y = 0 reflecting with 0.5 and 0.8 and the others
+        # absorbing: the response is the direct sound, the images across
+        # each wall (weights 0.5 and 0.8) and the corner image (0.5 x 0.8).
+        # The wall x = Lx alone puts the image at 2 Lx - x0.
+        label = 'walls'
+        Lr, s, r = (6., 5., 4.), (1.2, 2.0, 1.5), (4.1, 3.2, 2.6)
+        f, H = spectrum(rooms.image_source_ir(Lr, s, r, FS, 0.05, [0.5, 0, 0.8, 0, 0, 0]))
+        _, Hx = spectrum(rooms.image_source_ir(Lr, s, r, FS, 0.05, [0, 1, 0, 0, 0, 0]))
+        band = (f > 50) & (f < FS / 4)
+        w = 2 * np.pi * f[band]
+
+        def wave(pos, weight):
+            d = np.linalg.norm(np.subtract(pos, r))
+            return weight * np.exp(-1j * w * d / C) / (4 * np.pi * d)
+        x0, y0, z0 = s
+        expected = (wave(s, 1) + wave((-x0, y0, z0), 0.5) + wave((x0, -y0, z0), 0.8)
+                    + wave((-x0, -y0, z0), 0.4))
+        expected_x = wave(s, 1) + wave((2 * Lr[0] - x0, y0, z0), 1)
+        ref = np.abs(wave(s, 1)).max()
+        e = np.abs(H[band] - expected).max() / ref
+        ex = np.abs(Hx[band] - expected_x).max() / ref
+        print(f'\n  two walls with their corner {e:.1e}, wall x = Lx {ex:.1e}')
+        self.assertLess(e, 2e-3)
+        self.assertLess(ex, 2e-3)
+
+    def test_rejects_positions_outside(self):
+        with self.assertRaises(ValueError):
+            rooms.image_source_ir((5, 4, 3), (5.5, 1, 1), (1, 1, 1), FS, 0.1)
+
+
+if __name__ == '__main__':
+    unittest.main()