From f536de18c26982d957471d295ff51fd13be77c7f Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 01:19:22 +0100 Subject: [PATCH] [FEATURE] add image_source_ir to pytta.rooms image_source_ir gives the impulse response of a rectangular room by the image source method of Allen and Berkley (1979), as ita_roomacoustics_imagesources of the ITA-Toolbox: the free-field Green's function of every image source within the duration, weighted by the reflection factors of the six walls it goes through, each arrival placed between samples by a Hann-windowed sinc. Tests check against sums of spherical waves: the anechoic case, one reflecting wall, two walls with different factors and their corner image, and the opposite wall. Co-Authored-By: Claude Fable 5.1 --- pytta/rooms.py | 74 +++++++++++++++++++++++++++++++++++ tests/test_image_sources.py | 78 +++++++++++++++++++++++++++++++++++++ 2 files changed, 152 insertions(+) create mode 100644 tests/test_image_sources.py diff --git a/pytta/rooms.py b/pytta/rooms.py index 8ab8949..d147d83 100644 --- a/pytta/rooms.py +++ b/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') diff --git a/tests/test_image_sources.py b/tests/test_image_sources.py new file mode 100644 index 0000000..7804d27 --- /dev/null +++ b/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()