# Binaural rendering from microphone array signals of arbitrary geometry

Naoto Iijima  Shoichi Koyama  Hiroshi Saruwatari  The University of Tokyo,
Graduate School of Information Science and Technology, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-8656, Japan

###### Abstract

A method of binaural rendering from microphone array signals of arbitrary
geometry is proposed. To reproduce binaural signals from microphone array
recordings at a remote location, a spherical microphone array is generally
used for capturing a soundfield. However, owing to the lack of flexibility in
the microphone arrangement, the single spherical array is sometimes
impractical for estimating a large region of a soundfield. We propose a method
based on harmonic analysis of infinite order, which allows the use of
arbitrarily placed microphones. In the synthesis of the estimated soundfield,
a spherical-wave-decomposition-based binaural rendering is also formulated to
take into consideration the distance in measuring head-related transfer
functions. We develop and evaluate a composite microphone array consisting of
multiple small arrays. Experimental results including those of listening tests
indicate that our proposed method is robust against change in listening
position in the recording area.

##  I  Introduction

Spatial sound reproduction is essential to achieve virtual or augmented
reality (VR/AR) audio systems; therefore, a wide range of spatial audio
techniques have been developed. Soundfield reproduction/control using multiple
loudspeakers makes it possible to reproduce a spatial sound inside a
relatively large region Berkhout:JASA_J_1993 ; Daniel:AES114conv ;
Poletti:J_AES_2005 ; Spors:AES124conv ; Koyama:IEEE_J_ASLP2013 . On the other
hand, binaural rendering techniques reproduce sounds received in both the
listener’s left and right ears generally through headphones. Binaural signals
are typically synthesized using transfer functions from sources to the
eardrums in a free field, i.e., the so-called head-related transfer functions
(HRTFs), measured in advance. An HRTF includes the effects of diffraction and
reflection caused by the listener’s head, torso, and external pinna shape,
which contains important features for sound localization in humans HRTF .

When virtual source positions are given, it is possible to synthesize the
binaural signals in a free field by the direct convolution of transfer
functions and source signals; however, it is not easy to synthesize binaural
signals from signals captured by microphones in a reverberant environment,
where the synthesized signals should correspond to the signals received in
both ears when the listener is present in the recording area. It is also
preferable that the binaural signals vary in response to the listener’s head
movement from the perspective of the sound localization property in humans
headmove . One of the solutions will be the use of a dummy head that moves
synchronously with the listener’s movement telehead . However, it is difficult
to simultaneously bring the binaural signals to multiple listeners. One of the
promising approaches is binaural reproduction from microphone array recordings
senzi ; vah ; vah2 ; beam ; beam2 ; beam3 ; conv2 . Owing to the recent
developments of spatial soundfield recording techniques, it is becoming
possible to estimate a soundfield inside a region using multiple microphones
Park:JASA_J_2005 ; Poletti:J_AES_2005 ; micarray ; samarashinge_2014_ASLP ;
Ueno:IEEE_SPL2018 ; Nakanishi:WASPAA2019 . The reproduced binaural signals can
be highly accurate in a physical sense at low frequencies owing to the use of
the spatial soundfield information estimated by these techniques.

A typical microphone array used in the soundfield recording is a spherical
microphone array mounted on an acoustically rigid object conv2  since the
spherical harmonic analysis of a soundfield is effective for binaural
rendering conv ; plane_trans ; schultz2013data-based ; hrtf_pw_modal ;
hrtf_sph_prepro2 . The number of microphones and array size limit the
reproducible frequency range and size of the reproducing region. In
particular, for the binaural reproduction, the region of high estimation
accuracy should be larger than the head size. Therefore, the use of a single
spherical microphone array is sometimes impractical owing to the lack of
flexibility in the microphone arrangement.

We propose a method of binaural reproduction from distributed microphone array
signals. The soundfield estimation method based on harmonic analysis of
infinite order is applied to obtain harmonic coefficients of the soundfield in
the recording area. This estimation method makes it possible to use
arbitrarily placed microphones for capturing a sound field. For example,
multiple small microphone arrays, such as low-order spherical microphones
(ambisonics microphones) and acoustic vector sensors, will be more practically
feasible than a single large spherical microphone array. Moreover, such a
composite array will have flexibility and scalability in the array placement.
In the binaural rendering from the estimated harmonic coefficients, we use the
soundfield representation by spherical wave decomposition Iijima:IEEE_MMSP2020
. In contrast to a typically used plane wave decomposition method, our
rendering method takes into consideration the distance in HRTF measurement. In
numerical simulations, it is shown that the binaural signals can be accurately
reproduced by using our soundfield estimation method with the composite array
of microphones. Moreover, the peaks and dips in the amplitude response of the
binaural signals reproduced on the basis of the conventional plane wave
decomposition method can be shifted from the true ones. They are accurately
reproduced by our spherical-wave-decomposition-based binaural rendering
method. We also develop a practical microphone array consisting of multiple
small microphone arrays, and the binaural signals reproduced by using this
array system are also evaluated.

The rest of this paper is organized as follows. In Sect. II, we introduce
several notations and basic theories on sound field representation. The
problem statement and prior works are introduced in Sect. III. In Sect. IV, we
describe our proposed binaural rendering method for distributed microphone
array signals. Numerical simulation results for comparing binaural rendering
and soundfield estimation methods are shown in Sects. V and VI, respectively.
Experimental results in a practical environment are shown in Sect. VII.
Finally, our conclusion is presented in Sect. VIII.

##  II  Preliminaries

First, we introduce several notations and basic theories used throughout this
paper.

###  II.1  Notations

Italic letters denote scalars, lower case italic boldface letters denote
vectors, and upper case italic boldface letters denote matrices. The sets of
real and complex numbers are denoted by $\mathbb{R}$ and $\mathbb{C}$,
respectively. The unit sphere in $\mathbb{R}^{3}$ is denoted by
$\mathbb{S}_{2}$. Subscripts of scalars, vectors, and matrices indicate their
indexes. For example, $x_{i,j}$ is the $(i,j)$th entry of matrix $\bm{X}$.

The imaginary unit is denoted as $\mathrm{j}=\sqrt{-1}$. The complex
conjugate, transpose, conjugate transpose, and inverse are denoted by
superscripts $(\cdot)^{\ast}$, $(\cdot)^{\mathsf{T}}$, $(\cdot)^{\mathsf{H}}$,
and $(\cdot)^{-1}$, respectively. The absolute value of a scalar $x$ and the
Euclid norm of a vector $\bm{x}$ are denoted as $|x|$ and $\|\bm{x}\|$,
respectively. The inner product of vectors $\bm{x}$ and $\bm{y}$ is denoted as
$\langle\bm{x},\bm{y}\rangle$.

The functions $j_{n}(\cdot)$ and $h_{n}(\cdot)$ are the $n$th order spherical
Bessel function and $n$th-order spherical Hankel function of the second kind,
respectively. The spherical harmonic function of order $n$ and degree $m$ for
azimuth and zenith angles ($\phi$ and $\theta$, respectively) is defined as

$\displaystyle
Y_{n}^{m}(\theta,\phi)=\sqrt{\frac{2n+1}{4\pi}\frac{(n-m)!}{(n+m)!}}P_{n}^{m}(\cos\theta)e^{\mathrm{j}m\phi},$
(1)

where $P_{n}^{m}(\cdot)$ is the associated Legendre function of order $n$ and
degree $m$.

###  II.2  Sound field representation using spherical wavefunctions

The stationary soundfield of angular frequency $\omega$ inside a region not
including any sources $D$ satisfies the homogeneous Helmholtz equation as

$\displaystyle(\Delta+k^{2})u(\bm{r},k)=0,$ (2)

where $\Delta$ is the Laplacian and $u(\bm{r},k)$ is the sound pressure at the
position $\bm{r}\in D$ and wavenumber $k$ defined as $k=\omega/v$ with sound
velocity $v$. The solution of Eq. (2) for the source-free interior region $D$
can be expanded using spherical wavefunctions at the expansion center
$\bm{r}_{\mathrm{c}}$ as

$\displaystyle
u(\bm{r},k)=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\alpha_{n}^{m}(\bm{r}_{\mathrm{c}},k)\varphi_{n}^{m}(\bm{r}-\bm{r}_{\mathrm{c}},k).$
(3)

Here, $\varphi_{n}^{m}(\cdot)$ is the spherical wavefunction defined for
$\bm{r}=[r,\theta,\phi]^{\mathsf{T}}$ in the spherical coordinates as

$\displaystyle\varphi_{n}^{m}(\bm{r},k)=\sqrt{4\pi}j_{n}(kr)Y_{n}^{m}(\theta,\phi).$
(4)

We also represent Eq. $\eqref{eq:shd}$ in a matrix form as

$\displaystyle
u(\bm{r},k)=\bm{\alpha}(\bm{r}_{\mathrm{c}},k)^{\mathsf{T}}\bm{\varphi}(\bm{r}-\bm{r}_{\mathrm{c}},k),$
(5)

where $\bm{\alpha}(\bm{r}_{\mathrm{c}},k)\in\mathbb{C}^{\infty}$ and
$\bm{\varphi}(\bm{r}-\bm{r}_{\mathrm{c}},k)\in\mathbb{C}^{\infty}$ are the
infinite dimensional vectors of the expansion coefficients
$\alpha_{n}^{m}(\bm{r}_{\mathrm{c}},k)$ and spherical wavefunctions
$\varphi_{n}^{m}(\bm{r}-\bm{r}_{\mathrm{c}},k)$, respectively. Note that
$\alpha_{0}^{0}(\bm{r},k)$ corresponds to $u(\bm{r},k)$.

The expansion coefficients at two expansion centers $\bm{r}$ and
$\bm{r}^{\prime}$ ($\bm{\alpha}(\bm{r},k)$ and
$\bm{\alpha}(\bm{r}^{\prime},k)$, respectively) can be related as

$\displaystyle\bm{\alpha}(\bm{r},k)=\bm{T}(\bm{r}-\bm{r}^{\prime},k)\bm{\alpha}(\bm{r}^{\prime},k),$
(6)

where $\bm{T}(\bm{r}-\bm{r}^{\prime},k)\in\mathbb{C}^{\infty\times\infty}$ is
the translation operator defined as

$\displaystyle\left[\bm{T}(\bm{r})\bm{\alpha}\right]_{n^{\prime}}^{m^{\prime}}$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}[4\pi(-1)^{m}\mathrm{j}^{n^{\prime}-n}$
$\displaystyle\cdot\sum_{l=0}^{n^{\prime}+n}\mathrm{j}^{l}j_{l}(kr)Y_{l}^{m^{\prime}-m}(\theta,\phi)^{*}\mathcal{G}(n,m;n^{\prime},-m^{\prime};l)]\alpha_{n}^{m}.$
(7)

Here, $[\cdot]_{n^{\prime}}^{m^{\prime}}$ is the element of order $n^{\prime}$
and degree $m^{\prime}$, and $\mathcal{G}(\cdot)$ is the Gaunt coefficient
multiple . The translation operator satisfies the following equations:

$\displaystyle\bm{T}(-\bm{r},k)$
$\displaystyle=\bm{T}^{-1}(\bm{r},k)=\bm{T}^{\mathsf{H}}(\bm{r},k),$ (8)
$\displaystyle\bm{T}(\bm{r}+\bm{r}^{\prime},k)$
$\displaystyle=\bm{T}(\bm{r},k)\bm{T}(\bm{r}^{\prime},k).$ (9)

###  II.3  Head-related transfer function and its representation in spherical
harmonic domain

Frequency-domain transfer functions from sources to the listener’s left and
right eardrums are usually measured using loudspeakers on a spherical surface
in an anechoic chamber datahrtf2 . We refer to these transfer functions as
HRTFs although an HRTF is generally defined as that divided by the transfer
function from the loudspeaker to the center of the head without considering
the listener. As shown in Fig. 1, we assume that HRTFs are measured from
loudspeakers on a spherical surface $\partial\Omega$, whose radius is denoted
as $R_{\mathrm{s}}$, in a free field. The HRTFs for the left and right
eardrums at the wavenumber $k$ are denoted as
$h_{\mathrm{L}}(\bm{r}_{\mathrm{s}},k)$ and
$h_{\mathrm{R}}(\bm{r}_{\mathrm{s}},k)$, respectively, with the position of
the loudspeaker $\bm{r}_{\mathrm{s}}\in\partial\Omega$. Hereafter, the
subscripts $\mathrm{L}$ and $\mathrm{R}$ represent the signals for the left
and right ears including HRTFs, but they are described together as
$\mathrm{L,R}$ in the subscript to represent both the left and right ear
signals, e.g., $h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}},k)$.

The HRTF representation in the spherical harmonic domain is generally used for
interpolating binaural signals god ; god2 ; god3 . We use this representation
to obtain a simple formulation of binaural signals on the basis of the
harmonic analysis of the soundfield. The HRTFs are measured by loudspeakers at
$J$ positions on $\partial\Omega$, whose positions are denoted by
$\bm{r}_{\mathrm{s},j}=[R_{\mathrm{s}},\theta_{\mathrm{s},j},\phi_{\mathrm{s},j}]^{\mathsf{T}}$
($j\in\\{1,\ldots,J\\}$). We consider to expand
$h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}},k)$ using the complex conjugate of
spherical harmonic functions as

$\displaystyle
h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s},j},k)=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}H_{\mathrm{L,R},n}^{m}(k)Y_{n}^{m}(\theta_{\mathrm{s},j},\phi_{\mathrm{s},j})^{\ast}.$
(10)

By truncating the infinite sum up to $N$, we can describe Eq. (10) using the
linear equation

$\displaystyle\bm{h}_{\mathrm{L,R}}(k)=\bm{Y}\bm{\psi}_{\mathrm{L,R}}(k),$
(11)

where $\bm{h}_{\mathrm{L,R}}(k)\in\mathbb{C}^{J}$ and
$\bm{\psi}_{\mathrm{L,R}}(k)\in\mathbb{C}^{(N+1)^{2}}$ are the vectors
consisting of $h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s},j},k)$ and
$H_{\mathrm{L,R},n}^{m}(k)$, respectively, and
$\bm{Y}\in\mathbb{C}^{J\times(N+1)^{2}}$ is the matrix consisting of spherical
harmonic functions
$Y_{n}^{m}(\theta_{\mathrm{s},j},\phi_{\mathrm{s},j})^{\ast}$. When
$J\geq(N+1)^{2}$, $\bm{\psi}_{\mathrm{L,R}}(k)$ can be obtained as the least-
squares solution with Tikhonov regularization as

$\displaystyle\bm{\psi}_{\mathrm{L,R}}(k)=\left(\bm{Y}^{\mathsf{H}}\bm{Y}+\gamma\bm{Q}\right)^{-1}\bm{Y}^{\mathsf{H}}\bm{h}_{\mathrm{L,R}}(k),$
(12)

where $\gamma$ is a constant parameter, $\bm{Q}$ is the weighting diagonal
matrix, and the diagonal element for the $n$th-order coefficient of $\bm{Q}$
is $1+n(n+1)$ god ; god3 .

Figure 1: Geometry for HRTF measurement in spherical coordinates. The transfer
functions from the loudspeaker on a spherical surface $\partial\Omega$ to the
listener’s left and right eardrums are obtained.

##  III Problem statement and prior works

###  III.1 Binaural rendering from microphone array signals

Our objective is to reproduce binaural signals received by a listener at a
remote location by using measurements of multiple microphones placed in a
recording area. By using the listener’s HRTFs (or their proxies), the binaural
signals are synthesized so that they are identical to those received when the
listener is present in the recording area. The reproduced binaural signals
should vary depending on the listener’s head position and direction. To
reproduce the spatial sound including reverberations, it is necessary to
capture the soundfield inside a target listening region by using multiple
microphones. Then, the binaural signals are synthesized using HRTFs.

Suppose that a source-free region of interest $D$ is set in the recording area
and $I$ microphones are placed inside $D$. The microphone signals of
wavenumber $k$ are denoted as
$\bm{s}(k)=[s_{1}(k),\ldots,s_{I}(k)]^{\mathsf{T}}$ ($\in\mathbb{C}^{I}$). The
HRTFs from $J$ loudspeakers on a spherical surface $\partial\Omega$ with
respect to the listener’s left and right eardrums,
$h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}},k)$
($\bm{r}_{\mathrm{s}}\in\partial\Omega$), are measured beforehand. The goal is
to reproduce the binaural signals of the listener,
$y_{\mathrm{L,R}}(\bm{r},k)$, when the listener (the origin of the HRTF
measurement, i.e., the listener’s head center) is positioned at $\bm{r}\in D$.

This binaural reproduction problem is twofold: soundfield estimation by
multiple microphones and binaural rendering using the estimated soundfield. We
introduce some prior works in the following sections. Hereafter, the
wavenumber $k$ is omitted for notational simplicity.

###  III.2 Estimation of expansion coefficients using spherical microphone
array

Generally, expansion coefficients of the soundfield with spherical
wavefunctions are estimated in the recording area. A single spherical
microphone array is typically used for such estimation conv2 . To avoid the
forbidden frequency problem, microphones mounted on an acoustically rigid
baffle, an array of directional microphones, or multiple layers of microphone
arrays are usually employed Poletti:J_AES_2005 ; Rafaely:IEEE_J_ASLP_2005 ;
dualsphere ; Koyama:JASA_J_2016 . We here assume that the spherical microphone
array with a rigid baffle is placed in $D$ with its center at the origin
$\bm{r}_{\mathrm{o}}$.

Since the spherical baffle is assumed to be acoustically rigid, the radial
velocity on the surface becomes zero. Then, by denoting the radius of the
spherical array as $R_{\mathrm{m}}$, the sound pressure $u(\bm{r})$
($\bm{r}=[r,\theta,\phi]^{\mathsf{T}}\in D$, $r\geq R_{\mathrm{m}}$) is
described by spherical wavefunction expansion as

$\displaystyle u(\bm{r})$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})$
$\displaystyle\ \
\cdot\sqrt{4\pi}\left(j_{n}(kr)-\frac{j_{n}^{\prime}(kR_{\mathrm{m}})}{h^{\prime}_{n}(kR_{\mathrm{m}})}h_{n}(kr)\right)Y_{n}^{m}(\theta,\phi),$
(13)

where $\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})$ is the interior expansion
coefficient at the array center $\bm{r}_{\mathrm{o}}$. Thus, the observation
of the $i$th microphone $s_{i}$ at
$\bm{r}_{\mathrm{m},i}=[R_{\mathrm{m}},\theta_{\mathrm{m},i},\phi_{\mathrm{m},i}]^{\mathsf{T}}$
is represented as

$\displaystyle
s_{i}=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\frac{-\sqrt{4\pi}\mathrm{j}}{k^{2}R_{\mathrm{m}}^{2}h_{n}^{\prime}(kR_{\mathrm{m}})}Y_{n}^{m}(\theta_{\mathrm{m},i},\phi_{\mathrm{m},i})\alpha_{n}^{m}(\bm{r}_{\mathrm{o}}).$
(14)

Here, the following relation is used:

$\displaystyle
j_{n}(x)h_{n}^{\prime}(x)-j_{n}^{\prime}(x)h_{n}(x)=-\frac{\mathrm{j}}{x^{2}}.$
(15)

One approach to estimating $\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})$ up to a
predefined truncation order is to derive an analytical representation of
$\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})$ from Eq. (14) conv2 , where the integral
over the spherical microphone surface must be computed with the discrete
observation values $\bm{s}$. We here introduce the least-squares solution of
Eq. (14) to avoid numerical integration. By truncating the infinite sum up to
$N$, Eq. (14) can be approximated using the following linear equation:

$\displaystyle\bm{s}=\bm{\Pi}\bar{\bm{\alpha}}(\bm{r}_{\mathrm{o}}),$ (16)

where $\bar{\bm{\alpha}}(\bm{r}_{\mathrm{o}})\in\mathbb{C}^{(N+1)^{2}}$ is the
vector consisting of $\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})$ and
$\bm{\Pi}\in\mathbb{C}^{I\times(N+1)^{2}}$ is the matrix consisting of the
coefficients except $\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})$ in Eq. (14). When
$I\geq(N+1)^{2}$, $\bar{\bm{\alpha}}(\bm{r}_{\mathrm{o}})$ can be obtained as
the least-squares solution with Tikhonov regularization as

$\displaystyle\bar{\bm{\alpha}}(\bm{r}_{\mathrm{o}})=(\bm{\Pi}^{\mathsf{H}}\bm{\Pi}+\eta\bm{I})^{-1}\bm{\Pi}^{\mathsf{H}}\bm{s},$
(17)

where $\eta$ is a constant parameter. This method is applicable only to the
estimation of the expansion coefficients at the array center.

###  III.3 Binaural rendering based on plane wave decomposition

The binaural signals are reproduced from the expansion coefficients estimated
in the recording area. Many binaural rendering methods are based on the plane
wave decomposition of the soundfield conv ; conv2 ; schultz2013data-based .
First, the soundfield $u(\bm{r})$ is represented by the weighted integral of
plane waves as

$\displaystyle
u(\bm{r})=\frac{1}{4\pi}\int_{\mathbb{S}_{2}}w_{\mathrm{pw}}(\bm{\eta})e^{\mathrm{j}k\langle\bm{\eta},\bm{r}\rangle}\mathrm{d}\bm{\eta},$
(18)

where $w_{\mathrm{pw}}(\bm{\eta})$ is the weight for the plane wave of the
arrival direction
$\bm{\eta}=[1,\theta^{\prime},\phi^{\prime}]^{\mathsf{T}}\in\mathbb{S}_{2}$.
This representation is referred to as the Herglotz wavefunction singlelayer .
We here reformulate the plane wave expansion using expansion coefficients of
the spherical wavefunctions. The plane wave function can be expanded as

$\displaystyle
e^{\mathrm{j}k\langle\bm{\eta},\bm{r}\rangle}=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}4\pi\mathrm{j}^{n}j_{n}(kr)Y_{n}^{m}(\theta^{\prime},\phi^{\prime})^{\ast}Y_{n}^{m}(\theta,\phi).$
(19)

By substituting Eq. (19) into Eq. (18), one can obtain

$\displaystyle u(\bm{r})=$ $\displaystyle\ \
\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\mathrm{j}^{n}j_{n}(kr)Y_{n}^{m}(\theta,\phi)\int_{\mathbb{S}_{2}}w_{\mathrm{pw}}(\bm{\eta})Y_{n}^{m}(\theta^{\prime},\phi^{\prime})^{\ast}\mathrm{d}\bm{\eta}.$
(20)

The plane wave weight $w_{\mathrm{pw}}(\bm{\eta})$ can be related to the
expansion coefficients $\alpha_{n}^{m}(\bm{r})$ in the recording area as

$\displaystyle
w_{\mathrm{pw}}(\bm{\eta})=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\sqrt{4\pi}\mathrm{j}^{-n}\alpha_{n}^{m}(\bm{r})Y_{n}^{m}(\theta^{\prime},\phi^{\prime}).$
(21)

This relation can be confirmed by substituting Eq. (21) into Eq. (20) and by
using the orthogonality of the spherical harmonic functions.

The binaural signals are usually obtained by the summation of
$h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})$ multiplied by the weight
$w_{\mathrm{pw}}(\bm{\eta})$ of each angle. We here obtain a simple
formulation by regarding $h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})$ as functions
of continuous $\bm{r}_{\mathrm{s}}\in\partial\Omega$. We represent the
binaural signals at the position $\bm{r}_{\mathrm{o}}$ as

$\displaystyle
y_{\mathrm{L,R}}(\bm{r}_{\mathrm{o}})=\frac{1}{R_{\mathrm{s}}^{2}}\int_{\partial\Omega}w_{\mathrm{pw}}(\hat{\bm{r}}_{\mathrm{s}})h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})\mathrm{d}\bm{r}_{\mathrm{s}},$
(22)

where $\hat{\bm{r}}_{\mathrm{s}}=\bm{r}_{\mathrm{s}}/\|\bm{r}_{\mathrm{s}}\|$.
We expand $h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})$ using the complex conjugate
of spherical harmonic functions as in Eq. (10). Thus, the binaural signals at
$\bm{r}_{\mathrm{o}}$, $y_{\mathrm{L,R}}(\bm{r}_{\mathrm{o}})$ are represented
using $\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})$ and $H_{\mathrm{L,R},n}^{m}$ as

$\displaystyle
y_{\mathrm{L,R}}(\bm{r}_{\mathrm{o}})=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\sqrt{4\pi}\mathrm{j}^{-n}\alpha_{n}^{m}(\bm{r}_{\mathrm{o}})H_{\mathrm{L,R},n}^{m}.$
(23)

The infinite sum is truncated up to a predefined truncation order in practice.
The distance of the loudspeakers in measuring
$h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})$, i.e., $R_{\mathrm{s}}$, is assumed to
be sufficiently large, and $h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})$ is
approximated as the transfer function from plane wave sources.

##  IV  Proposed method of binaural rendering from distributed microphone
array signals

###  IV.1  Estimation of expansion coefficients using distributed microphone
array

We estimate the expansion coefficients of spherical wavefunctions at the
arbitrary position $\bm{r}\in D$ using microphones distributed over $D$. To
achieve this, we apply the harmonic analysis of infinite order for the
soundfield. The microphone directivities are assumed to be known, and the
expansion coefficients of the directivity pattern of the $i$th microphone are
denoted by $c_{i,n}^{m}$. Then, the observation of the $i$th microphone at
$\bm{r}_{\mathrm{m},i}$, $s_{i}$, is described using the expansion
coefficients of the soundfield as

$\displaystyle s_{i}$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}c_{i,n}^{m\ast}\alpha_{n}^{m}(\bm{r}_{\mathrm{m},i})$
$\displaystyle=\bm{c}_{i}^{\mathsf{H}}\bm{\alpha}(\bm{r}_{\mathrm{m},i})$
$\displaystyle=\bm{c}_{i}^{\mathsf{H}}\bm{T}(\bm{r}_{\mathrm{m},i}-\bm{r})\bm{\alpha}(\bm{r}),$
(24)

where $\bm{\alpha}(\bm{r})\in\mathbb{C}^{\infty}$ and
$\bm{c}_{i}\in\mathbb{C}^{\infty}$ are the infinite-dimensional vectors
consisting of $\alpha_{n}^{m}(\bm{r})$ and $c_{i,n}^{m}$, respectively.
Equation (24) can be rewritten as

$\displaystyle\bm{s}=\bm{\Xi}(\bm{r})^{\mathsf{H}}\bm{\alpha}(\bm{r}),$ (25)

where $\bm{\Xi}(\bm{r})\in\mathbb{C}^{\infty\times I}$ is obtained as

$\displaystyle\bm{\Xi}(\bm{r})$
$\displaystyle=[(\bm{c}_{1}^{\mathsf{H}}\bm{T}(\bm{r}_{\mathrm{m},1}-\bm{r}))^{\mathsf{H}}\cdots(\bm{c}_{I}^{\mathsf{H}}\bm{T}(\bm{r}_{\mathrm{m},I}-\bm{r}))^{\mathsf{H}}]$
$\displaystyle=[\bm{T}(\bm{r}-\bm{r}_{\mathrm{m},1})\bm{c}_{1}\cdots\bm{T}(\bm{r}-\bm{r}_{\mathrm{m},I})\bm{c}_{I}].$
(26)

We used the property of the translation operator (9) to derive Eq. (26). The
expansion coefficients $\bm{\alpha}(\bm{r})$ is estimated as Ueno:IEEE_SPL2018

$\displaystyle\bm{\alpha}(\bm{r})=\bm{\Xi}(\bm{r})(\bm{\Psi}+\lambda\bm{I})^{-1}\bm{s},$
(27)

where $\lambda$ is a constant parameter and
$\bm{\Psi}=\bm{\Xi}(\bm{r})^{\mathsf{H}}\bm{\Xi}(\bm{r})\in\mathbb{C}^{I\times
I}$. From the property in Eq. (9), the $(i,i^{\prime})$th element of
$\bm{\Psi}$ becomes

$\displaystyle(\bm{\Psi})_{i,i^{\prime}}$
$\displaystyle=\bm{c}_{i}^{\mathsf{H}}\bm{T}(\bm{r}_{\mathrm{m},i}-\bm{r})\bm{T}(\bm{r}-\bm{r}_{\mathrm{m},i^{\prime}})\bm{c}_{i^{\prime}}$
$\displaystyle=\bm{c}_{i}^{\mathsf{H}}\bm{T}(\bm{r}_{\mathrm{m},i}-\bm{r}_{\mathrm{m},i^{\prime}})\bm{c}_{i^{\prime}}.$
(28)

Therefore, $\bm{\Psi}$ does not depend on the position $\bm{r}$ and depends
only on the microphone positions and directivities. Note that $\bm{c}_{i}$ is
typically modeled by low-order coefficients.

###  IV.2 Binaural rendering based on spherical wave decomposition

To consider the loudspeaker distance in the HRTF measurement, we represent the
binaural signals assuming that $h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})$ is the
transfer function from the point source. Thus, the soundfield $u(\bm{r})$ is
represented as the weighted integral of spherical waves from a point source on
$\partial\Omega$ as

$\displaystyle
u(\bm{r})=\int_{\partial\Omega}w_{\mathrm{sph}}(\bm{r}_{\mathrm{s}})G(\bm{r}-\bm{r}_{\mathrm{s}})\mathrm{d}\bm{r}_{\mathrm{s}},$
(29)

where $G(\bm{r}-\bm{r}_{\mathrm{s}})$ is the transfer function of the point
source, which is equivalent to the three dimensional (3D) free-field Green’s
function defined as acous

$\displaystyle
G(\bm{r}-\bm{r}_{\mathrm{s}})=\frac{e^{-\mathrm{j}k\|\bm{r}-\bm{r}_{\mathrm{s}}\|}}{4\pi\|\bm{r}-\bm{r}_{\mathrm{s}}\|}$
$\displaystyle\ \
=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}-\mathrm{j}kj_{n}(kr)h_{n}(kR_{\mathrm{s}})Y^{m}_{n}(\theta,\phi)Y_{n}^{m}(\theta_{\mathrm{s}},\phi_{\mathrm{s}})^{\ast}.$
(30)

We here use the same spherical surface $\partial\Omega$ as in the HRTF
measurements, and $\bm{r}_{\mathrm{s}}$ is the position vector on
$\partial\Omega$. The representation in Eq. (29) is referred to as the single-
layer potential singlelayer  or simple source formulationacous . The spherical
wave weight $w_{\mathrm{sph}}(\bm{r}_{\mathrm{s}})$ can be related to the
expansion coefficients $\alpha_{n}^{m}(\bm{r})$ as

$\displaystyle
w_{\mathrm{sph}}(\bm{r}_{\mathrm{s}})=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\frac{\sqrt{4\pi}\mathrm{j}}{kR_{\mathrm{s}}^{2}h_{n}(kR_{\mathrm{s}})}\alpha_{n}^{m}(\bm{r})Y_{n}^{m}(\theta_{\mathrm{s}},\phi_{\mathrm{s}}).$
(31)

This relation can be confirmed by substituting Eqs. (30) and (31) into Eq.
(29) and using the orthogonality of the spherical harmonic function.

Finally, the binaural signals at $\bm{r}$, $y_{\mathrm{L,R}}(\bm{r})$, are
obtained using $\alpha_{n}^{m}(\bm{r})$ and $H_{\mathrm{L,R},n}^{m}$ as

$\displaystyle y_{\mathrm{L,R}}(\bm{r})$
$\displaystyle=\int_{\partial\Omega}w_{\mathrm{sph}}(\bm{r}_{\mathrm{s}})h_{\mathrm{L,R}}(\bm{r}_{\mathrm{s}})\mathrm{d}\bm{r}_{\mathrm{s}}$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\frac{\sqrt{4\pi}\mathrm{j}}{kh_{n}(kR_{\mathrm{s}})}\alpha_{n}^{m}(\bm{r})H_{\mathrm{L,R},n}^{m}.$
(32)

The infinite sum is truncated up to a finite order in practice. It can be
confirmed that Eq. (32) corresponds to Eq. (23) in the plane-wave-
decomposition-based method when $R_{\mathrm{s}}$ becomes infinity because
$h_{n}(kR_{\mathrm{s}})$ can be approximated as $\mathrm{j}^{n+1}/k$ for far-
field directivity ($R_{\mathrm{s}}\to\infty$) acous . Similar formulations
have been derived to take into consideration the loudspeaker distance in the
context of near-field compensation of higher order Ambisonics
Daniel:AES116conv  and HRTF interpolation god ; Pollow:Acustica2012 . In Sect.
V, it is shown that this distance-compensated rendering method significantly
improves the amplitude response of the reproduced binaural signals.

###  IV.3 Adaptation to rotation and translation of listener’s head

The binaural signals should be reproduced with adaptation to the listener’s
head movement. In the plane-wave-decomposition-based methods, the binaural
signals are adapted to the listening position by the phase shift of plane wave
weights conv ; plane_trans . When the expansion center is moved from the
origin to $\bm{r}^{\prime}$, Eq. (18) is transformed to

$\displaystyle
u(\bm{r})=\frac{1}{4\pi}\int_{\mathbb{S}_{2}}w_{\mathrm{pw}}(\bm{\eta})e^{\mathrm{j}k\langle\bm{\eta},\bm{r}^{\prime}\rangle}e^{\mathrm{j}k\langle\bm{\eta},\bm{r}-\bm{r}^{\prime}\rangle}\mathrm{d}\bm{\eta}.$
(33)

Thus, the plane wave weight at the shifted position $\bm{r}^{\prime}$ is
obtained as

$\displaystyle
w_{\mathrm{pw}}(\bm{\eta};\bm{r}^{\prime})=w_{\mathrm{pw}}(\bm{\eta})e^{\mathrm{j}k\langle\bm{\eta},\bm{r}^{\prime}\rangle}.$
(34)

To rotate the binaural signals, Eq. (22) is computed using HRTFs for rotated
angles. When the HRTFs of the rotated angles are not available, some
interpolation techniques are necessary god ; god2 ; god3 .

In the methods based on the spherical wavefunction expansion, the translation
and rotation of the listening position are achieved by transforming the
estimated expansion coefficients. The center of the expansion coefficients is
changed from $\bm{r}$ to $\bm{r}^{\prime}$ using the translation operator in
Eq. (6). In the truncation-based methods as described in Sect. III.2, this
translation operator is calculated by truncating $\bm{T}(\cdot)$ up to a
predefined order. On the other hand, in the proposed method, the translation
is calculated as

$\displaystyle\bm{\alpha}(\bm{r}^{\prime})$
$\displaystyle=\bm{T}(\bm{r}^{\prime}-\bm{r})\bm{\alpha}(\bm{r})$
$\displaystyle=\bm{T}(\bm{r}^{\prime}-\bm{r})\bm{\Xi}(\bm{r})(\bm{\Psi}+\lambda\bm{I})^{-1}\bm{s}$
$\displaystyle=\bm{\Xi}(\bm{r}^{\prime})(\bm{\Psi}+\lambda\bm{I})^{-1}\bm{s}.$
(35)

Since $\bm{\Psi}$ does not depend on the position $\bm{r}$ or
$\bm{r}^{\prime}$, the expansion coefficients at the listening position can be
obtained by changing $\bm{\Xi}(\bm{r})$ to $\bm{\Xi}(\bm{r}^{\prime})$ with
fixed $(\bm{\Psi}+\lambda\bm{I})^{-1}\bm{s}$. Thus, the spherical wave weight
at the shifted position $\bm{r}^{\prime}$,
$w_{\mathrm{sph}}(\bm{r}_{\mathrm{s}};\bm{r}^{\prime})$, is obtained by
replacing $\bm{\alpha}(\bm{r})$ in Eq. (32) with
$\bm{\alpha}(\bm{r}^{\prime})$.

Adaptation to the head rotation is achieved using the rotated expansion
coefficients described with the Euler angles
$(\mathring{\alpha},\mathring{\beta},\mathring{\gamma})$ as

$\displaystyle\alpha_{n}^{m^{\prime}}(\bm{r};\mathring{\alpha},\mathring{\beta},\mathring{\gamma})=\sum_{m=-n}^{n}D_{m^{\prime},m}^{n}(\mathring{\alpha},\mathring{\beta},\mathring{\gamma})\alpha_{n}^{m}(\bm{r}),$
(36)

where $D_{m,m^{\prime}}^{n}(\cdot)$ is the element of the Wigner-$D$ matrix
Edmonds . By using the Wigner $D$-matrix
$\bm{R}(\mathring{\alpha},\mathring{\beta},\mathring{\gamma})\in\mathbb{C}^{\infty\times\infty}$
consisting of $D_{m,m^{\prime}}^{n}(\cdot)$, we can represent Eq. (36) in a
matrix form as

$\displaystyle\bm{\alpha}(\bm{r};\mathring{\alpha},\mathring{\beta},\mathring{\gamma})=\bm{R}(\mathring{\alpha},\mathring{\beta},\mathring{\gamma})\bm{\alpha}(\bm{r}).$
(37)

The expansion coefficient $\bm{\alpha}(\bm{r})$ is estimated up to a finite
order in practice. When the truncation order of $\bm{\alpha}(\bm{r})$ is $N$,
$\bm{R}(\mathring{\alpha},\mathring{\beta},\mathring{\gamma})$ becomes a
$(N+1)^{2}\times(N+1)^{2}$ block diagonal matrix Edmonds .

Our proposed method is summarized as follows. When the estimated expansion
coefficients are truncated up to $N$, the binaural signals at the listening
position $\bm{r}$ with the rotation angles
$(\mathring{\alpha},\mathring{\beta},\mathring{\gamma})$ are represented as

$\displaystyle
y_{\mathrm{L,R}}(\bm{r};\mathring{\alpha},\mathring{\beta},\mathring{\gamma})=\bm{\psi}_{\mathrm{L,R}}^{\mathsf{T}}\bm{\Lambda}\bar{\bm{R}}(\mathring{\alpha},\mathring{\beta},\mathring{\gamma})\bar{\bm{\Xi}}(\bm{r})(\bm{\Psi}+\lambda\bm{I})^{-1}\bm{s},$
(38)

where $\bar{\bm{\Xi}}\in\mathbb{C}^{(N+1)^{2}\times I}$ and
$\bar{\bm{R}}\in\mathbb{C}^{(N+1)^{2}\times(N+1)^{2}}$ are constructed by
truncating $\bm{\Xi}$ and $\bm{R}$, respectively, and
$\bm{\Lambda}\in\mathbb{C}^{(N+1)^{2}\times(N+1)^{2}}$ is a diagonal matrix
consisting of $\sqrt{4\pi}\mathrm{j}/kh_{n}(kR_{\mathrm{s}})$. Since the
conversion from microphone observations to binaural signals is a linear
operation, $y_{\mathrm{L,R}}$ is obtained by applying the
multiple–input–multiple–output (MIMO) finite impulse response (FIR) filter to
$\bm{s}$ in practical implementations.

##  V  Comparison of binaural rendering methods

We first compare binaural rendering accuracy between plane-wave-decomposition-
based and spherical-wave-decomposition-based methods described in Sects. III.3
and IV.2, which are hereafter denoted as PLN and SPH, respectively. The
expansion coefficient of the soundfield in the recording area
$\alpha_{n}^{m}(\bm{r})$ is assumed to be given in this section.

HRTFs were obtained by numerical computation using the boundary element method
(BEM). The software Mesh2HRTF mesh2hrtf ; mesh  was used. The head shape of a
subject was scanned by magnetic resonance imaging (MRI), and the mesh data of
$25968$ surface elements were used to generate HRTFs. On the basis of
reciprocity, HRTFs were obtained as transfer functions from the point source
at the entrance of the ear canal to the sampling points on the spherical
surface $\partial\Omega$. An acoustically rigid surface of the head was
assumed. The radius of $\partial\Omega$, $R_{\mathrm{s}}$, was
$1.5~{}\mathrm{m}$, which is the typically used distance for HRTF measurements
datahrtf ; datahrtf2 . The sampling points on $\partial\Omega$ were from
$0^{\circ}$ to $355^{\circ}$ for the azimuth angle and from $10^{\circ}$ to
$160^{\circ}$ for the zenith angle at intervals of $5^{\circ}$; therefore, the
total number of sampling points was $2232$. The sound velocity was
$346.2~{}\mathrm{m/s}$. The frequency range of interest was set as
$100$–$15000~{}\mathrm{Hz}$ at intervals of $100~{}\mathrm{Hz}$.

We assumed that a single point source was positioned at
$\bm{r}_{\mathrm{ps}}=[r_{\mathrm{ps}},\theta_{\mathrm{ps}},\phi_{\mathrm{ps}}]^{\mathsf{T}}$
in the free-field recording area, and the expansion coefficient
$\alpha_{n}^{m}(\bm{r})$ was analytically calculated. Thus, true binaural
signals were also obtained as the HRTF from the point source at
$\bm{r}_{\mathrm{ps}}$ using BEM. The truncation order for HRTFs in Eq. (10)
was set as $35$. The truncation order $N$ for the binaural rendering in Eqs.
(23) and (32) was set as $N=\min(\lceil ekR_{\mathrm{w}}/2\rceil,35)$
Ward:IEEE_J_SAP2001 , where $e$ is Napier’s constant and
$R_{\mathrm{w}}=0.45~{}\mathrm{m}$ is the average shoulder width of Japanese
males.

Binaural signals were evaluated for the left ear from the point sources on the
horizontal plane ($\theta_{\mathrm{ps}}=90^{\circ}$). We here added the
additional arguments $r_{\mathrm{ps}}$ and $\phi_{\mathrm{ps}}$ to
$y_{\mathrm{L}}(k)$ to indicate the location of the virtual point source of
each binaural signal. For evaluation, we introduce two measures. One is the
normalized mean square error (NMSE) defined as

$\displaystyle\mathrm{NMSE}(r_{\mathrm{ps}},\phi_{\mathrm{ps}},k_{f})=$
$\displaystyle\hskip
20.0pt10\log_{10}\frac{|{y}_{\mathrm{L}}^{\mathrm{est}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})-{y}_{\mathrm{L}}^{\mathrm{true}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})|^{2}}{|{y}_{\mathrm{L}}^{\mathrm{true}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})|^{2}},$
(39)

where ${y}_{\mathrm{L}}^{\mathrm{true}}$ and ${y}_{\mathrm{L}}^{\mathrm{est}}$
are the true and estimated binaural signals, respectively, and the subscript
$f$ denotes the index of frequency bins. The other is the spectral distortion
(SD) defined as SD

$\mathrm{SD}(r_{\mathrm{ps}},\phi_{\mathrm{ps}})=\sqrt{\frac{1}{F}\sum_{f}\left(20\log_{10}\frac{|{y}_{\mathrm{L}}^{\mathrm{est}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})|}{|{y}_{\mathrm{L}}^{\mathrm{true}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})|}\right)^{2}},$
(40)

where only the difference in amplitude response is evaluated. To compare the
methods based on plane wave and spherical wave decomposition, the amplitudes
of the binaural signals generated by PLN and SPH were normalized to remove the
amplitude bias as

$\displaystyle\hat{y}_{\mathrm{L}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})=\frac{y_{\mathrm{L}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})}{\sqrt{\frac{1}{F}\sum_{f}|y_{\mathrm{L}}(k_{f};r_{\mathrm{ps}},\phi_{\mathrm{ps}})|^{2}}},$
(41)

where $F$ is the number of frequency bins. The phase difference was also
corrected in the time domain using the cross correlation between the true and
reproduced signals.

Figure2a.eps.4(a) PLN Figure2b.eps.4(b) SPH (Proposed)

Figure 2: NMSEs with respect to azimuth angle of point source at
$r_{\mathrm{ps}}=2.0~{}\mathrm{m}$ and frequency. Figure 3: Amplitude response
of binaural signal of left ear for point source at
$(2.0~{}\mathrm{m},90^{\circ},-90^{\circ})$ with respect to frequency.

Figure 2 shows the NMSEs of the binaural signals generated by PLN and SPH with
respect to the azimuth angle of the point source at
$r_{\mathrm{ps}}=2.0~{}\mathrm{m}$ and frequency. A generally lower NMSE was
achieved by SPH than by PLN. In particular, the NMSE of PLN was relatively
large at high frequencies. As an example, the amplitude response of the
binaural signal for the point source at
$(2.0~{}\mathrm{m},90^{\circ},-90^{\circ})$ with respect to frequency is
plotted in Fig. 3. The dips in the amplitude response were accurately
reproduced by SPH. On the other hand, the frequencies having dips for PLN were
shifted from those for the true response at several frequencies, e.g., 7.7
kHz.

Figure 4: SD with respect to azimuth angle of point source at
$r_{\mathrm{ps}}=2.0~{}\mathrm{m}$. Figure 5: SD with respect to point source
distance.

The SD is plotted for $r_{\mathrm{ps}}=2.0~{}\mathrm{m}$ in Fig. 4 to evaluate
general accuracy with respect to the azimuth angle. Again, a low SD was
achieved by SPH. The SD of PLN was particularly high at around
$\phi_{\mathrm{ps}}=-90^{\circ}$, which is the point source direction of the
opposite side of the head from the left ear where the direct sound is
relatively weak. The SD averaged for the azimuth angle with respect to the
source distance $r_{\mathrm{ps}}$ is shown in Fig. 5. The average SD of SPH
was lower than that of PLN for all distances.

##  VI Design and evaluation of composite microphone array

Figure6a.eps2.2in(a) Single array Figure6b.eps2.2in(b) Composite array
Figure6c.eps2.2in(c) Spherical array

Figure 6: Geometry of microphone arrays for (a) single, (b) composite, and (c)
spherical arrays.

Although the array of microphones (nearly) uniformly arranged on a single
sphere can efficiently capture the soundfield, particularly for estimating the
expansion coefficients of the spherical wavefunctions at the array center, it
does not have flexibility in the microphone arrangement. One of the benefits
of our proposed soundfield estimation method based on infinite dimensional
analysis is its applicability to a wide range of microphone array geometries.
Therefore, we consider using multiple small microphone arrays. Such a
composite array system will have flexibility in the array placement and
scalability for adding/removing small microphone arrays, depending on the
setting of the region of interest $D$. We evaluate our soundfield estimation
method using a composite microphone array by numerical simulations.

For each small microphone array, we used one array of eight unidirectional
microphones as shown in Fig. 6(a), whose arrangement was identical to that of
commercially available 2nd-order ambisonics microphones (the shape of
tetragonal trapezohedron) OctoMic ; ambisonic_3Dprint . The composite
microphone array consisted of eight small microphone arrays as shown in Fig.
6(b), so the total number of microphones was 64. To cover the listening area,
which should be larger than the average size of the listener’s head, four
small arrays were equiangularly placed on the circle of $0.145~{}\mathrm{m}$
radius at two heights $z=\pm 0.025~{}\mathrm{m}$. For comparison, we also
employed a spherical array of omnidirectional microphones mounted on a rigid
baffle of $0.145~{}\mathrm{m}$ radius, where the arrangement of 64 microphones
was determined by the spherical $t$-design for $t=7$ t-design  (Fig. 6(c)).

The directivity pattern of the $i$th unidirectional microphone is represented
as

$\displaystyle
c_{i}(\theta_{\mathrm{m},i},\phi_{\mathrm{m},i})=\beta+(1-\beta)\langle\bm{\eta},\bm{\eta}_{\mathrm{m},i}\rangle,$
(42)

where $\beta\in[0,1]$ is a constant parameter, $\bm{\eta}_{\mathrm{m},i}$ and
$\bm{\eta}$ are the directions of the $i$th unidirectional microphone (peak of
directivity) and incident soundwave, respectively. Thus, the expansion
coefficients of the directivity pattern are derived as

$\displaystyle\bm{c}_{i}$
$\displaystyle=[c_{i,0}^{0},c_{i,1}^{-1},c_{i,1}^{0},c_{i,1}^{1}]^{\mathsf{T}}$
$\displaystyle=\left[\beta,\frac{\sqrt{4\pi}\mathrm{j}}{3}(1-\beta)Y_{1}^{-1}(\theta_{\mathrm{m},i},\phi_{\mathrm{m},i})^{*},\right.$
$\displaystyle\hskip
14.22636pt\frac{\sqrt{4\pi}\mathrm{j}}{3}(1-\beta)Y_{1}^{0}(\theta_{\mathrm{m},i},\phi_{\mathrm{m},i})^{\ast},$
$\displaystyle\hskip
14.22636pt\left.\frac{\sqrt{4\pi}\mathrm{j}}{3}(1-\beta)Y_{1}^{1}(\theta_{\mathrm{m},i},\phi_{\mathrm{m},i})^{\ast}\right]^{\mathsf{T}}.$
(43)

We set $\beta=1/2$ in the experiments described in this section.

We evaluated the binaural reproduction accuracy using the composite microphone
array (Fig. 6(b)) and spherical microphone array (Fig. 6(c)). Our proposed
soundfield estimation method based on infinite-dimensional analysis described
in Sect. IV.1 was used for the composite array. For the spherical array, the
truncation-based estimation method in Sect. III.2 was applied. The binaural
rendering method from the estimated expansion coefficients was the SPH for
both arrays. The HRTF dataset of NEUMANN KU-100 lebedev_hrtf  was used, which
was measured for $2702$ loudspeaker positions on a sphere of
$1.5~{}\mathrm{m}$ radius. The truncation order for the binaural rendering in
Eq. (10) was set as $7$ since the number of microphones was 64.

A point source was set at
$(1.5~{}\mathrm{m},0.0~{}\mathrm{m},0.0~{}\mathrm{m})$ in the recording area.
We generated binaural signals at the origin using the observed signals of the
microphone arrays, changing their center positions. The positions were
$21\times 21$ grid points inside a square region of $0.8~{}\mathrm{m}\times
0.8~{}\mathrm{m}$, and 60 directions were sampled at $\theta=\pi/2$ from the
directions of HRTFs. The NMSEs defined in Eq. $\eqref{eq:NMSE}$ were averaged
for directions and frequencies up to $1.6~{}\mathrm{kHz}$, and plotted for
each position in Figs. 8 and 8. Those along the $x$\- and $z$-axes are shown
in Fig. 9, compared with the case of the single small array (Fig. 6(a)). The
NMSE of the spherical array was particularly low around the origin, i.e., the
listening position. The composite array achieved a comparable or lower NMSE
than the spherical array in the off-center region on the $x$-$y$ plane at
$z=0$. The NMSEs at
$(-0.071~{}\mathrm{m},0.071~{}\mathrm{m},0.0~{}\mathrm{m})$ with respect to
the azimuth angle and frequency are plotted in Fig. 10.

Figure7a.eps.25(a) Composite array Figure7b.eps.25(b) Spherical array

Figure 7: Distribution of average NMSE on $x$-$y$-plane at $z=0$.

Figure8a.eps.25(a) Composite array Figure8b.eps.25(b) Spherical array

Figure 8: Distribution of average NMSE on $y$-$z$-plane at $x=0$.

Figure9a.eps.35(a) $x$-axis Figure9b.eps.35(b) $z$-axis

Figure 9: Average NMSEs along the $x$\- and $z$-axes.

Figure10a.eps.4(a) Composite array Figure10b.eps.4(b) Spherical array

Figure 10: NMSEs at
$(-0.071~{}\mathrm{m},0.071~{}\mathrm{m},0.0~{}\mathrm{m})$ with respect to
azimuth angle and frequency.

##  VII  Experiments in practical environment

###  VII.1  Evaluation of practical system

We developed a practical composite microphone array as shown in Fig. 11. Eight
small arrays consisting of eight unidirectional microphones were used, and
their arrangements were identical to those in Fig. 6(b). The directivity
pattern of several microphones was measured, and the parameter $\beta$ in Eq.
(43) was adapted to the measurements. We first evaluated the binaural
reproduction accuracy using this microphone array.

We reproduced the binaural signals at the origin, whose position was
$1.3~{}\mathrm{m}$ above the ground. We placed the composite array with its
center at positions A and B, and measured impulse responses of a loudspeaker
at $(1.5~{}\mathrm{m},0.0~{}\mathrm{m},0.0~{}\mathrm{m})$ using swept-sine
signals swept_sine .

  * •

Position A: $(0.0~{}\mathrm{m},0.0~{}\mathrm{m},0.0~{}\mathrm{m})$

  * •

Position B: $(-0.071~{}\mathrm{m},0.071~{}\mathrm{m},0.0~{}\mathrm{m})$

Again, we used the HRTF dataset of the NEUMANN KU-100 dummy head lebedev_hrtf
for binaural rendering. We compared the reproduced binaural signals with the
true signals measured using the NEUMANN KU-100 dummy head placed at the
origin. We also compared the binaural signals reproduced by a single small
microphone array (Fig. 6(a)) placed at positions A and B. The sampling
frequency was set at $48~{}\mathrm{kHz}$.

As evaluation measures, we used the interaural time difference (ITD) and
interaural level difference (ILD) defined as

$\displaystyle\mathrm{ITD}$ $\displaystyle=\mathop{\rm
arg~{}max}\limits_{\tau}\frac{\int_{0}^{T}\bar{y}_{\mathrm{L}}(t)\bar{y}_{\mathrm{R}}(t+\tau)\mathrm{d}t}{\sqrt{\int_{0}^{T}\bar{y}_{\mathrm{L}}(t)^{2}\bar{y}_{\mathrm{R}}(t+\tau)^{2}\mathrm{d}t}},$
(44) $\displaystyle\mathrm{ILD}$
$\displaystyle=10\log_{10}\frac{\int_{0}^{T}\bar{y}_{\mathrm{L}}(t)^{2}\mathrm{d}t}{\int_{0}^{T}\bar{y}_{\mathrm{R}}(t)^{2}\mathrm{d}t},$
(45)

where $\bar{y}_{\mathrm{L,R}}(t)$ is the binaural signal in the time domain
and $T$ is set as the signal length. These measures are important features of
sound localization in humans ITD_1600 . We removed reflections from the
signals as preprocessing and applied a low-pass filter of $1.6~{}\mathrm{kHz}$
cutoff frequency. This is because the accurate reproduction is limited to low
frequencies, as shown in the numerical simulation (Fig. 10), and ITD
contributes to the frequency range below $1.6~{}\mathrm{kHz}$ ITD_1600 .

Figures 12 and 13 show the ITD and ILD of the binaural signals reproduced by
the proposed composite array and single small array, compared with those of
the true binaural signals. The ITD and ILD are plotted with respect to the
azimuth angle at positions A and B. The ITD of both methods was close to the
true one, but that of the single array was slightly smaller than the true ITD
at $90^{\circ}$. The ILD at position A was accurately reproduced; however,
that of the single array at position B significantly deteriorated. In the
proposed composite array, the reproduced ILD at the both positions was close
to the true one.

Figure 11: Developed system using composite microphone array.

Figure12a.eps6.5cm(a) Position A Figure12b.eps6.5cm(b) Position B

Figure 12: True and reproduced ITDs at positions A and B.

Figure13a.eps6.5cm(a) Position A Figure13b.eps6.5cm(b) Position B

Figure 13: True and reproduced ILDs at positions A and B.

###  VII.2  Listening experiments

We evaluated our binaural reproduction method by listening experiments.
Perceptual experiments using multiple stimulus with hidden reference and
anchor (MUSHRA) mushra  were conducted to compare the quality of the binaural
signals generated using the composite array and single small array. Again, the
HRTF dataset of the NEUMANN KU-100 dummy head was used for the binaural
rendering. The sampling frequency was $48~{}\mathrm{kHz}$. The binaural
signals at the origin were reproduced by each microphone array at positions A
and B. The reference signals were the observation signals of the dummy head.
As an anchor, the signals obtained by filtering the reference signals using a
low-pass filter of $1.6~{}\mathrm{kHz}$ cutoff frequency were used. All the
test signals are summarized as follows.

  1. 1.

Reference (hidden reference): binaural signal observed with the dummy head.

  2. 2.

C1/composite array (A): binaural signals reproduced using the composite array
signals at position A

  3. 3.

C2/composite array (B): binaural signals reproduced using the composite array
signals at position B

  4. 4.

C3/single array (A): binaural signals reproduced using the single small array
signals at position A

  5. 5.

C4/single array (B): binaural signals reproduced using the single small array
signals at position B

  6. 6.

C5/anchor: signal obtained using a low-pass filter of $1.6~{}\mathrm{kHz}$
cutoff frequency with respect to the reference signal

Test participants were asked to rate the difference between the reference and
test signals on a scale from $0$ to $100$. The higher score means the test
signal is close to the reference signal. The test signals were repeatedly
played back and the participants listened via a headphone, and the
participants were able to seamlessly switch signals. Thirteen male subjects in
their 20s and 30s were included, and those who scored less than 90 on the
Reference more than twice were excluded from the evaluation. We used two
signals, music (acoustic guitar sound) and male speech, for sound sources,
taken from a dataset Araki:LVAICA2021 . The binaural signals of three source
directions, $\phi=45^{\circ}$, $90^{\circ}$, and $135^{\circ}$, obtained by
rotational operation were included in one sequence.

Figure 14 shows the box plot of perceptual experiment scores. The median
scores indicated that the proposed composite array (C1 and C2) and single
array (C3 and C4) differ with 95% confidence for both sound sources music and
male speech. It was validated that the evaluation data followed the normal
distribution by using the Lilliefors test ($p$-values $>$ 0.01). Then, the
$p$-values of the Tukey–Kramer multiple comparison test were obtained as shown
in Tables 2 (music) and 2 (male speech). The significant difference in mean
score between C1 and C3 (i.e., position A) and that between C2 and C4
(position B) were found in both sound sources. These results indicate that the
perceptual quality of the binaural signals reproduced with the composite array
was higher than that with the small array. Furthermore, there was no
significant difference between C1 and C2; therefore, our proposed method is
robust against the changes in listening position.

Figure14a.eps7cm(a) Music Figure14b.eps7cm(b) Male speech

Figure 14: Box plot of listening test scores. The center bar of each box
indicates the median score. The boxes and whiskers show the interquartile
range (IQR) and 1.5 times IQR, respectively. The notches in the box refer to
the 95% confidence interval of the median. Table 1: $p$-values (Tukey–Kramer)
for music.

C5 | 1 |  |  |  |
---|---|---|---|---|---
C4 | 0.0027 | 1 |  |  |
C3 | 0 | 0.4183 | 1 |  |
C2 | 0 | 0 | 0.0041 | 1 |
C1 | 0 | 0 | 0.0016 | 0.9997 | 1
Method | C5 | C4 | C3 | C2 | C1

Table 2: $p$-values (Tukey–Kramer) for male speech.

C5 | 1 |  |  |  |
---|---|---|---|---|---
C4 | 0.0003 | 1 |  |  |
C3 | 0.0002 | 1 | 1 |  |
C2 | 0 | 0 | 0.0001 | 1 |
C1 | 0 | 0 | 0 | 0.6550 | 1
Method | C5 | C4 | C3 | C2 | C1

##  VIII Conclusion

We proposed a method of binaural rendering from distributed microphone array
signals. For capturing a soundfield inside a spatial region, the soundfield
recording method based on harmonic analysis of infinite order is applied,
which allows the use of arbitrarily placed microphones in the recording area.
Typically used spherical microphone arrays are sometimes impractical because
the array size should be inherently larger than the head size, and the
microphone arrangement should be uniform on a sphere as much as possible.
Since our proposed method has flexibility and scalability in the microphone
arrangement, we investigated a combination of small microphone arrays, which
will be useful for capturing a soundfield in a large region. We also proposed
a spherical-wave-decomposition-based binaural rendering method taking into
consideration the distance in HRTF measurement. Compared with the conventional
plane-wave-decomposition-based rendering method, the peaks and dips in the
amplitude response of the binaural signals are accurately reproduced. We
developed a practical composite microphone array system and evaluated the
system by measuring the reproduced ITD and ILD, and conducting listening
experiments. The experimental results indicated that our proposed method is
robust against changes in listening position inside a range of the region in
the recording area.

##  IX Acknowledgments

We are grateful to Professor Shuichi Sakamoto for providing us with head mesh
data. This work was supported by JST, PRESTO Grant Number JPMJPR18J4.

## References

  * (1) A. J. Berkhout, D. de Vries, and P. Vogel, “Acoustic control by wave field synthesis,” J. Acoust. Soc. Amer. 93(5), 2764–2778 (1993).
  * (2) J. Daniel, S. Moureau, and R. Nicol, “Further investigations of high-order ambisonics and wavefield synthesis for holophonic sound imaging,” in _Proc. 114th AES Conv._ , Amsterdam (2003).
  * (3) M. A. Poletti, “Three-dimensional surround sound systems based on spherical harmonics,” J. Audio Eng. Soc. 53(11), 1004–1025 (2005).
  * (4) S. Spors, R. Rabenstein, and J. Ahrens, “The theory of wave field synthesis revisited,” in _Proc. 124th AES Conv._ , Amsterdam (2008).
  * (5) S. Koyama, K. Furuya, Y. Hiwasaki, and Y. Haneda, “Analytical approach to wave field reconstruction filtering in spatio-temporal frequency domain,” IEEE Trans. Audio, Speech, Lang. Process. 21(4), 685–696 (2013).
  * (6) J. Blauert, _Spatial Hearing, Revised edition: The Psychophysics of Human Sound Localization_ (The MIT Press, Cambridge, 1996).
  * (7) T. R. Thurlow and P. S. Runge, “Effect of induced head movement on localization of direction of sounds,” J. Acoust. Soc. Amer. 42, 480–488 (1967).
  * (8) I. Toshima, H. Uematsu, and T. Hirahara, “A steerable dummy head that tracks three-dimensional head movement: TeleHead,” Acoust. Sci. Tech. 24(5), 327–329 (2003).
  * (9) S. Sakamoto, S. Hongo, T. Okamoto, Y. Iwaya, and Y. Suzuki, “Sound-space recording and binaural presentation system based on a 252-channel microphone array,” Acoust. Sci. Tech. 36(6), 516–526 (2015).
  * (10) E. Rasumow, M. Blau, S. Doclo, M. Hansen, S. V. de Par, D. Püschel, and V. Mellert, “Least squares versus non-linear cost functions for a virtual artificial head,” J. Acoust. Soc. Amer. 133(5) (2013).
  * (11) E. Rasumow, M. Hansen, S. van de Par, D. Püschel, V. Mellert, S. Doclo, and M. Blau, “Regularization approaches for synthesizing HRTF directivity patterns,” IEEE/ACM Trans. Audio, Speech, Lang. Process. 24(2), 215–225 (2016).
  * (12) W. Song, W. Ellermeier, and J. Hald, “Using beamforming and binaural synthesis for the psychoacoustical evaluation of target sources in noise,” J. Acoust. Soc. Amer. 123(2), 910–924 (2007).
  * (13) N. R. Shabtai and B. Rafaely, “Generalized spherical array beamforming for binaural speech reproduction,” IEEE/ACM Trans. Audio, Speech, Lang. Process. 22(1), 238–247 (2014).
  * (14) J. Atkins, “Robust beamforming and steering of arbitrary beam patterns using spherical arrays,” in _Proc. IEEE Int. Workshop Appl. Signal Process. Audio Acoust. (WASPAA)_ , New Paltz (2011), pp. 237–240.
  * (15) R. Duraiswami, D. N. Zotkin, Z. Li, E. Grassi, N. Gumerov, and L. S. Davis, “High order spatial capture and its binaural head-tracked playback over headphones with HRTF cues,” in _Proc. 119th AES Conv_ , New York (2005).
  * (16) M. Park and B. Rafaely, “Sound-field analysis by plane-wave decomposition using spherical microphone array,” J. Acoust. Soc. Amer. 118(5), 3094–3103 (2005).
  * (17) D. N. Zotkin, R. Duraiswami, and N. A. Gumerov, “Sound field decomposition using spherical microphone arrays,” in _Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)_ , Las Vegas (2008), pp. 277–280.
  * (18) P. Samarasinghe, T. D. Abhayapala, and M. A. Poletti, “Wavefield analysis over large areas using distributed higher order microphones,” IEEE/ACM Trans. Audio, Speech, Lang. Process. 22(3), 647–658 (2014).
  * (19) N. Ueno, S. Koyama, and H. Saruwatari, “Sound field recording using distributed microphones based on harmonic analysis of infinite order,” IEEE Signal Process. Lett. 25(1), 135–139 (2018).
  * (20) M. Nakanishi, N. Ueno, S. Koyama, and H. Saruwatari, “Two-dimensional sound field recording with multiple circular microphone arrays considering multiple scattering,” in _Proc. IEEE Int. Workshop Appl. Signal Process. Audio Acoust. (WASPAA)_ , New Paltz (2019), pp. 363–367.
  * (21) D. Menzies and M. Al-Akaidi, “Nearfield binaural synthesis and ambisonics,” J. Acoust. Soc. Amer. 121(3), 1559–1563 (2007).
  * (22) L. Birnie, T. D. Abhayapala, P. Samarasinghe, and V. Tourbabin, “Sound field translation methods for binaural reproduction,” in _Proc. IEEE Int. Workshop Appl. Signal Process. Audio Acoust. (WASPAA)_ , New Paltz (2019), pp. 135–139.
  * (23) F. Schultz and S. Spors, “Data-based binaural synthesis including rotational and translatory head-movements,” in _Proc. 52nd AES Int. Conf._ (2013).
  * (24) B. Bernschütz, A. Giner, C. Pörschmann, and J. Arend, “Binaural reproduction of plane waves with reduced modal order,” Acta Acustica United Acustica 100 (2014).
  * (25) M. Zaunschirm, C. Schörkhuber, and R. Höldrich, “Binaural rendering of ambisonic signals by head-related impulse response time alignment and a diffuseness contraint,” J. Acoust. Soc. Amer. 143(6), 3616–3627 (2018).
  * (26) N. Iijima, S. Koyama, and H. Saruwatari, “Binaural rendering from distributed microphone signals considering loudspeaker distance in measurements,” in _Proc. IEEE Int. Workshop Multimedia Signal Process. (MMSP)_ , Tampere (2020).
  * (27) P. A. Martin, _Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles_ (Cambridge University Press, New York, 2006).
  * (28) K. Watanabe, Y. Iwaya, Y. Suzuki, S. Takane, and S. Sato, “Dataset of head-related transfer functions measured with a circular loudspeaker array,” Acoust. Sci. Tech. 35(3), 159–165 (2014).
  * (29) D. N. Zotkin, R. Duraiswami, and N. A. Gumerov, “Interpolation and range expolation of HRTFs,” in _Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)_ , Montreal (2004), pp. 45–48.
  * (30) K.-V. Nguyen, T. Carpentier, M. Noisternig, and O. Warusfel, “Calculation of head related transfer functions in the proximity recording using spherical harmonic decomposition: Comparison with measurements and evaluations,” in _Proc. Int. Symp. Ambisonics, Spherical Acoust._ , Paris (2010).
  * (31) W. Zhang, T. D. Abhayapala, R. A. Kennedy, and R. Duraiswami, “Insights into head-related transfer function: Spatial dimensionality and continuous representation,” J. Acoust. Soc. Amer. 127(4), 2347–2357 (2010).
  * (32) B. Rafaely, “Analysis and design of spherical microphone arrays,” IEEE Trans. Audio, Speech, Lang. Process. 13(1), 135–143 (2005).
  * (33) I. Balmages and B. Rafaely, “Open-sphere designs for spherical microphone arrays,” IEEE Trans. Audio, Speech, Lang. Process. 15(2), 727–732 (2007).
  * (34) S. Koyama, K. Furuya, K. Wakayama, S. Shimauchi, and H. Saruwatari, “Analytical approach to transforming filter design for sound field recording and reproduction using circular arrays with a spherical baffle,” J. Acoust. Soc. Amer. 139(3), 1024–1036 (2016).
  * (35) D. Colton and R. Kress, _Inverse Acoustic and Electromagnetic Scattering Theory_ , Vol. 93 (Springer, New York, 2013).
  * (36) E. G. Williams, _Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography_ (Academic Press, London, 1999).
  * (37) J. Daniel and S. Moreau, “Further study of sound field coding with higher order ambisonics,” in _Proc. 116th AES Conv._ , Berlin (2004).
  * (38) M. Pollow, K.-V. Nguyen, O. Warusfel, T. Capentier, M. Müller-Trapet, M. Vorländer, and M. Noisternig, “Calculation of head-related transfer functions for arbitrary field points using spherical harmonic decomposition,” Acta Acust. united with Acust. 98, 72–82 (2012).
  * (39) A. R. Edmonds, _Angular Momentum in Quantum Mechanics_ (Princeton University Press, Princeton, 1957).
  * (40) H. Ziegelwanger, P. Majdak, and W. Kreuzer, “Numerical calculation of listener-specific head-related transfer functions and sound localization: Microphone model and mesh discretization,” J. Acoust. Soc. Amer. 138(1), 208–222 (2015).
  * (41) H. Ziegelwanger, W. Kreuzer, and P. Majdak, “Mesh2HRTF: Open-source software package for the numerical calculation of head-related transfer functions,” in _Proc. Int. Congr. Sound Vibr. (ICSV)_ , Florence (2015).
  * (42) C. Armstrong, L. Thresh, D. Murphy, and G. Kearney, “A perceptual evaluation of individual and non-individual HRTFs: A case study of the SADIE II database,” Appl. Sci. 8(11) (2018).
  * (43) D. B. Ward and T. D. Abhayapala, “Reproduction of a plane-wave sound field using an array of loudspeakers,” IEEE Trans. Speech Audio Process. 9(6), 697–707 (2001).
  * (44) F. C. Tommasini, O. A. Ramos, M. X. Hug, and F. R. Bermejo, “Usage of spectral distortion for objective evaluation of personalized HRTF in the median plane,” Int. J. Acoust. Vibr. (IJAV) 20(2), 81–89 (2015).
  * (45) “Core sound octomic™ 2nd-order ambisonic microphone” , https://www.core-sound.com/products/octomic.
  * (46) F. Lopez-Lezcano, “The sphear project update: Refining the octasphear, a 2nd order ambisonics microphone,” in _Proc. EAA Spatial Audio Signal Process. Symp._ , Paris (2019), pp. 103–108.
  * (47) X. Chen and R. S. Womersley, “Existence of solutions to systems of underdetermined equations and spherical designs,” SIAM J. Numer. Anal. 44(6), 2326–2341 (2006).
  * (48) J. M. Arend, A. Neidhardt, and C. Pörschmann, “Measurement and perceptual evaluation of a spherical near-field HRTF set,” in _Proc. 29th Tonmeistertagung - VDT Int. Conv._ , Cologne (2016), pp. 356–363.
  * (49) Y. Suzuki, F. Asano, H.-Y. Kim, and T. Sone, “An optimum computer-generated pulse signal suitable for the measurement of very long impulse responses,” J. Acoust. Soc. Amer. 97, 1119–1123 (1995).
  * (50) W. A. Yost and R. R. Fay, eds., _Human Psychophysics_ (Springer-Verlag New York, 1993).
  * (51) International Telecommunication Union, Radiocommunication Sector, “Recommendation ITU-R BS – 1534-3 method for the subjective assessment of intermediate quality level of audio systems,” (2014).
  * (52) S. Araki, F. Nesta, E. Vincent, Z. Koldovsky, G. Nolte, A. Ziehe, and A. Benichoux, “The 2011 signal separation evaluation campaign (SiSEC2011) – audio source separation –,” in _Proc. Int. Conf. Latent Variable Analysis Signal Separation (LVA/ICA_ (2012), pp. 414–422.