diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..57877a0 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -57,6 +57,9 @@ Modules: * pytta.iso3741: calculations according to the standard; + * pytta.vbap: + vector base amplitude panning for 2D and 3D loudspeaker layouts; + .. toctree:: classes/index utils/index @@ -64,4 +67,5 @@ Modules: functions generate rooms + vbap properties diff --git a/docs/pyttadoc/vbap.rst b/docs/pyttadoc/vbap.rst new file mode 100644 index 0000000..3d1b0d0 --- /dev/null +++ b/docs/pyttadoc/vbap.rst @@ -0,0 +1,7 @@ +.. currentmodule:: pytta + +VBAP +==== + +.. automodule:: pytta.vbap + :members: diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..908e6ab 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -72,6 +72,9 @@ * pytta.iso3741: calculations according to the standard; + * pytta.vbap: + vector base amplitude panning for 2D and 3D loudspeaker layouts; + You can find usage examples in the toolbox's examples folder. @@ -108,6 +111,7 @@ from . import utils from . import rooms from . import iso3741 +from . import vbap from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +134,7 @@ 'generate', 'utils', 'iso3741', + 'vbap', 'rooms', # Functions diff --git a/pytta/vbap.py b/pytta/vbap.py new file mode 100644 index 0000000..6448614 --- /dev/null +++ b/pytta/vbap.py @@ -0,0 +1,254 @@ +""" +PyTTa VBAP: +----------- + +Vector base amplitude panning (VBAP) of Pulkki (1997) for 2D and 3D +loudspeaker layouts, and its spread version, multiple-direction amplitude +panning (MDAP). + +Available functions: + + >>> groups = pytta.vbap.loudspeaker_pairs(azimuths) # 2D + >>> groups = pytta.vbap.loudspeaker_triplets(directions) # 3D + >>> inverses = pytta.vbap.inverse_matrices(directions, groups) + >>> G = pytta.vbap.gains(source_directions, groups, inverses, spread=0) + +Directions are in degrees: azimuth for 2D layouts, [azimuth, elevation] for +3D ones, with x = cos(azi) cos(elev), y = sin(azi) cos(elev) and z = +sin(elev) (0 deg in front, 90 deg to the left). Loudspeaker indices are +0-based. + +Origin of the code +------------------ + +This module is a port to Python of the MATLAB library +Vector-Base-Amplitude-Panning by Archontis Politis and Ville Pulkki (2015), +https://github.com/polarch/Vector-Base-Amplitude-Panning (findLsPairs.m, +findLsTriplets.m, invertLsMtx.m, vbap.m, getSpreadSrcDirs.m), under the BSD +3-Clause licence reproduced below. The gains follow the MATLAB code +operation by operation. + +References +---------- + +Pulkki, V. (1997). Virtual sound source positioning using vector base +amplitude panning. Journal of the Audio Engineering Society, 45(6), 456-466. + +Pulkki, V. (1999). Uniform spreading of amplitude panned virtual sources. +Proceedings of the IEEE Workshop on Applications of Signal Processing to +Audio and Acoustics, New Paltz, NY. + +Licence of the ported code +-------------------------- + + Copyright (c) 2015, Archontis Politis + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are + met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + * Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED + TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +""" + +import numpy as np +from scipy.spatial import ConvexHull + + +def _cartesian(directions): + d = np.radians(np.asarray(directions, dtype=float)) + if d.ndim == 1: + return np.column_stack([np.cos(d), np.sin(d)]) + azi, elev = d[:, 0], d[:, 1] + return np.column_stack([np.cos(azi) * np.cos(elev), + np.sin(azi) * np.cos(elev), np.sin(elev)]) + + +def loudspeaker_pairs(azimuths, max_aperture=None): + """ + Adjacent loudspeaker pairs of a 2D layout (findLsPairs). + + Parameters + ---------- + azimuths : array_like + Loudspeaker azimuths in degrees. + max_aperture : float, optional + Pairs that span this angle in degrees or more are left out. + + Returns + ------- + numpy.ndarray + (pairs, 2) loudspeaker indices, around the circle by azimuth. + """ + azimuths = np.asarray(azimuths, dtype=float) + order = np.argsort(azimuths, kind='stable') + groups = np.column_stack([order, np.roll(order, -1)]) + if max_aperture is not None: + s = azimuths[np.append(order, order[0])] + apertures = np.diff(s) + apertures[-1] = 360 - abs(apertures[-1]) + groups = groups[apertures < max_aperture] + return groups + + +def loudspeaker_triplets(directions, max_aperture=None): + """ + Loudspeaker triplets of a 3D layout (findLsTriplets). + + The triplets are the faces of the convex hull of the loudspeaker + directions, which on the unit sphere is their Delaunay triangulation. + Faces whose outward normal is at 90 deg or more from their centroid are + left out: those are the faces that close a layout with no loudspeakers + below or above, through which no source direction passes. + + Parameters + ---------- + directions : array_like + (loudspeakers, 2) [azimuth, elevation] in degrees. + max_aperture : float, optional + Triplets with a side longer than this angle in degrees are left out. + + Returns + ------- + numpy.ndarray + (triplets, 3) loudspeaker indices, each triplet counter-clockwise + seen from outside and starting at its smallest index, sorted. + """ + U = _cartesian(directions) + faces = [] + for face in ConvexHull(U).simplices: + v = U[face] + n = np.cross(v[1] - v[0], v[2] - v[1]) + if np.dot(n, v.mean(axis=0) - U.mean(axis=0)) < 0: + face, v, n = face[::-1], v[::-1], -n + # MATLAB keeps acos(dot(centroid, n)) < pi/2; a dot product of + # rounding size, as on a flat bottom face, rounds to pi/2 there. + if np.dot(v.mean(axis=0), n) > 1e-12: + faces.append(np.roll(face, -int(np.argmin(face)))) + faces = np.array(faces, dtype=int) + if max_aperture is not None: + sides = np.stack([np.sum(U[faces[:, a]] * U[faces[:, b]], axis=1) + for a, b in ((0, 1), (1, 2), (2, 0))], axis=1) + faces = faces[np.all(np.arccos(np.clip(sides, -1, 1)) + <= np.radians(max_aperture), axis=1)] + return faces[np.lexsort((faces[:, 1], faces[:, 0]))] + + +def inverse_matrices(directions, groups): + """ + Inverse of the matrix of loudspeaker vectors of each group (invertLsMtx). + + Parameters + ---------- + directions : array_like + Loudspeaker azimuths (2D) or [azimuth, elevation] (3D) in degrees. + groups : numpy.ndarray + Pairs or triplets from loudspeaker_pairs or loudspeaker_triplets. + + Returns + ------- + numpy.ndarray + (groups, dim, dim): the inverse of the matrix whose rows are the + unit vectors of the group, so that the gains of a source u are + u @ inverse. + """ + U = _cartesian(directions) + return np.linalg.inv(U[np.asarray(groups)]) + + +def _spread_directions(direction, spread, num_src=8): + # getSpreadSrcDirs with its defaults: num_src sources, one ring in 3D. + if np.ndim(direction) == 0: + d = np.radians(np.linspace(direction - spread / 2, direction + spread / 2, + num_src + 1)) + return np.column_stack([np.cos(d), np.sin(d)]) + u = _cartesian(np.atleast_2d(direction))[0] + u_x = np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]]) + theta = 2 * np.pi / num_src + R = (np.eye(3) * np.cos(theta) + np.sin(theta) * u_x + + (1 - np.cos(theta)) * np.outer(u, u)) + base = np.zeros((num_src, 3)) + if direction[1] > 89 or direction[1] < -89: + base[0] = [1, 0, 0] + else: + c = np.cross(u, [0, 0, 1]) + base[0] = c / np.sqrt(np.sum(c**2)) + for n in range(1, num_src): + base[n] = R @ base[n - 1] + ring = u + base * np.tan(np.radians(spread / 2)) + ring = ring / np.sqrt(np.sum(ring**2, axis=1))[:, None] + return np.vstack([ring, u]) + + +def _group_gains(u, groups, inverses, first_only): + g = np.zeros(int(groups.max()) + 1) + for group, inv in zip(groups, inverses): + g_tmp = u @ inv + if np.min(g_tmp) > -0.001: + g[group] = g_tmp / np.sqrt(np.sum(g_tmp**2)) + if first_only: + break + return g + + +def gains(directions, groups, inverses, spread=0.): + """ + VBAP gains of a set of source directions (vbap). + + Parameters + ---------- + directions : array_like + Source azimuths (2D) or (sources, 2) [azimuth, elevation] (3D) in + degrees. + groups : numpy.ndarray + Pairs or triplets from loudspeaker_pairs or loudspeaker_triplets. + inverses : numpy.ndarray + From inverse_matrices. + spread : float, optional + Spread in degrees for MDAP: the gains of 9 directions across the + spread (8 on a ring and the source in 3D) are summed. The default is + 0, plain VBAP. + + Returns + ------- + numpy.ndarray + (sources, loudspeakers) gains, non-negative, of norm 1 per source. + A source outside every group (below a hemisphere layout, say) gets + all gains 0, as in the MATLAB code. + """ + groups = np.asarray(groups) + three_d = groups.shape[1] == 3 + directions = np.asarray(directions, dtype=float) + directions = directions.reshape(-1, 2) if three_d else directions.reshape(-1) + out = np.zeros((len(directions), int(groups.max()) + 1)) + for n, d in enumerate(directions): + if spread: + g = sum(_group_gains(u, groups, inverses, three_d) + for u in _spread_directions(d, spread)) + else: + g = _group_gains(_cartesian(np.atleast_2d(d) if three_d else [d])[0], + groups, inverses, three_d) + with np.errstate(invalid='ignore'): + out[n] = np.fmax(0, g / np.sqrt(np.sum(g**2))) + return out diff --git a/tests/golden/make_vbap_polarch.m b/tests/golden/make_vbap_polarch.m new file mode 100644 index 0000000..e6abb89 --- /dev/null +++ b/tests/golden/make_vbap_polarch.m @@ -0,0 +1,20 @@ +% Golden data for tests/test_vbap.py: polarch Vector-Base-Amplitude-Panning +% (Archontis Politis, BSD 3-Clause), the copy that Iceberg installs in +% Toolboxes/polarch/vbap. Run from tests/golden with that folder on the path. +% Indices are 1-based. +ls2 = [0 30 -30 110 -110]; % 5.0, 2D +ls3 = [0 0; 90 0; 180 0; -90 0; 45 45; 135 45; -135 45; -45 45; 0 90]; % hemisphere +src2 = (-180:7.5:180)'; +[a, e] = meshgrid(-180:20:160, -10:20:90); src3 = [a(:) e(:)]; +pairs = findLsPairs(ls2); +inv2 = invertLsMtx(ls2, pairs); +g2 = vbap(src2, pairs, inv2); +g2s = vbap(src2, pairs, inv2, 20); +tri = findLsTriplets(ls3); +inv3 = invertLsMtx(ls3, tri); +g3 = vbap(src3, tri, inv3); +g3s = vbap(src3, tri, inv3, 20); +save('vbap_polarch.mat', 'ls2', 'ls3', 'src2', 'src3', 'pairs', 'inv2', 'g2', ... + 'g2s', 'tri', 'inv3', 'g3', 'g3s'); +fprintf('%d pairs, %d triplets\n', size(pairs, 1), size(tri, 1)); +disp(tri) diff --git a/tests/golden/vbap_polarch.mat b/tests/golden/vbap_polarch.mat new file mode 100644 index 0000000..8a568a5 Binary files /dev/null and b/tests/golden/vbap_polarch.mat differ diff --git a/tests/test_vbap.py b/tests/test_vbap.py new file mode 100644 index 0000000..00e5d13 --- /dev/null +++ b/tests/test_vbap.py @@ -0,0 +1,90 @@ +""" +Tests for pytta.vbap, vector base amplitude panning (Pulkki 1997), a port of +the VBAP library of Archontis Politis. Each test pins one rule. + +The golden data come from golden/make_vbap_polarch.m (MATLAB R2026b with +polarch Vector-Base-Amplitude-Panning, the copy that Iceberg installs); +their loudspeaker indices are 1-based. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +from pytta import vbap + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden') + + +def unit(azi, elev=0.): + a, e = np.radians(azi), np.radians(elev) + return np.array([np.cos(a) * np.cos(e), np.sin(a) * np.cos(e), np.sin(e)]) + + +class TestVBAP(unittest.TestCase): + + def test_stereo_pair_follows_the_tangent_law(self): + # Loudspeakers at +-30 deg: (g1 - g2) / (g1 + g2) = tan(theta) / + # tan(30 deg), with g1^2 + g2^2 = 1 (Pulkki 1997, eq. 4 and 6). + ls = [30, -30] + groups = vbap.loudspeaker_pairs(ls) + inv = vbap.inverse_matrices(ls, groups) + theta = np.linspace(-29, 29, 59) + g = vbap.gains(theta, groups, inv) + law = (g[:, 0] - g[:, 1]) / (g[:, 0] + g[:, 1]) + err = np.max(np.abs(law - np.tan(np.radians(theta)) / np.tan(np.radians(30)))) + print(f'\n max |tangent law error| = {err:.1e}, ' + f'max |norm - 1| = {np.max(np.abs(np.linalg.norm(g, axis=1) - 1)):.1e}') + self.assertLess(err, 1e-12) + np.testing.assert_allclose(np.linalg.norm(g, axis=1), 1, atol=1e-12) + + def test_3d_gains_point_at_the_source(self): + # VBAP definition: the gains of the active triplet weight the + # loudspeaker vectors into a vector along the source direction, all + # gains are non-negative and their norm is 1. A source on a + # loudspeaker gets that loudspeaker alone. + ls = np.array([[0, 0], [90, 0], [180, 0], [-90, 0], [45, 45], + [135, 45], [-135, 45], [-45, 45], [0, 90]]) + groups = vbap.loudspeaker_triplets(ls) + inv = vbap.inverse_matrices(ls, groups) + rng = np.random.default_rng(1) + src = np.column_stack([rng.uniform(-180, 180, 200), rng.uniform(1, 89, 200)]) + g = vbap.gains(src, groups, inv) + L = np.array([unit(*d) for d in ls]) + p = g @ L + u = np.array([unit(*s) for s in src]) + cross = np.max(np.linalg.norm(np.cross(p, u), axis=1)) + print(f'\n max |p x u| = {cross:.1e} over 200 sources') + self.assertLess(cross, 1e-12) + self.assertTrue(np.all(g >= 0)) + np.testing.assert_allclose(np.linalg.norm(g, axis=1), 1, atol=1e-12) + on = vbap.gains(ls[[4]], groups, inv)[0] + self.assertAlmostEqual(on[4], 1, places=12) + self.assertAlmostEqual(np.sum(np.delete(on, 4)**2), 0, places=12) + + def test_matches_polarch(self): + # Same layouts and sources: pairs and triplets (1-based there), the + # inverse matrices (stacked by columns there) and the gains, with and + # without a spread of 20 deg (MDAP). + d = sio.loadmat(os.path.join(GOLDEN, 'vbap_polarch.mat'), squeeze_me=True) + pairs = vbap.loudspeaker_pairs(d['ls2']) + tri = vbap.loudspeaker_triplets(d['ls3']) + np.testing.assert_array_equal(pairs + 1, d['pairs']) + np.testing.assert_array_equal(tri + 1, d['tri']) + inv2 = vbap.inverse_matrices(d['ls2'], pairs) + inv3 = vbap.inverse_matrices(d['ls3'], tri) + diffs = { + 'inv2': np.abs(inv2.transpose(0, 2, 1).reshape(len(pairs), -1) - d['inv2']).max(), + 'inv3': np.abs(inv3.transpose(0, 2, 1).reshape(len(tri), -1) - d['inv3']).max(), + 'g2': np.abs(vbap.gains(d['src2'], pairs, inv2) - d['g2']).max(), + 'g2s': np.abs(vbap.gains(d['src2'], pairs, inv2, 20) - d['g2s']).max(), + 'g3': np.abs(vbap.gains(d['src3'], tri, inv3) - d['g3']).max(), + 'g3s': np.abs(vbap.gains(d['src3'], tri, inv3, 20) - d['g3s']).max()} + print('\n ' + ', '.join(f'{k} {v:.1e}' for k, v in diffs.items())) + for v in diffs.values(): + self.assertLess(v, 1e-12) + + +if __name__ == '__main__': + unittest.main()