From 46b7da82faa2619607a67987bd53325896d654c0 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 22:57:04 +0100 Subject: [PATCH] [FEATURE] add pytta.vbap with vector base amplitude panning pytta.vbap pans sources on 2D and 3D loudspeaker layouts (Pulkki 1997): loudspeaker_pairs and loudspeaker_triplets find the groups (the convex hull of the directions in 3D, without the faces that close a hemisphere), inverse_matrices inverts the matrix of each group and gains gives the gains of a set of source directions, with an optional spread (MDAP). It is a port of the MATLAB library Vector-Base-Amplitude-Panning of Archontis Politis, whose BSD 3-Clause licence is reproduced in the module header. Tests check the tangent law of a stereo pair, that the 3D gains weight the loudspeaker vectors into the source direction, and parity with the MATLAB library on golden data from tests/golden/make_vbap_polarch.m (groups identical, inverses and gains within 4.4e-16). Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/index.rst | 4 + docs/pyttadoc/vbap.rst | 7 + pytta/__init__.py | 5 + pytta/vbap.py | 254 +++++++++++++++++++++++++++++++ tests/golden/make_vbap_polarch.m | 20 +++ tests/golden/vbap_polarch.mat | Bin 0 -> 4942 bytes tests/test_vbap.py | 90 +++++++++++ 7 files changed, 380 insertions(+) create mode 100644 docs/pyttadoc/vbap.rst create mode 100644 pytta/vbap.py create mode 100644 tests/golden/make_vbap_polarch.m create mode 100644 tests/golden/vbap_polarch.mat create mode 100644 tests/test_vbap.py 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 0000000000000000000000000000000000000000..8a568a5180c3cbe940ca8bffad17c9682911e470 GIT binary patch literal 4942 zcma)9XHXN`wx&o85CsGj4G{~yNUxzOO3@&qpmdNTEfF}<2_Q&O1msBXsEAZS1cVTZ z(jn9X1PGvn5A8IK@BFYOGw|pzp;c%tAjiOk z5l0nTgV||yc5{f_u7^L%IbsFX*28PXJU8j*rkp+2uQ)AC6Q&>`ZH^H4Jt5ji`@`{5 zd}%v!0My3q+9>a(i+F|Y5lKHbZ<1YI^Slo__ybZoq;r8INNRXd8PR#nL6MqTs)xrU0OSi?RyYBD zv8gvK(>ZU2bc_Wbg7OsUrJJ9H@&a5L^Cxc2%#8&*vE8E7?{Tr@D8)qn7RtN-Qz)t? zU+-M~w@}8|vOgAZ26RE*f~4M%I5gji3#Zzh`$ zC+Y+2u{~7fFpCdI`Brr015S){4}Y+hp0BnFa-`_7uU%!pSTwC@<1lN|wsyYwS+B_~ zAqs-YU&7RVaZ$^hqhcA8`*9}r_fgHo2x4vIXlw|%sA6@nY!!w-C#S>QoUPifxmWyPr{Tl=4X`7~0Etq~xs=A^YPwWDPKTG#wj-q{JU2}U z$CkX9k+BNtTHpx^;h&PAS|Lsw(Q~;;n)2t`OIO$cQg$s}-KRuqy-HRV`1wQZ^n#6l zdX)6UI>1Hnb@^=TUZR?iCQ*lhpr1(;AzGn%$dD5v`VlO6G?KA6oJf9whc*9?@EtL=YRrt82Cwlzxj4$9vg+OSiLts|0ptN)o z;!VzGl!fOTdcY8H2&E&@5OR_qn%_Y!qc4nbF^R#84M6%>8{BcogjVC(Y`eY=X)>71 zVMwkX0cf_S3xQ}+( zFa#2x_vm@IbawNI1k{zTa06yyUoAnrb$k=v>t6BJInD-}W#^3r>qz7U@Rx3ea;Z52 z3@54l-d_`^y;zETMizdHRpG>cTrCIO_usCjZ8*&H?IoWj-)8U0r*Y?ck~nUBjEKG( z^}rumedM-p@v>b!V7_YbhA#f`+piZl7{lLzz@6AB{5RBoDtfkU{hC=GCew3EO51*7?+2^l-|-CDPaM&=|T zH{WRpJ;4u63|qsO_C}+t8u@km2o+?Z?xN=FFpPydNtq&Mv?yfUsn1o)z+IU@bzgce zRf$_M|4c|PGz@Dw_Ym|{=mBWWWXR-<ZOL`fwy<>YKgIN1{*g-GQI(kOC<=uwWqlofy!Bd&ymkG%n(abmY!&m<=T3+$ z>2*{gS$O*6&$=C8l!I)+$R5_|`|=9T;v96xKSZ&hKN=}m3vnmkp6Q$v-8mmWW16rO z?gQ)SaC-+s%hB`|B6b}kY5ffAKkcG9m948tp%jM zw?8p-S^I7=$Gf4?`R37}zx^rd z9*AK$*|W8Ncl7hwkuEDByKl+n4nz%kZ|nN*`d7U(72(0@QM!@8lH)bXrkyo9?zDP( zh>YWJob`t9K{L(kF6cxGCwS2`5cr--C(gVg`;tv5$AxhW1PTr#GSY!5Z_rDQib#@8wq{4vwE!~$>96Y^ClXW3IvfqdUye<-q}q;EOj zmFyAX{E{i={4Eb!cPOcTY?kW$uunAnh>NA^?m*{n(^%s8Z5lgSVU{+w!o9^M)_^$c zcOSUf*vzlVV&YC{-=Bz^dS~>w$*jgsAuC0iJ?X9}UM(xhR3?Ml``oqXH6_(?Mh~uA zkc2F zE3R@PU}1!_lvVYVJ6?M0fTAd8SbzJIC)Y|^y*@_qyhWfkzR6^s>mQMR7~CIrd5q|T z#^O#|)fPR)lykO!k*?)$6R{?di-&GglkyB3Br9drNgzu?5pDKqU-7a`pi04(Ust;(t_Ab{O$hHXdfKxft+0oPhXMh-HA4U^L{%}N{Dgtr{r4|v79V4u3y zcHGTF>rnX(9p^Jppc@RnzxpJ$8gt7bU3N^qYUcG4h&X9YSU(@imw|?uLndWt(?lkP4guq} z-;@`je@cL!f`J^upxY)OUTzStc4%ncGS!q9m?&Y8s=NF`Z{b^J;)DL4H@PrQ;L@Lm zr)l}AbjQ&uRN%=yBHDN`?AAdA>k8A3h**}}73UM`%S<%ULe*-a-k&?rWv4F3NA_Qb ze7JL%Zn;X`xaG|=k~FEl3zRu32VPUjtB4OqMYZT zpS{Ibj-Cq8Epye(^$zYTJ2iUj`80|N3r9_l;iulKXu_MZ?@%saE;s*!c8IuajA9C- z#iqr_+zMhDaMQx~<5!-Nlhk_=_l-15l_nhnf#HA)6Ve$S_3w@lxFDiN2`MNZCZ?D~cDMzzotZ^$%U&Ojq^R1xl{Y3g7y&pM(rxaVD zVl^?^@8^y6+|N2~%pEs+Efp+?NI%0wK(*otFK%46&tC~7$&&h4+RTsmpR=~L zk6!5%9VmlZFT!GO!)FI9+a5~ylVLOu>s0Hg*?J*dg#%+*%@VGV?QzACHS;6RUxP98 zf#l^rUu0cZm0jLt7j?`wQnfkOHEB@V>w#+=-M&vSxk{EhNBbjnW9Fm>{U)*J8u8i1 zquBRK#wZ07)b;WuPHjSSI)u`A>ZsNEd)S?OOPbvzq#g;k!BtC&;z5n$A3ux!%hVvh ztk}%pIOd=sD<|7ROX)2`&jU0B2bm!`J+L)!ti-ACv(Wf{A#3Ls?U>tpc-__0l}de7 z!Jm->@Q!#b=STXub6Wl*tl8MrQn4_S0wKi#^=FGj)6f!IK$VZT#JaeCj`xF35Vi#Qs?j{KCrW!noDM>r*bQf2vuY+7!7NBQs(A z<`s1LIBxrd>n71{Xnmcc=q`zAN)Ep)(En5SYeCwesyy_C&gX})I0~g~({gC|YEG=YuvZ6W0+0r0xZhhw6b9YS(mL>B=dA$k7UiPCkhWGSvV z+tSTseZ==$n=dmJ%P}+WBs$o_m)%Dwp2Z!JAp>4-4X!P-L$8qR(-Njhfc46aA%#e* z-ABZ>kKEH(b6Y=9H<#pwOR>DV=?=4}-3PeC6*)t(UNONX60XkxBJMK`i-|vJMM~dL z76HF6LFJ!EfwQ=2fY)Lsz@5A}RpYsre8YU&-%Sl)$HoGxjbE*dGxeUEF)!d*-6}b< zej|kdGKlS}&`DJlWz|(-adtCK&eQcN($5w6h8SBapH*!bsN-T*f#6bQ-eGnY zlO>+C-$YmcA=>*dqG#S$z@G69072CkDWz@>;H#ddS42*FCTy^?cnY(P{syA{1yoPd z7l<=FnUH-gRRa_YXA@G8ANUs69Vh0NbRiz072^~$@%m!u33H2Aou~9t`Xd8Ws863_ zBih>3d)pPAJvj4 zlHUl*H0C$(4W^(7moUllD;Swx(c0v=}b@ApNr(oo>6KB$6xCN&S<1}t6v+%qJcky1D zeRe_=+y!H{imi&h6k8B0dwyD+1&1gR`knLUzj7+PwBU^~#qBrQ6GwODnGiW#%B#=m=#-eXIL>V<3qJ4<%r7qqJ?ehe*hX# B%ozXx literal 0 HcmV?d00001 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()