Acoustic Holography and Holophony
Franz Zotter, 2016.
The examples below show the directions and widths indicated by the $\boldsymbol{r}_\mathrm{V}$ and $\boldsymbol{r}_\mathrm{E}$ measures for pairwise vector-base amplitude panning along the horizon.
Assume that you are given the desired panning direction as
\begin{equation} \boldsymbol\theta=\begin{bmatrix} \cos\varphi_\mathrm{s}\\ \sin\varphi_\mathrm{s} \end{bmatrix} \end{equation}and the coordinates of your loudspeakers as
\begin{equation} \mathbf{X}=\begin{bmatrix} \cos\varphi_1, &\dots, &\cos\varphi_\mathrm{L}\\ \sin\varphi_1, &\dots, &\sin\varphi_\mathrm{L} \end{bmatrix}. \end{equation}Then VBAP searches the loudspeaker pair that encloses the direction $\boldsymbol\theta$. This is done by constructing the convex hull, which defines the lines (simplices) connecting all the loudspeaker pairs by the loudspeaker indices involved. For a ring of $\mathrm{L}$ loudspeakers with loudspeakers sorted in azimuth, it looks like this
\begin{equation} Q=\mathrm{convhull}\{\mathbf{X}\}=\begin{bmatrix} 1& 2\\ 2& 3\\ \vdots &\vdots\\ a& b\\ \vdots &\vdots\\ \mathrm{L}-1& \mathrm{L}\\ \mathrm{L} & 1 \end{bmatrix}. \end{equation}Now VBAP searches through the convex hull to find the loudspeaker pair (simplex of the indices $a,b$) of which the wheighted superposition of direction vectors allows for a solution with numerically all-positive weights $\boldsymbol{g}_q=[g_a,g_b]$ with $g_a,g_b>0$ to represent $\boldsymbol\theta$:
$\boldsymbol g=\boldsymbol 0$ with length $\mathrm{L}$
For all $q$ indexing through all pairs:
$\boldsymbol{g}_q=\mathbf{X}[Q[q,:]]^{-1}\boldsymbol{\theta}$
if $g_a,g_b>-10^{-3}$:
$\boldsymbol{g}[[a,b],:]=\boldsymbol{g}_q$
break loop
This algorithm is defined below.
In [1]:
import numpy as np
from numpy.linalg import inv
def vectorpan(xys,xyls,simplices):
g=np.zeros(phils.shape[0])
for n in range(0,simplices.shape[0]):
na=simplices[n,0]
nb=simplices[n,1]
M=np.array([xyls[:,na],xyls[:,nb]]).T
gnm=np.dot(inv(M),xys)
if np.sum(gnm<-1e-3)==0:
g[na]=gnm[0]
g[nb]=gnm[1]
break
return g
The vector measures are easily extended to more than 2 loudspeakers, i.e., to
\begin{equation} \boldsymbol{r}_\mathrm{V}= \frac{\sum_{l=1}^\mathrm{L}g_l\,\boldsymbol\theta_l}{\sum_{l=1}^\mathrm{L}g_l}, \qquad\qquad \boldsymbol{r}_\mathrm{E}= \frac{\sum_{l=1}^\mathrm{L}g_l^2\,\boldsymbol\theta_l}{\sum_{l=1}^\mathrm{L}g_l^2}. \end{equation}
In [2]:
def r_vector(g,theta):
L=theta.size
thx=np.sin(theta[:])
thy=np.cos(theta[:])
rx=np.dot(thx,g).T
ry=np.dot(thy,g).T
normalizer=sum(g,0)
rx/=normalizer
ry/=normalizer
return np.array([rx,ry])
This vector model is now applied on an example using 5 equally spaced loudpeaker on the horizon and panning for all directions from $-180^\circ\dots180^\circ$. It shows again the directions and widths displayed by $\boldsymbol r_\mathrm{V,E}$.
In [3]:
from scipy.spatial import ConvexHull
from bokeh.plotting import figure, output_file, show
from bokeh.io import output_notebook
L=5
Npts=200
phils=np.arange(0,L)*2*np.pi/L
xyls=np.array([np.cos(phils),np.sin(phils)])
phis=np.linspace(-np.pi*0.99,np.pi,Npts)
xys=np.array([np.cos(phis),np.sin(phis)])
g=np.zeros([L,Npts])
qh=ConvexHull(xyls.T)
for n in range(0,Npts):
gn=vectorpan(xys[:,n],xyls,qh.simplices)
gn=np.abs(gn)
gn/=np.sqrt(np.sum(gn))
g[:,n]=gn
output_notebook()
p1=figure(title="r vector direction VBAP",plot_width=300,plot_height=250)
p2=figure(title="r vector width VBAP",plot_width=300,plot_height=250,x_range=[-180,180],y_range=[0,100])
r=r_vector(g,phils)
mr=np.sqrt(np.sum(r**2,0))
p1.line(phis*180/np.pi,np.arctan2(r[0,:],r[1,:])*180/np.pi,legend="rV",line_width=2)
p2.line(phis*180/np.pi,2*np.arccos(mr)*180/np.pi,legend="2acos||rV||",line_width=2)
r=r_vector(g**2,phils)
mr=np.sqrt(np.sum(r**2,0))
p1.line(phis*180/np.pi,np.arctan2(r[0,:],r[1,:])*180/np.pi,legend="rE",color="red",line_width=2)
p2.line(phis*180/np.pi,2*np.arccos(mr)*180/np.pi,legend="2acos||rE||",color="red",line_width=2)
p1.legend.location="top_left"
p2.legend.background_fill_alpha = 0.5
show(p1)
show(p2)
It seems that the vectors do not indicate the same directions: only the $\boldsymbol r_\mathrm{V}$ vector is perfectly controlled. Alternatively to VBAP, VBIP could be done, which just takes the square roots of the panning gains obtained from VBAP, before normalization.
In [4]:
for n in range(0,Npts):
gn=vectorpan(xys[:,n],xyls,qh.simplices)
gn=np.abs(gn)
# new line inserted just for VBIP:
gn=np.sqrt(gn)
gn/=np.sqrt(np.sum(gn))
g[:,n]=gn
p3=figure(title="r vector direction VBIP",plot_width=300,plot_height=250)
p4=figure(title="r vector width VBIP",plot_width=300,plot_height=250,x_range=[-180,180],y_range=[0,100])
r=r_vector(g,phils)
mr=np.sqrt(np.sum(r**2,0))
p3.line(phis*180/np.pi,np.arctan2(r[0,:],r[1,:])*180/np.pi,legend="rV",line_width=2)
p4.line(phis*180/np.pi,2*np.arccos(mr)*180/np.pi,legend="2acos||rV||",line_width=2)
r=r_vector(g**2,phils)
mr=np.sqrt(np.sum(r**2,0))
p3.line(phis*180/np.pi,np.arctan2(r[0,:],r[1,:])*180/np.pi,legend="rE",color="red",line_width=2)
p4.line(phis*180/np.pi,2*np.arccos(mr)*180/np.pi,legend="2acos||rE||",color="red",line_width=2)
p3.legend.location="top_left"
p4.legend.background_fill_alpha = 0.5
show(p3)
show(p4)
Now the $\boldsymbol r_\mathrm{E}$ vector is perfectly controlled and the widths tend to have narrower notches.
MDAP does the same as VBAP, but it always superimposes several sources to obtain the panning gains. One could describe it as
\begin{equation} \boldsymbol{\tilde g}=\mathrm{VBAP}\{\varphi_\mathrm{s}-\textstyle\frac{\alpha}{2},\mathbf{X}\}+\mathrm{VBAP}\{\varphi_\mathrm{s},\mathbf{X}\}+\mathrm{VBAP}\{\varphi_\mathrm{s}+\textstyle\frac{\alpha}{2},\mathbf{X}\}\qquad\qquad \boldsymbol g=\frac{\boldsymbol{\tilde g}}{\|\boldsymbol{\tilde g}\|} \end{equation}
In [5]:
alpha=phils[1]-phils[0]
xys2=np.array([np.cos(phis-alpha/2),np.sin(phis-alpha/2)])
xys3=np.array([np.cos(phis+alpha/2),np.sin(phis+alpha/2)])
for n in range(0,Npts):
gn=vectorpan(xys[:,n],xyls,qh.simplices)
# lines inserted for MDAP
gn+=vectorpan(xys2[:,n],xyls,qh.simplices)
gn+=vectorpan(xys3[:,n],xyls,qh.simplices)
gn=np.abs(gn)
gn/=np.sqrt(np.sum(gn))
g[:,n]=gn
p5=figure(title="r vector direction MDAP",plot_width=300,plot_height=250)
p6=figure(title="r vector width MDAP",plot_width=300,plot_height=250,x_range=[-180,180],y_range=[0,100])
r=r_vector(g,phils)
mr=np.sqrt(np.sum(r**2,0))
p5.line(phis*180/np.pi,np.arctan2(r[0,:],r[1,:])*180/np.pi,legend="rV",line_width=2)
p6.line(phis*180/np.pi,2*np.arccos(mr)*180/np.pi,legend="2acos||rV||",line_width=2)
r=r_vector(g**2,phils)
mr=np.sqrt(np.sum(r**2,0))
p5.line(phis*180/np.pi,np.arctan2(r[0,:],r[1,:])*180/np.pi,legend="rE",color="red",line_width=2)
p6.line(phis*180/np.pi,2*np.arccos(mr)*180/np.pi,legend="2acos||rE||",color="red",line_width=2)
p5.legend.location="top_left"
p6.legend.background_fill_alpha = 0.5
show(p5)
show(p6)
Obviously MDAP could smooth out much of the variations in $\boldsymbol{r}_\mathrm{V,E}$ and also the difference between both of the measures.
In [ ]: