#  Weighted pressure matching based on kernel interpolation
for sound field reproduction

Shoichi KOYAMA and Kazuyuki ARIKAWA(1)

### ABSTRACT

A sound field reproduction method called weighted pressure matching is
proposed. Sound field reproduction is aimed at synthesizing the desired sound
field using multiple loudspeakers inside a target region. Optimization-based
methods are derived from the minimization of errors between synthesized and
desired sound fields, which enable the use of an arbitrary array geometry in
contrast with integral-equation-based methods. Pressure matching is widely
used in the optimization-based sound field reproduction methods because of its
simplicity of implementation. Its cost function is defined as the synthesis
errors at multiple control points inside the target region; then, the driving
signals of the loudspeakers are obtained by solving a least-squares problem.
However, in pressure matching, the region between the control points is not
taken into consideration. We define the cost function as the regional
integration of the synthesis error over the target region. On the basis of the
kernel interpolation of the sound field, this cost function is represented as
the weighted square error of the synthesized pressures at the control points.
Experimental results indicate that the proposed weighted pressure matching
outperforms conventional pressure matching.

Keywords: Sound field reproduction, Kernel interpolation, Weighted pressure
matching

## 1 INTRODUCTION

Sound field reproduction is aimed at synthesizing spatial sound using multiple
loudspeakers (or secondary sources). Such a technique can be applied to
virtual/augmented reality audio, generation of multiple sound zones for
personal audio, and noise cancellation in a spatial region.

Sound field reproduction methods can be classified into two major categories:
integral-equation-based and optimization-based methods. The integral-equation-
based methods are developed from the boundary integral representations derived
from the Helmholtz equation, such as wave field synthesis and higher-order
ambisonics [3, 19, 17, 1, 24, 12]. The optimization-based methods are derived
from the minimization of a certain cost function defined for synthesized and
desired sound fields in a target region, such as pressure matching and mode
matching [16, 9, 7, 17, 4, 21, 11]. Many integral-equation-based methods
require the array geometry of loudspeakers to have a simple shape, such as a
sphere, plane, circle, or line, and driving signals are obtained from a
discrete approximation of an integral equation. In optimization-based methods,
the loudspeakers can be placed arbitrarily, and driving signals are generally
derived as a closed-form least-squares solution. In particular, pressure
matching is widely used among the optimization-based methods because of its
simplicity of implementation. Pressure matching is based on synthesizing the
desired pressures at a discrete set of control points placed over the target
region.

An issue of pressure matching is that the region between the control points is
not taken into consideration because of the discrete approximation. Therefore,
its reproduction accuracy can deteriorate when the distribution of the control
points is not sufficiently dense. However, the smaller the number of control
points, the better in practice, because the transfer functions between the
loudspeakers and control points are generally measured in advance. We propose
an optimization-based sound field reproduction method called weighted pressure
matching. We define the cost function as the regional integration of the
synthesis error over the target region. On the basis of the kernel
interpolation of the sound field [20, 22], this cost function is represented
by the pressures at the control points with the regional integration of kernel
functions. When the same kernel functions are used for interpolating primary
and secondary sound fields, i.e., the desired field and the sound field of
each loudspeaker, respectively, the driving signal is obtained as the solution
of weighted least squares problem with the weighting matrix defined by using
the kernel functions. Experimental evaluation comparing pressure matching and
weighted pressure matching is performed.

## 2 PROBLEM STATEMENT AND PRIOR WORKS

### 2.1 Problem formulation

Figure 1: The desired sound field is synthesized in the target region $\Omega$
using multiple secondary sources.

Suppose that $L$ secondary sources (loudspeakers) are placed around a target
region $\Omega\subset\mathbb{R}^{3}$ as shown in Fig. 1. The sound field
$u_{\mathrm{syn}}(\bm{r},\omega)$ at the position $\bm{r}\in\mathbb{R}^{3}$
and angular frequency $\omega\in\mathbb{R}$ synthesized using the secondary
sources is represented as

$\displaystyle
u_{\mathrm{syn}}(\bm{r},\omega)=\sum_{l=1}^{L}d_{l}(\omega)g_{l}(\bm{r},\omega),$
(1)

where $d_{l}(\omega)$ is the driving signal of the $l$th secondary source
($l\in\\{1,\ldots,L\\}$), and $g_{l}(\bm{r},\omega)$ is the transfer function
from the $l$th secondary source to the position $\bm{r}$. The transfer
functions $g_{l}(\bm{r},\omega)$ are assumed to be known by measuring or
modeling them in advance. Hereafter, the angular frequency $\omega$ is omitted
for notational simplicity.

The goal of sound field reproduction is to obtain $\\{d_{l}\\}_{l=1}^{L}$ of
the $L$ secondary sources so that $u_{\mathrm{syn}}(\bm{r})$ coincides with
the desired sound field, denoted by $u_{\mathrm{des}}(\bm{r})$, inside
$\Omega$. We define the cost function to determine the driving signal
$\\{d_{l}\\}_{l=1}^{L}$ as

$\displaystyle J$
$\displaystyle=\int_{\Omega}\left|\sum_{l=1}^{L}d_{l}g_{l}(\bm{r})-u_{\mathrm{des}}(\bm{r})\right|^{2}\mathrm{d}\bm{r}.$
(2)

The optimal driving signal can be obtained by solving the minimization problem
of $J$.

### 2.2 Pressure matching

Since it is difficult to solve the minimization problem of $J$ owing to the
regional integration over $\Omega$, several methods based on the approximation
of $J$ have been proposed. A simple strategy for solving it is to discretize
the target region $\Omega$ into multiple control points, which is referred to
as the pressure matching. Assume that $N$ control points are placed over
$\Omega$ and their positions are denoted by $\bm{r}_{\mathrm{c},n}$
($n\in\\{1,\ldots,N\\}$). The cost function $J$ is approximated as the error
between the synthesized and desired pressures at the control points. The
optimization problem of pressure matching is written as

$\displaystyle\operatorname*{minimize}_{\bm{d}\in\mathbb{C}^{L}}\|\bm{Gd}-\bm{u}^{\mathrm{des}}\|^{2}+\eta\|\bm{d}\|^{2},$
(3)

where $\bm{d}=[d_{1},\ldots,d_{L}]^{\mathsf{T}}\in\mathbb{C}^{L}$ is the
vector of the driving signals,
$\bm{u}^{\mathrm{des}}=[u_{\mathrm{des}}(\bm{r}_{\mathrm{c},1}),\ldots,u_{\mathrm{des}}(\bm{r}_{\mathrm{c},N})]^{\mathsf{T}}\in\mathbb{C}^{N}$
is the vector of the desired sound pressures, and
$\bm{G}\in\mathbb{C}^{N\times L}$ is the matrix consisting of the transfer
functions $g_{l}(\bm{r}_{\mathrm{c},n})$ between $L$ secondary sources and $N$
control points. The second term is the regularization term to prevent an
excessively large amplitude of $\bm{d}$, and $\eta$ is a constant parameter.
The closed-form solution of Eq. (3) is obtained as

$\displaystyle\hat{\bm{d}}=\left(\bm{G}^{\mathsf{H}}\bm{G}+\eta\bm{I}\right)^{-1}\bm{G}^{\mathsf{H}}\bm{u}^{\mathrm{des}}.$
(4)

Another strategy to approximately solve Eq. (2) is to represent the sound
field by spherical wavefunction expansion [23, 14], which is referred to as
mode matching [17]. The driving signal is obtained so that the expansion
coefficients of the synthesized sound field coincide with those of the desired
sound field. The mode matching method is generalized by introducing a
weighting matrix consisting of the regional integration of the spherical
wavefunctions, which is called weighted mode matching [21]. Instead of
empirical truncation of the expansion order in mode matching, optimal weights
on the expansion coefficients are applied in weighted mode matching.

## 3 WEIGHTED PRESSURE MATCHING

Since pressure matching is based on the discrete approximation of the target
region, the region between the control points is taken into consideration. We
consider incorporating a sound field interpolation method into pressure
matching. First, we introduce the kernel interpolation method for sound
fields. Then, the weighted pressure matching is formulated by representing the
sound field using the kernel interpolation.

### 3.1 Kernel interpolation of sound field

The goal of the sound field interpolation is to estimate the pressure
distribution $u(\bm{r})$ from the discrete set of pressure measurements
$s_{m}$ at $\bm{r}_{m}$ ($m\in\\{1,\ldots,M\\}$). This interpolation problem
is formulated as

$\displaystyle\operatorname*{minimize}_{u\in\mathcal{H}}\sum_{m=1}^{M}|u(\bm{r}_{m})-s_{m}|^{2}+\lambda\|u\|_{\mathcal{H}}^{2},$
(5)

where $\mathcal{H}$ is a function space for which we seek a solution,
$\|\cdot\|_{\mathcal{H}}$ is a norm on $\mathcal{H}$, and $\lambda$ is a
constant parameter. When $\mathcal{H}$ is a function space called a
reproducing kernel Hilbert space, the optimization problem (5) corresponds to
kernel ridge regression, for which a closed-form solution can be obtained
[15]. Here, $\mathcal{H}$ is assumed to be a reproducing kernel Hilbert space
with the inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}}$ and the
positive-definite reproducing kernel
$\kappa:\mathcal{H}\times\mathcal{H}\to\mathbb{C}$. On the basis of the
representer theorem [18], the solution for Eq. (5) is obtained as

$\displaystyle
u(\bm{r})=\bm{\kappa}(\bm{r})^{\mathsf{T}}\left(\bm{K}+\lambda\bm{I}\right)^{-1}\bm{s},$
(6)

where

$\displaystyle\bm{\kappa}(\bm{r})$
$\displaystyle=\begin{bmatrix}\kappa(\bm{r},\bm{r}_{1})&\ldots&\kappa(\bm{r},\bm{r}_{M})\end{bmatrix}^{\mathsf{T}}$
(7) $\displaystyle\bm{K}$
$\displaystyle=\begin{bmatrix}\kappa(\bm{r}_{1},\bm{r}_{1})&\cdots&\kappa(\bm{r}_{1},\bm{r}_{M})\\\
\vdots&\ddots&\vdots\\\
\kappa(\bm{r}_{M},\bm{r}_{1})&\cdots&\kappa(\bm{r}_{M},\bm{r}_{M})\end{bmatrix}$
(8) $\displaystyle\bm{s}$
$\displaystyle=\begin{bmatrix}s_{1}&\ldots&s_{M}\end{bmatrix}^{\mathsf{T}}.$
(9)

Next, it is necessary to define appropriate $\mathcal{H}$ and
$\langle\cdot,\cdot\rangle_{\mathcal{H}}$, as well as $\kappa$. The pressure
field $u$ inside the source-free simply connected region $\Omega$ can be
modeled as a solution of the homogeneous Helmholtz equation as

$\displaystyle(\nabla^{2}+k^{2})u=0,$ (10)

where $\nabla^{2}$ is the Laplacian and $k=\omega/c$ is the wavenumber defined
with the sound velocity $c$. Any solution of Eq. (10) can be well approximated
by the superposition of plane waves, i.e., the Herglotz wavefunction [6, 22],
as

$\displaystyle\hat{u}(\bm{r})=\frac{1}{4\pi}\int_{\mathbb{S}_{2}}\tilde{u}(\bm{\xi})\mathrm{e}^{-\mathrm{j}\bm{k}^{\mathsf{T}}\bm{r}}\mathrm{d}\bm{\xi},$
(11)

where $\mathbb{S}_{2}$ is the unit sphere, $\tilde{u}(\bm{\xi})$ is the
(square-integrable) complex amplitude of the plane wave of arrival direction
$\bm{\xi}\in\mathbb{S}_{2}$, and $\bm{k}=-k\bm{\xi}$ is the wave vector. Using
this representation, we define the inner product and norm over the Hilbert
space $\mathcal{H}$ as

$\displaystyle\langle u_{1},u_{2}\rangle_{\mathcal{H}}$
$\displaystyle=\frac{1}{4\pi}\int_{\mathbb{S}_{2}}\frac{1}{\gamma(\bm{\xi})}\tilde{u}_{1}(\bm{\xi})^{\ast}\tilde{u}_{2}(\bm{\xi})\mathrm{d}\bm{\xi},$
(12) $\displaystyle\|u\|_{\mathcal{H}}$ $\displaystyle=\sqrt{\langle
u,u\rangle_{\mathcal{H}}},$ (13)

respectively. Here, $\gamma(\bm{\xi})$ is a directional weighting function,
which is introduced to incorporate prior knowledge on source directions. We
set the kernel function $\kappa(\bm{r}_{1},\bm{r}_{2})$ as

$\displaystyle\kappa(\bm{r}_{1},\bm{r}_{2})=\frac{1}{4\pi}\int_{\mathbb{S}_{2}}\gamma(\bm{\xi})\mathrm{e}^{-\mathrm{j}\bm{k}^{\mathsf{T}}(\bm{r}_{1}-\bm{r_{2}})}\mathrm{d}\bm{\xi}.$
(14)

It can be confirmed that $\kappa(\bm{r}_{1},\bm{r}_{2})$ is the reproducing
kernel of $\mathcal{H}$ because the inner product of
$\kappa(\bm{r},\bm{r}^{\prime})$ and $u(\bm{r})$ equals $u(\bm{r}^{\prime})$.

We define the directional weighting function as

$\displaystyle\gamma(\bm{\xi})=\mathrm{e}^{\rho\bm{\xi}^{\mathsf{T}}\hat{\bm{r}}},$
(15)

where $\rho\geq 0$ is a constant parameter and $\hat{\bm{r}}\in\mathbb{S}_{2}$
represents the prior arrival direction of the source. This function is derived
from the von Mises–Fisher distribution in directional statistics [13]. For
$\rho>0$, one can find that the smaller the
$\rho\bm{\xi}^{\mathsf{T}}\hat{\bm{r}}$ is, the larger the norm
$\|u\|_{\mathcal{H}}$ becomes, and vice versa. Therefore, the regularization
term in Eq. (5) becomes larger when the difference between the prior arrival
direction $\hat{\bm{r}}$ and the direction of $\bm{\xi}$ becomes larger. When
$\rho=0$, the weighting function becomes uniform ($\gamma(\bm{\xi})=1$). The
shape of $\gamma(\bm{\xi})$ becomes sharper with increasing $\rho$.

By substituting Eq. (15) into Eq. (14), we can derive the kernel function with
directional weighting as

$\displaystyle\kappa(\bm{r}_{1},\bm{r}_{2})$
$\displaystyle=\frac{1}{4\pi}\int_{\mathbb{S}_{2}}\mathrm{e}^{\rho\bm{\xi}^{\mathsf{T}}\hat{\bm{r}}}\cdot\mathrm{e}^{-\mathrm{j}\bm{k}^{\mathsf{T}}(\bm{r}_{1}-\bm{r_{2}})}\mathrm{d}\bm{\xi}$
$\displaystyle=j_{0}\left(\left[\left(\mathrm{j}\rho\sin\theta\cos\phi-
kx_{12}\right)^{2}+\left(\mathrm{j}\rho\sin\theta\sin\phi-
ky_{12}\right)^{2}+\left(\mathrm{j}\rho\cos\theta-
kz_{12}\right)^{2}\right]^{\frac{1}{2}}\right),$ (16)

where $j_{0}(\cdot)$ is the 0th-order spherical Bessel function of the first
kind, $\phi$ and $\theta$ are the azimuth and zenith angles of $\hat{\bm{r}}$,
and $\bm{r}_{1}-\bm{r}_{2}=[x_{12},y_{12},z_{12}]^{\mathsf{T}}$. By setting
$\rho=0$, we can simplify Eq. (16), and obtain the kernel function of uniform
weighting as

$\displaystyle\kappa(\bm{r}_{1},\bm{r}_{2})$
$\displaystyle=\frac{1}{4\pi}\int_{\mathbb{S}_{2}}\mathrm{e}^{-\mathrm{j}\bm{k}^{\mathsf{T}}(\bm{r}_{1}-\bm{r_{2}})}\mathrm{d}\bm{\xi}$
$\displaystyle=j_{0}(k\|\bm{r}_{1}-\bm{r}_{2}\|).$ (17)

Thus, the kernel interpolation of a sound field is achieved by using Eq. (6)
with the kernel function in Eq. (16) or (17). These kernel functions are
derived for a three-dimensional (3D) sound field, but similar kernel functions
can be derived for a two-dimensional (2D) sound field [10].

### 3.2 Weighted pressure matching based on kernel interpolation

We apply the kernel interpolation in the sound field reproduction to
accurately approximate the cost function $J$ in Eq. (2) from the pressures at
the control points. A similar idea has been applied in the context of spatial
active noise control [8, 10, 2] and multizone sound field control [5]. The
transfer functions $g_{l}(\bm{r})$ and desired sound field
$u_{\mathrm{des}}(\bm{r})$ are interpolated from those at the control points
as

$\displaystyle\hat{g}_{l}(\bm{r})$
$\displaystyle=\bm{\kappa}_{l}(\bm{r})^{\mathsf{T}}\left(\bm{K}_{l}+\lambda\bm{I}\right)^{-1}\bm{g}_{l}:=\bm{z}_{l}(\bm{r})^{\mathsf{T}}\bm{g}_{l}$
(18) $\displaystyle\hat{u}_{\mathrm{des}}$
$\displaystyle=\bm{\kappa}^{\mathrm{des}}(\bm{r})^{\mathsf{T}}\left(\bm{K}^{\mathrm{des}}+\lambda\bm{I}\right)^{-1}\bm{u}^{\mathrm{des}}:=\bm{z}^{\mathrm{des}}(\bm{r})^{\mathsf{T}}\bm{u}^{\mathrm{des}},$
(19)

where $\bm{g}_{l}\in\mathbb{C}^{L}$ is the $l$th column vector of $\bm{G}$,
and
$\bm{\kappa}_{l}(\bm{r}),\bm{\kappa}^{\mathrm{des}}(\bm{r})\in\mathbb{C}^{N}$
and $\bm{K}_{l},\bm{K}^{\mathrm{des}}\in\mathbb{C}^{N\times N}$ are
respectively the vectors and matrices consisting of the kernel function
defined with the positions $\\{\bm{r}_{\mathrm{c},n}\\}_{n=1}^{N}$. Then, the
cost function $J$ can be approximated using
$\hat{\bm{g}}(\bm{r})=[\bm{z}_{1}(\bm{r})^{\mathsf{T}}\bm{g}_{1},\ldots,\bm{z}_{L}(\bm{r})^{\mathsf{T}}\bm{g}_{L}]^{\mathsf{T}}$
as

$\displaystyle J$
$\displaystyle\approx\int_{\Omega}\left|\sum_{l=1}^{L}d_{l}\hat{g}_{l}(\bm{r})-\hat{u}_{\mathrm{des}}(\bm{r})\right|^{2}\mathrm{d}\bm{r}$
$\displaystyle=\int_{\Omega}\left|\hat{\bm{g}}(\bm{r})^{\mathsf{T}}\bm{d}-\bm{z}^{\mathrm{des}}(\bm{r})^{\mathsf{T}}\bm{u}^{\mathrm{des}}\right|^{2}\mathrm{d}\bm{r}$
$\displaystyle=\bm{d}^{\mathsf{H}}\bm{W}_{gg}\bm{d}-\bm{d}^{\mathsf{H}}\bm{W}_{gu}\bm{u}^{\mathrm{des}}+C,$
(20)

where

$\displaystyle\bm{W}_{gg}$
$\displaystyle=\int_{\Omega}\hat{\bm{g}}(\bm{r})^{\ast}\bm{g}(\bm{r})^{\mathsf{T}}\mathrm{d}\bm{r}$
(21) $\displaystyle\bm{W}_{gu}$
$\displaystyle=\int_{\Omega}\hat{\bm{g}}(\bm{r})^{\ast}\bm{z}^{\mathrm{des}}(\bm{r}^{\mathsf{T}})\mathrm{d}\bm{r},$
(22)

and $C$ is the term not including $\bm{d}$. Therefore, the optimal driving
signal $\bm{d}$ is obtained by solving

$\displaystyle\operatorname*{minimize}_{\bm{d}\in\mathbb{C}^{L}}\bm{d}^{\mathsf{H}}\bm{W}_{gg}\bm{d}-\bm{d}^{\mathsf{H}}\bm{W}_{gu}\bm{u}^{\mathrm{des}}+\eta\|\bm{d}\|^{2}.$
(23)

Again, the regularization term is added. This minimization problem also has
the closed-form solution as

$\displaystyle\hat{\bm{d}}=\left(\bm{W}_{gg}+\eta\bm{I}\right)^{-1}\bm{W}_{gu}\bm{u}^{\mathrm{des}}.$
(24)

Since the directions of the secondary sources are generally known, the kernel
functions for interpolating $g_{l}(\bm{r})$ as well as $\bm{z}_{l}(\bm{r})$
can be defined in advance by using Eq. (16). When the desired sound field is
set with the parameters of source directions and positions, $\hat{\bm{r}}$ for
$u_{\mathrm{des}}(\bm{r})$ can also be given. When it is difficult to set
$\hat{\bm{r}}$ in advance, the kernel function of uniform weighting Eq. (17)
will be appropriate for defining $\bm{z}_{l}(\bm{r})$ and/or
$\bm{z}^{\mathrm{des}}(\bm{r})$.

When using the same kernel function for interpolating $g_{l}(\bm{r})$ and
$u_{\mathrm{des}}(\bm{r})$, e.g., the kernel function of uniform weighting Eq.
(17), the cost function $J$ in (20) can be further simplified. By defining
$\bm{z}(\bm{r})^{\mathsf{T}}=\bm{\kappa}(\bm{r})^{\mathsf{T}}(\bm{K}+\lambda\bm{I})^{-1}$,
we can obtain $J$ as

$\displaystyle J$
$\displaystyle\approx\int_{\Omega}\left|\bm{z}(\bm{r})^{\mathsf{T}}\left(\bm{Gd}-\bm{u}^{\mathrm{des}}\right)\right|^{2}$
$\displaystyle=\left(\bm{Gd}-\bm{u}^{\mathrm{des}}\right)^{\mathsf{H}}\bm{W}\left(\bm{Gd}-\bm{u}^{\mathrm{des}}\right),$
(25)

where

$\displaystyle\bm{W}=\int_{\Omega}\bm{z}(\bm{r})^{\ast}\bm{z}(\bm{r})^{\mathsf{T}}\mathrm{d}\bm{r}.$
(26)

Thus, the optimal driving signal is derived as

$\displaystyle\hat{\bm{d}}=\left(\bm{G}^{\mathsf{H}}\bm{W}\bm{G}+\eta\bm{I}\right)^{-1}\bm{G}^{\mathsf{H}}\bm{W}\bm{u}^{\mathrm{des}}.$
(27)

This driving signal can be regarded as the solution of the weighted mean
square error between synthesized and desired pressures at the control points.
Therefore, we refer to the proposed method as weighted pressure matching. The
weighted pressure matching enables the increase in the reproduction accuracy
of pressure matching only by introducing the weighting matrix. Note that the
matrices $\bm{W}$, $\bm{W}_{gg}$, and $\bm{W}_{gu}$ can be computed only with
the positions of the control points and the target region $\Omega$ by defining
the kernel function using Eq. (16) and/or (17).

## 4 EXPERIMENTS

Figure 2: Experimental setup. The target region was set as a 2D square region.
Red dots and blue crosses indicate loudspeakers and control points,
respectively.

We conducted numerical experiments to evaluate the proposed method. Pressure
matching and weighted pressure matching are hereafter denoted as PM and WPM,
respectively. We also evaluated WPM with directional weighting, which is
hereafter denoted as WPM (directional). In the experiments, we assumed a 2D
free field.

As shown in Fig. 2, 12 loudspeakers were placed at the same intervals along
the border of a square with dimensions of $2.0\leavevmode\nobreak\
\mathrm{m}\times 2.0\leavevmode\nobreak\ \mathrm{m}$. The target region
$\Omega$ was set as a square region of $1.0\leavevmode\nobreak\
\mathrm{m}\times 1.0\leavevmode\nobreak\ \mathrm{m}$. The centers of these
squares were at the origin. Sixteen control points were placed at the same
intervals over the target region. In Fig. 2, the loudspeakers and control
points are indicated by red dots and blue crosses, respectively. Each
loudspeaker was assumed to be a point source. The desired sound field was set
to be a single plane-wave field, whose propagation direction was
$\pi/4\leavevmode\nobreak\ \mathrm{rad}$.

$\bm{u}^{\mathrm{des}}$ and $\bm{G}$ were given as pressure measurements at
the control points. In WPM, the kernel function was set to be uniform, i.e.,
Eq. (17). The kernel function in WPM (directional) was the directional kernel
given by Eq. (16), where the direction $\hat{\bm{r}}$ was set as the true
directions of primary and secondary sources, and the parameter $\rho$ was
$5.0$. The regularization parameter $\lambda$ in Eq. (6) and $\eta$ in Eqs.
(4), (24), and (27) were set as $10^{-6}$. For evaluation measure, we define
the signal-to-distortion ratio (SDR) as

$\displaystyle\mathrm{SDR}(\omega)=\frac{\int_{\Omega}|u_{\mathrm{des}}(\bm{r},\omega)|^{2}\mathrm{d}\bm{r}}{\int_{\Omega}|u_{\mathrm{syn}}(\bm{r},\omega)-u_{\mathrm{des}}(\bm{r},\omega)|^{2}\mathrm{d}\bm{r}},$
(28)

where the integration was computed at the evaluation points regularly
distributed over the target region. The intervals of the evaluation points
were $10^{-2}\leavevmode\nobreak\ \mathrm{m}$.

Figure 3: SDR with respect to frequency.

(a) PM (b) WPM (c) WPM (Directional)

Figure 4: Synthesized pressure distributions at $450\leavevmode\nobreak\
\mathrm{Hz}$.

(a) PM (b) WPM (c) WPM (Directional)

Figure 5: Square error distributions at $450\leavevmode\nobreak\ \mathrm{Hz}$.
The SDRs of PM, WPM, and WPM (directional) were $11.9$, $17.3$, and
$18.3\leavevmode\nobreak\ \mathrm{dB}$, respectively.

The SDR with respect to the frequency is plotted in Fig. 5. The SDRs were
larger than $20\leavevmode\nobreak\ \mathrm{dB}$ at frequencies below
$390\leavevmode\nobreak\ \mathrm{Hz}$. Above $400\leavevmode\nobreak\
\mathrm{Hz}$, the SDRs of WPM and WPM (directional) were larger than that of
PM. In particular, WPM (directional) maintained a large SDR at high
frequencies. In Figs. 5 and 5, the synthesized pressure distribution and
square error distribution at $450\leavevmode\nobreak\ \mathrm{Hz}$ of each
method are shown. The region of small errors in WPM was larger than that of
PM. The SDRs of PM, WPM, and WPM (directional) were $11.9$, $17.3$, and
$18.3\leavevmode\nobreak\ \mathrm{dB}$, respectively. It can be considered
that these differences are the effect of the interpolation. By taking into
consideration the region between the control points, we can improve the
reproduction accuracy. Since the interpolation accuracy of the directional
kernel is high because of the use of prior information on the source
directions, WPM (directional) achieved the highest reproduction accuracy.

## 5 CONCLUSION

We proposed a sound field reproduction method called weighted pressure
matching. Pressure matching is a widely used optimization-based sound field
reproduction method because of its simplicity. Since pressure matching is
based on the synthesis of desired pressures at a discrete set of control
points distributed over the target region, the region between the control
points is not taken into consideration. Our cost function is defined as the
regional integration of the synthesis error over the target region. On the
basis of the kernel interpolation of sound fields, the driving signal of
weighted pressure matching is obtained as the weighted least squares solution
with the weighting matrix consisting of the regional integration of the kernel
functions. In the numerical experiments, the weighted pressure matching
achieved high reproduction accuracy compared with conventional pressure
matching, especially at high frequencies. By introducing the directional
weighting for interpolating the sound field, we can further increase the
reproduction accuracy.

## ACKNOWLEDGMENTS

This work was supported by JST FOREST Program (Grant Number JPMJFR216M,
Japan).

## REFERENCES

##

  * [1] J. Ahrens and S. Spors.  An analytical approach to sound field reproduction using circular and spherical loudspeaker distributions.  Acta Acust. united Acust., 94:988–999, 2008.
  * [2] K. Arikawa, S. Koyama, and H. Saruwatari.  Spatial active noise control based on individual kernel interpolation of primary and secondary sound fields.  In Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), pages 1056–1060, Singapore, May 2022.
  * [3] A. J. Berkhout, D. de Vries, and P. Vogel.  Acoustic control by wave field synthesis.  J. Acoust. Soc. Amer., 93(5):2764–2778, 1993.
  * [4] T. Betlehem and T. D. Abhayapala.  Theory and design of sound field reproduction in reverberant environment.  J. Acoust. Soc. Amer., 117(4):2100–2111, 2005.
  * [5] J. Brunnström, S. Koyama, and M. Moonen.  Variable span trade-off filter for sound zone control with kernel interpolation weighting.  In Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), pages 1071–1075, Singapore, May 2022.
  * [6] D. Colton and R. Kress.  Inverse Acoustic and Electromagnetic Scattering Theory.  Springer, New York, NY, USA, 2013.
  * [7] 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, Netherlands, 2003.
  * [8] H. Ito, S. Koyama, N. Ueno, and H. Saruwatari.  Feedforward spatial active noise control based on kernel interpolation of sound field.  In Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), pages 511–515, Brighton, May 2019.
  * [9] O. Kirkeby, P. A. Nelson, F. O. Bustamante, and H. Hamada.  Local sound field reproduction using digital signal processing.  J. Acoust. Soc. Amer., 100(3):1584–1593, 1996.
  * [10] S. Koyama, J. Brunnström, H. Ito, N. Ueno, and H. Saruwatari.  Spatial active noise control based on kernel interpolation of sound field.  IEEE/ACM Trans. Audio, Speech, Lang. Process., 29:3052–3063, 2021\.
  * [11] S. Koyama, G. Chardon, and L. Daudet.  Optimizing source and sensor placement for sound field control: An overview.  IEEE/ACM Trans. Audio, Speech, Lang. Process., 28:686–714, 2020\.
  * [12] 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\.
  * [13] K. V. Mardia and P. E. Jupp.  Directional Statistics.  John Wiley & Sons, Chichester, UK, 2009.
  * [14] P. A. Martin.  Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles.  Cambridge University Press, New York, 2006.
  * [15] K. P. Murphy.  Machine Learning: A Probabilistic Perspective.  MIT Press, Cambridge, UK, 2012.
  * [16] P. A. Nelson.  Active control of acoustic fields and the reproduction of sound.  J. Sound Vibr., 177(4):447–477, 1993.
  * [17] M. A. Poletti.  Three-dimensional surround sound systems based on spherical harmonics.  J. Audio Eng. Soc., 53(11):1004–1025, 2005.
  * [18] B. Schölkopf, R. Herbrich, and A. J. Smola.  A generalized representer theorem.  In Proc. Int. Conf. Comput. Learn. Theory (COLT), pages 416–426, Amsterdam, Netherlands, Jul. 2001.
  * [19] S. Spors, R. Rabenstein, and J. Ahrens.  The theory of wave field synthesis revisited.  In Proc. 124th AES Conv., Amsterdam, Netherlands, 2008.
  * [20] 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.
  * [21] N. Ueno, S. Koyama, and H. Saruwatari.  Three-dimensional sound field reproduction based on weighted mode-matching method.  IEEE/ACM Trans. Audio, Speech, Lang. Process., 27(12):1852–1867, 2019.
  * [22] N. Ueno, S. Koyama, and H. Saruwatari.  Directionally weighted wave field estimation exploiting prior information on source direction.  IEEE Trans. Signal Process., 69:2383–2395, 2021.
  * [23] E. G. Williams.  Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography.  Academic Press, London, UK, 1999.
  * [24] Y. J. Wu and T. D. Abhayapala.  Theory and design of soundfield reproduction using continuous loudspeaker concept.  IEEE Trans. Audio, Speech, Lang. Process., 17(1):107–116, 2009\.