# LOCALIZING ACOUSTIC ENERGY IN SOUND FIELD SYNTHESIS
BY DIRECTIONALLY WEIGHTED EXTERIOR RADIATION SUPPRESSION

###### Abstract

A method for synthesizing the desired sound field while suppressing the
exterior radiation power with directional weighting is proposed. The exterior
radiation from the loudspeakers in sound field synthesis systems can be
problematic in practical situations. Although several methods to suppress the
exterior radiation have been proposed, suppression in all outward directions
is generally difficult, especially when the number of loudspeakers is not
sufficiently large. We propose the directionally weighted exterior radiation
representation to prioritize the suppression directions by incorporating it
into the optimization problem of sound field synthesis. By using the proposed
representation, the exterior radiation in the prioritized directions can be
significantly reduced while maintaining high interior synthesis accuracy,
owing to the relaxed constraint on the exterior radiation. Its performance is
evaluated with the application of the proposed representation to amplitude
matching in numerical experiments.

Index Terms—  Sound field synthesis, multizone sound field control, personal
audio, exterior radiation, directional weighting

##  1 Introduction

Sound field synthesis/reproduction/control aims at synthesizing the desired
spatial sound inside a target region using multiple loudspeakers (or secondary
sources). One or more finite interior regions of the target are surrounded by
secondary sources in typical sound field synthesis systems.

A large number of sound field synthesis techniques have been proposed for
various applications. In spatial audio for virtual/augmented reality, the
desired pressure (i.e., magnitude and phase) distribution is synthesized in
the target region. In contrast to methods based on analytical representations
derived from boundary integral equations, such as wave field synthesis and
higher-order ambisonics [1, 2, 3, 4, 5, 6, 7], methods based on the
minimization of the squared error between synthesized and desired sound
fields, such as pressure matching and mode matching [8, 4, 9, 10], have
practical advantages since the array geometry of the secondary sources can be
arbitrary. These methods can also be extended to spatial active noise control
by combining with sound field estimation methods [11, 12]. In multizone sound
field control for personal audio systems [13, 14, 15, 16], the desired sound
field is sometimes set using the magnitude distribution or acoustic potential
energy since the desired propagation direction of the synthesized sound can be
unspecified.

Most sound field synthesis methods described above are based on solving the
minimization problem of synthesis error inside the target region; however, the
exterior region of the target is usually not taken into consideration. Thus,
the exterior radiation power of the secondary sources can be significantly
large. Several techniques to synthesize the desired sound field in the target
region while suppressing the exterior radiation power, thus localizing the
acoustic energy, have been proposed [17, 9, 18, 19, 20]. Basically, these
methods are based on imposing a constraint on the exterior radiation power by
a penalty term or (in)equality constraint in addition to the cost function of
the interior synthesis error. However, suppressing the exterior radiation in
all the outward directions is very difficult particularly when the number of
secondary sources is not sufficiently large. An excessive constraint on the
exterior radiation can also affect the synthesis accuracy in the interior
region.

In some practical situations, allowable directions of exterior radiation can
be given in advance, e.g., directions where there are no people. To exploit
this prior information on the priority directions of the suppression and relax
the constraint on the exterior radiation power, we propose the directionally
weighted exterior radiation suppression in sound field synthesis. Since the
proposed directionally weighted exterior radiation power is represented by a
quadratic form of the driving signal of the secondary sources, it is simply
applicable to various sound field synthesis methods. We also demonstrate its
application to amplitude matching, i.e., a method to synthesize the desired
amplitude (or magnitude) distribution over the target region, and its
performance evaluation by numerical experiments.

##  2 Sound Field Synthesis With Exterior Radiation Suppression

Suppose that the desired sound field is synthesized in a target region
$\Omega\subset\mathbb{R}^{3}$ by $L$ secondary sources at
$\\{\bm{r}_{l}\\}_{l=1}^{L}$. All the secondary sources are assumed to exist
inside a spherical region $\Omega_{\mathrm{S}}$ ($\supset\Omega$), with the
center $\bm{r}_{\mathrm{S}}$ (coordinate origin) and radius $R_{\mathrm{S}}$,
and outside $\Omega$, i.e., $\bm{r}_{l}\in\Omega_{\mathrm{S}}\backslash\Omega$
(see Fig. 1). The synthesized sound field $u_{\mathrm{syn}}$ at the position
$\bm{r}$ and angular frequency $\omega$ is represented by a linear combination
of the transfer functions of secondary sources
$\\{g_{l}(\bm{r},\omega)\\}_{l=1}^{L}$ as

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

where $d_{l}\in\mathbb{C}$ is the driving signal of $l$th secondary source. We
also define a vector of the driving signals
$\bm{d}=[d_{1},\ldots,d_{L}]^{\mathsf{T}}\in\mathbb{C}^{L}$. Hereafter,
$\omega$ is omitted for notational simplicity.

Synthesizing the desired sound field inside $\Omega$ while suppressing the
exterior radiation power is basically achieved by solving the following
optimization problem:

$\displaystyle\mathop{\rm
minimize}\limits_{\bm{d}\in\mathbb{C}^{L}}\mathcal{J}(\bm{d})+\gamma\mathcal{E}(\bm{d})+\alpha\mathcal{R}(\bm{d}),$
(2)

where $\mathcal{J}(\bm{d})$ evaluates the synthesis error inside $\Omega$,
$\mathcal{E}(\bm{d})$ is a measure for the exterior radiation power of the
secondary sources, $\mathcal{R}(\bm{d})$ is the regularization term for
$\bm{d}$ (typically, $\mathcal{R}(\bm{d}):=\|\bm{d}\|^{2}$), and $\gamma$ and
$\alpha$ are positive constants. In (weighted) pressure matching,
$\mathcal{J}(\bm{d})$ is formulated as the squared synthesis error of pressure
at control points distributed over $\Omega$ [8, 10]. $\mathcal{J}(\bm{d})$ can
also be represented by using spherical wavefunction expansion coefficients of
$\\{g_{l}(\bm{r})\\}_{l=1}^{L}$ and the desired sound field as in (weighted)
mode matching [4, 9]. In amplitude matching [16], squared synthesis error of
magnitude at the control points is used for $\mathcal{J}(\bm{d})$ for
multizone sound field control. Although (2) is defined by adding
$\mathcal{E}(\bm{d})$ as a penalty term, it is also possible to impose an
(in)equality constraint on $\mathcal{E}(\bm{d})$ [19, 20]. Another approach is
to maximize the ratio between interior and exterior acoustic potential
energies as in acoustic contrast control [15]. In any case, it is necessary to
appropriately define the exterior radiation power $\mathcal{E}(\bm{d})$ as a
simple function of $\bm{d}$ as possible.

Fig. 1: Sound field synthesis using multiple secondary sources.

##  3 Directionally Weighted Exterior Radiation Suppression

In [18], the exterior radiation power $\mathcal{E}(\bm{d})$ is formulated as
an angular integral of the acoustic intensity in the outward normal direction.
However, it is sometimes difficult to suppress exterior radiation power in all
directions without sufficient secondary sources. To enable relaxing the
constraint on the exterior radiation power by using prior knowledge on
allowable directions of exterior radiation, we formulate directionally
weighted exterior radiation power.

###  3.1 Formulation of Directionally Weighted Exterior Radiation Power

The exterior radiation power from $\partial\Omega_{\mathrm{S}}$, i.e., the
boundary of $\Omega_{\mathrm{S}}$, is defined as

$\displaystyle\mathcal{E}_{\mathrm{uni}}(\bm{d})=\int_{\partial\Omega_{\mathrm{S}}}I_{r}(\bm{r})\mathrm{d}\bm{r},$
(3)

where $I_{r}(\bm{r})$ is the acoustic intensity in the outward normal
direction of $\partial\Omega_{\mathrm{S}}$, represented by $r$:

$\displaystyle
I_{r}(\bm{r})=\frac{1}{2}\mathrm{Re}\left[u_{\mathrm{syn}}(\bm{r})\frac{\mathrm{j}}{\rho
ck}\frac{\partial}{\partial r}u_{\mathrm{syn}}(\bm{r})^{\ast}\right].$ (4)

Here, $c$ is the sound speed, $k:=\omega/c$ is the wave number, $\rho$ is the
density of air, $\mathrm{j}$ is the imaginary unit, and $\mathrm{Re}[\cdot]$
represents the real part of a complex number. We refer to
$\mathcal{E}_{\mathrm{uni}}(\bm{d})$ defined in (3) as uniformly weighted
exterior radiation power. In [18], it is shown that
$\mathcal{E}_{\mathrm{uni}}(\bm{d})$ is represented by a quadratic form of
$\bm{d}$ with a matrix consisting of spherical wavefunction expansion
coefficients of $\\{g_{l}(\bm{r})\\}_{l=1}^{L}$ and translation operator.

We define directionally weighted exterior radiation power
$\mathcal{E}_{\mathrm{dir}}(\bm{d})$ as

$\displaystyle\mathcal{E}_{\mathrm{dir}}(\bm{d})=\int_{\partial\Omega_{\mathrm{S}}}I_{r}(\bm{r})w(\theta,\phi)\mathrm{d}\bm{r},$
(5)

where $w:\mathbb{S}_{2}\to\mathbb{R}_{\geq 0}$ is the directional weighting
function, and $\theta$ and $\phi$ are respectively zenith and azimuth angles.
By using $w$, one can design a directional priority of exterior radiation
penalized in the cost function (2).

We show that $\mathcal{E}_{\mathrm{dir}}(\bm{d})$ is also represented as a
quadratic form of $\bm{d}$ in a similar manner to
$\mathcal{E}_{\mathrm{uni}}(\bm{d})$. First, $w$ is expanded by spherical
harmonic function $Y_{n}^{m}$ with order $n$ and degree $m$ around
$\bm{r}_{\mathrm{S}}$ as

$\displaystyle
w(\theta,\phi)=\sum_{n,m}\tilde{w}_{n}^{m}(\bm{r}_{\mathrm{S}})Y_{n}^{m}(\theta,\phi),$
(6)

where $\tilde{w}_{n}^{m}(\bm{r}_{\mathrm{S}})$ is the expansion coefficients
and $\sum_{n,m}:=\sum_{n=1}^{\infty}\sum_{m=-n}^{n}$. Besides,
$u_{\mathrm{syn}}$ is expanded by spherical wavefunction for exterior field
around $\bm{r}_{\mathrm{S}}$ for
$\bm{r}=(r,\theta,\phi)\in\mathbb{R}^{3}\backslash\Omega_{\mathrm{S}}$ as

$\displaystyle
u_{\mathrm{syn}}(\bm{r})=\sum_{n,m}\mathring{u}_{\mathrm{syn},n}^{m}(\bm{r}_{\mathrm{S}})h_{n}(kr)Y_{n}^{m}(\theta,\phi),$
(7)

where $\mathring{u}_{\mathrm{syn},n}^{m}(\bm{r}_{\mathrm{S}})$ is the
expansion coefficients, and $h_{n}$ is the $n$th-order spherical Hankel
function of the first kind. By substituting (6) and (7) into (5),
$\mathcal{E}_{\mathrm{dir}}(\bm{d})$ is written as

$\displaystyle\mathcal{E}_{\mathrm{dir}}(\bm{d})=$
$\displaystyle\quad\frac{1}{2\rho
ck}\mathrm{Re}\Big{[}\mathrm{j}\sum_{\nu,\mu}\tilde{w}_{\nu}^{\mu}(\bm{r}_{\mathrm{S}})\sum_{n,m}\sum_{n^{\prime},m^{\prime}}\mathring{u}_{\mathrm{syn},n}^{m}(\bm{r}_{\mathrm{S}})\mathring{u}_{\mathrm{syn},n^{\prime}}^{m^{\prime}}(\bm{r}_{\mathrm{S}})^{\ast}$
$\displaystyle\quad\cdot h_{n}(kR_{\mathrm{S}})\frac{\partial}{\partial
R_{\mathrm{S}}}h_{n^{\prime}}(kR_{\mathrm{S}})^{\ast}$
$\displaystyle\quad\cdot
R_{\mathrm{S}}^{2}\int_{0}^{2\pi}\mathrm{d}\phi\int_{0}^{\pi}\mathrm{d}\theta
Y_{n}^{m}(\theta,\phi)Y_{n^{\prime}}^{m^{\prime}}(\theta,\phi)^{\ast}Y_{\nu}^{\mu}(\theta,\phi)\sin\theta\Big{]}.$
(8)

The integral of the triple product of spherical harmonic functions in (8) is
rewritten as

$\displaystyle\int_{0}^{2\pi}\mathrm{d}\phi\int_{0}^{\pi}\mathrm{d}\theta
Y_{n}^{m}(\theta,\phi)Y_{n^{\prime}}^{m^{\prime}}(\theta,\phi)^{\ast}Y_{\nu}^{\mu}(\theta,\phi)\sin\theta$
$\displaystyle=\sum_{q=|n-\nu|}^{n+\nu}\mathcal{G}(n,m;\nu,\mu;q)\delta_{q,n^{\prime}}\delta_{m+\mu,m^{\prime}},$
(9)

where $\mathcal{G}$ is referred to as Gaunt coefficients [21]. Thus,
$\mathcal{E}_{\mathrm{dir}}(\bm{d})$ is written as

$\displaystyle\mathcal{E}_{\mathrm{dir}}(\bm{d})$
$\displaystyle=-\frac{R_{\mathrm{S}}^{2}}{2\rho
c}\mathrm{Im}\Big{[}\sum_{\nu,\mu}\sum_{n,m}\sum_{n^{\prime},m^{\prime}}\tilde{w}_{\nu,\mu}(\bm{r}_{\mathrm{S}})\mathring{u}_{\mathrm{syn},n}^{m}(\bm{r}_{\mathrm{S}})\mathring{u}_{\mathrm{syn},n^{\prime}}^{m^{\prime}}(\bm{r}_{\mathrm{S}})^{\ast}$
$\displaystyle\quad\cdot
h_{n}(kR_{\mathrm{S}})h_{n^{\prime}}^{\prime}(kR_{\mathrm{S}})^{\ast}\sum_{q=|n-\nu|}^{n+\nu}\mathcal{G}(n,m;\nu,\mu;q)\delta_{q,n^{\prime}}\delta_{m+\mu,m^{\prime}}\Big{]}$
$\displaystyle=-\frac{R_{\mathrm{S}}^{2}}{2\rho
c}\mathrm{Im}\Big{[}\sum_{\nu,\mu}\sum_{n,m}\sum_{n^{\prime}=|n-\nu|}^{n+\nu}\tilde{w}_{\nu,\mu}(\bm{r}_{\mathrm{S}})\mathring{u}_{\mathrm{syn},n}^{m}(\bm{r}_{\mathrm{S}})\mathring{u}_{\mathrm{syn},n^{\prime}}^{m^{\prime}}(\bm{r}_{\mathrm{S}})^{\ast}$
$\displaystyle\quad\hskip 20.0pt\cdot
h_{n}(kR_{\mathrm{S}})h_{n^{\prime}}^{\prime}(kR_{\mathrm{S}})^{\ast}\mathcal{G}(n,m;\nu,\mu;n^{\prime})\Big{]}.$
(10)

Algorithm 1 ADMM algorithm for amplitude matching with exterior radiation
suppression

0:
$\bm{d}^{(0)},\bm{\lambda}^{(0)},\bm{a}_{\mathrm{des}},\bm{G},\bm{A},\gamma,\alpha,\xi$

1: $t=0$

2: while stopping criterion not satisfied do

3: $\bm{h}^{(t)}=\bm{Gd}^{(t)}+\bm{\lambda}^{(t)}/\xi$

4: $\bm{\theta}^{(t+1)}=\arg(\bm{h}^{(t)})$

5: $\bm{a}^{(t+1)}=(\xi|\bm{h}^{(t)}|+2\bm{a}_{\mathrm{des}})/(\xi+2)$

6:
$\bm{d}^{(t+1)}=\left(\bm{G}^{\mathsf{H}}\bm{G}+(2\gamma/\xi)\bm{A}+(2\alpha/\xi)\bm{I}\right)^{-1}$$\quad\quad\quad\quad\cdot\bm{G}^{\mathsf{H}}\left(\bm{a}^{(t+1)}\odot\mathrm{e}^{\mathrm{j}\bm{\theta}^{(t+1)}}-\bm{\lambda}^{(t)}/\xi\right)$

7:
$\bm{\lambda}^{(t+1)}=\bm{\lambda}^{(t)}+\xi(\bm{Gd}^{(t+1)}-\bm{a}^{(t+1)}\odot\mathrm{e}^{\mathrm{j}\bm{\theta}^{(t+1)}})$

8: $t=t+1$

9: end while

10: return $\bm{d}^{(t)}$

Next, $\mathring{u}_{\mathrm{syn},n}^{m}(\bm{r}_{\mathrm{S}})$ is represented
by a weighted sum of the expansion coefficients of $g_{l}(\bm{r})$ around
$\bm{r}_{\mathrm{S}}$, i.e., $\mathring{g}_{l,n}^{m}(\bm{r}_{\mathrm{S}})$,
with the weight of $d_{l}$ [9] as

$\displaystyle\mathring{u}_{\mathrm{syn},n}^{m}(\bm{r}_{\mathrm{S}})=\sum_{l=1}^{L}d_{l}\mathring{g}_{l,n}^{m}(\bm{r}_{\mathrm{S}}).$
(11)

By using the translation operator $\hat{S}_{\nu,n}^{\mu,m}$,
$\mathring{g}_{l,n}^{m}(\bm{r}_{\mathrm{S}})$ is represented by using the
expansion coefficients around the secondary source position $\bm{r}_{l}$ as

$\displaystyle\mathring{g}_{l,n}^{m}(\bm{r}_{\mathrm{S}})=\sum_{\nu,\mu}\mathring{g}_{l,\nu}^{\mu}(\bm{r}_{l})\hat{S}_{\nu,n}^{\mu,m}(\bm{r}_{\mathrm{S}}-\bm{r}_{l}),$
(12)

when $\|\bm{r}-\bm{r}_{\mathrm{S}}\|>\|\bm{r}_{\mathrm{S}}-\bm{r}_{l}\|$.
Here, the translation operator $\hat{S}_{\nu,n}^{\mu,m}$ is defined as [21]

$\displaystyle\hat{S}_{\nu,n}^{\mu,m}(\bm{r})=$ $\displaystyle
4\pi\mathrm{j}^{n-\nu}\sum_{q=|\nu-n|}^{\nu+n}\mathrm{j}^{q}(-1)^{\mu}j_{q}(kr)Y_{q}^{m-\mu}(\theta,\phi)^{\ast}$
$\displaystyle\quad\quad\cdot\mathcal{G}(\nu,\mu;n,-m;q).$ (13)

By substituting (11) and (12) into (10), $\mathcal{E}_{\mathrm{dir}}(\bm{d})$
is represented by a quadratic form of $\bm{d}$ as

$\displaystyle\mathcal{E}_{\mathrm{dir}}(\bm{d})=\bm{d}^{\mathsf{H}}\bm{A}\bm{d}$
(14)

with

$\displaystyle\bm{A}=-\frac{R_{\mathrm{S}}^{2}}{2\rho
c}\mathrm{Im}\left[\sum_{n,m}\tilde{w}_{n}^{m}(\bm{r}_{\mathrm{S}})\bm{B}_{n}^{m}\right]$
(15)

and

$\displaystyle(\bm{B}_{n}^{m})_{l_{1},l_{2}}=\sum_{p,q}\mathring{g}_{l_{1},p}^{q}(\bm{r}_{l_{1}})^{\ast}\sum_{p^{\prime},q^{\prime}}\mathring{g}_{l_{2},p^{\prime}}^{q^{\prime}}(\bm{r}_{l_{2}})$
$\displaystyle\quad\quad\cdot\sum_{n,m}\sum_{n^{\prime}=|n-\nu|}^{n+\nu}\hat{S}_{p,n^{\prime}}^{q,(m+\mu)}(\bm{r}_{\mathrm{S}}-\bm{r}_{l_{1}})^{\ast}\hat{S}_{p^{\prime},n}^{q^{\prime},m}(\bm{r}_{\mathrm{S}}-\bm{r}_{l_{2}})$
$\displaystyle\quad\quad\cdot
h_{n}(kR_{\mathrm{S}})h_{n^{\prime}}^{\prime}(kR_{\mathrm{S}})^{\ast}\mathcal{G}(n,m;\nu,\mu;n^{\prime}).$
(16)

Here, $\mathrm{Im}[\cdot]$ denotes the imaginary part of a complex number. It
is necessary to truncate the expansion order to compute the summations in (15)
and (16) in practice, which is a difference from the computation of $\bm{A}$
for $\mathcal{E}_{\mathrm{uni}}(\bm{d})$ based on an analytical representation
given in [18]. It is preferable that $w$ can be represented by low-order
coefficients of $\tilde{w}_{n}^{m}$ to reduce the computational cost of
$\bm{A}$.

###  3.2 Application to Amplitude Matching

As discussed in Sect. 2, $\mathcal{E}_{\mathrm{dir}}(\bm{d})$ formulated in
(14) can be used in various sound field synthesis problems. As an example, we
demonstrate an application of the directionally weighted exterior radiation
suppression in amplitude matching.

The control points are set at $\\{\bm{r}_{i}\\}_{i=1}^{I}$ over $\Omega$. The
desired amplitude at the control points is denoted as
$\bm{a}_{\mathrm{des}}\in\mathbb{R}_{\geq 0}^{I}$. The transfer function
matrix between the secondary sources and control points is defined as
$\bm{G}\in\mathbb{C}^{I\times L}$. The goal is to synthesize the desired
amplitude $\bm{a}_{\mathrm{des}}$ at the control points while suppressing the
exterior radiation power outside the secondary sources. Thus, the cost
function is formulated as

$\displaystyle\mathop{\rm
minimize}\limits_{\bm{d}\in\mathbb{C}^{L}}\left\||\bm{Gd}|-\bm{a}_{\mathrm{des}}\right\|^{2}+\gamma\bm{d}^{\mathsf{H}}\bm{A}\bm{d}+\alpha\|\bm{d}\|^{2},$
(17)

where $|\cdot|$ represents the element-wise absolute value. Since the cost
function in (17) is neither convex nor differentiable, (17) does not have a
closed-form solution. In a similar manner to amplitude matching without
exterior radiation suppression [16], an algorithm based on alternating
direction method of multipliers (ADMM) for solving (17) can be derived, owing
to the simple form of $\mathcal{E}_{\mathrm{dir}}(\bm{d})$, as in Algorithm 1.
Here, $\bm{\lambda}^{(t)}\in\mathbb{C}^{I}$ is the Lagrange multiplier,
$\xi>0$ is the penalty parameter, and $t$ denotes the iteration index. Since
(17) is formulated for the single frequency, the filter to obtain the driving
signal in the time domain can be unnecessarily long. It is also possible to
apply the differential-norm penalty proposed in [16] for $\mathcal{R}(\bm{d})$
to make the filter length small by inducing smoothness of $\bm{d}$ between
frequency bins.

##  4 Experiments

Fig. 2: Experimental setup. Blue dots and green cylinder represent
loudspeakers and target region $\Omega$, respectively.

We performed numerical experiments to evaluate the proposed exterior radiation
suppression in amplitude matching. We compare amplitude matching (AM), AM with
uniformly weighted exterior radiation suppression (AM-Rad), and AM with
directionally weighted exterior radiation suppression (AM-Rad-Dir).

Fig. 2 shows the experimental setup. The target region $\Omega$ was a
cylindrical region of radius $0.2~{}\mathrm{m}$ and height $0.1~{}\mathrm{m}$.
147 control points were uniformly distributed over $\Omega$ at intervals of
$0.05~{}\mathrm{m}$. 6 loudspeakers were equiangularly placed on 4 circles of
radii $0.453~{}\mathrm{m}$ and $0.653~{}\mathrm{m}$ at $z=\pm
0.2~{}\mathrm{m}$. Therefore, the total number of loudspeakers was 24. Each
loudspeaker was assumed to be a point source.

The directional weighting function $w$ was defined as

$\displaystyle w(\theta,\phi)=1+\cos\phi\sin\theta,$ (18)

which allows exterior radiation in the negative direction of the $x$-axis and
suppresses it in the positive direction. The matrix $\bm{A}$ in (15) was
computed by truncating the expansion coefficients up to $\lceil
kR_{\mathrm{S}}\rceil$ with $R_{\mathrm{S}}=0.8~{}\mathrm{m}$. The desired
amplitude at the control points $\bm{a}_{\mathrm{des}}$ was set to $1$, i.e.,
uniform amplitude over $\Omega$. The time-domain driving signal was obtained
by using the differential-norm penalty. The parameter for the differential-
norm penalty was empirically set at $10$. The balancing parameter $\gamma$ was
set at $2.0\times 10^{4}$, which was the largest value so that the interior
synthesis accuracy was kept high below $1000~{}\mathrm{Hz}$.

As evaluation measures, we define the mean square error of the amplitude
synthesized at the control points for each frequency ($\mathrm{MSE}(\omega)$)
and exterior radiation power in the suppression directions
($P_{\mathrm{rad}}$) as

$\displaystyle\mathrm{MSE}(\omega)=\frac{1}{M}\left\|\bm{a}_{\mathrm{syn}}(\omega)-\bm{a}_{\mathrm{des}}(\omega)\right\|^{2}$
(19) $\displaystyle
P_{\mathrm{rad}}(\omega)=\int_{\phi_{1}}^{\phi_{2}}\mathrm{d}\phi\int_{\theta_{1}}^{\theta_{2}}\mathrm{d}\phi
I_{r}(\theta,\phi)r^{2}\sin\theta,$ (20)

where the range of the suppression direction is set by
$(\theta_{1},\theta_{2})=(0,\pi)$ and $(\phi_{1},\phi_{2})=(-\pi/2,\pi/2)$,
i.e., the directions of $w>1$.

Fig. 3: MSE of AM, AM-Rad, and AM-Rad-Dir with respect to frequency. Fig. 4:
$P_{\mathrm{rad}}$ of AM, AM-Rad, and AM-Rad-Dir with respect to frequency.

$\mathrm{MSE}(\omega)$ and $P_{\mathrm{rad}}(\omega)$ are plotted with respect
to frequency in Figs. 3 and 4. Although $\mathrm{MSE}(\omega)$ of AM-Rad and
AM-Rad-Dir was higher than that of AM, their $P_{\mathrm{rad}}(\omega)$ was
significantly smaller than that of AM. Comparing AM-Rad and AM-Rad-Dir,
$P_{\mathrm{rad}}(\omega)$ of AM-Rad-Dir was smaller than that of AM-Rad below
$400~{}\mathrm{Hz}$, owing to the directional weighting. Furthermore,
$\mathrm{MSE}(\omega)$ of AM-Rad-Dir was smaller than that of AM-Rad above
$600~{}\mathrm{Hz}$, which will be owing to the relaxation of the constraint
on the exterior radiation power by directional weighting.

(a) AM

(b) AM-Rad (c) AM-Rad-Dir

Fig. 5: Power distribution of synthesized sound field on $x$-$y$-plane at
$z=0$.

Fig. 5 shows the power distribution of the synthesized sound field up to
$1000~{}\mathrm{Hz}$ on $x$-$y$-plane at $z=0$. The exterior radiation of AM
was particularly larger than that of AM-Rad and AM-Rad-Dir because the
exterior radiation suppression is not taken into consideration in AM. In AM-
Rad, the exterior radiation was suppressed in all directions. Because of the
design of the directional weighting function $w$, the exterior radiation in
the negative $x$ direction was larger than that in the positive $x$ direction
in AM-Rad-Dir.

##  5 Conclusion

We proposed a method for synthesizing the desired sound field while
suppressing the exterior radiation power with directional weighting. The
proposed directionally weighted exterior radiation suppression enables
relaxing the constraint on the exterior radiation compared with uniformly
weighted one, thus the influence on the interior synthesis accuracy owing to
the constraint on the exterior radiation can be mitigated. Since the
directionally weighted exterior radiation power is represented by a quadratic
form of the driving signal of the secondary sources, it is simply applicable
to various sound field synthesis methods. In the numerical experiments using
amplitude matching, the exterior radiation was significantly suppressed in the
prioritized directions while maintaining the interior synthesis accuracy high.

##  6 Acknowledgment

This work was supported by JST FOREST Program Grant Number JPMJFR216M and JSPS
KAKENHI Grant Number 22H03608.

## References

  * [1] A. J. Berkhout, D. de Vries, and P. Vogel,  “Acoustic control by wave field synthesis,”  J. Acoust. Soc. Amer., vol. 93, no. 5, pp. 2764–2778, 1993.
  * [2] S. Spors, R. Rabenstein, and J. Ahrens,  “The theory of wave field synthesis revisited,”  in Proc. 124th AES Conv., Amsterdam, Netherlands, Oct. 2008.
  * [3] 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, Mar. 2003.
  * [4] M. A. Poletti,  “Three-dimensional surround sound systems based on spherical harmonics,”  J. Audio Eng. Soc., vol. 53, no. 11, pp. 1004–1025, 2005.
  * [5] J. Ahrens and S. Spors,  “An analytical approach to sound field reproduction using circular and spherical loudspeaker distributions,”  Acta Acust. united Ac., vol. 94, pp. 988–999, 2008.
  * [6] Y. J. Wu and T. D. Abhayapala,  “Theory and design of soundfield reproduction using continuous loudspeaker concept,”  IEEE Trans. Audio, Speech, Lang. Process., vol. 17, no. 1, pp. 107–116, 2009.
  * [7] 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., vol. 21, no. 4, pp. 685–696, 2013.
  * [8] O. Kirkeby and P. A. Nelson,  “Reproduction of plane wave sound fields,”  J. Acoust. Soc. Amer., vol. 94, no. 5, pp. 2992–3000, 1993.
  * [9] 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., vol. 27, no. 12, pp. 1852–1867, 2019.
  * [10] S. Koyama, K. Kimura, and N. Ueno,  “Weighted pressure and mode matching for sound field reproduction: Theoretical and experimental comparisons,”  J. Audio Eng. Soc., vol. 71, no. 4, pp. 173–185, 2023.
  * [11] 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., vol. 29, pp. 3052–3063, 2021.
  * [12] 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), Singapore, May 2022, pp. 1056–1060.
  * [13] T. Betlehem, W. Zhang, M. A. Poletti, and T. D. Abhayapala,  “Personal sound zones: Delivering interface-free audio to multiple listeners,”  IEEE Signal Process. Mag., vol. 32, no. 2, pp. 81–91, 2015.
  * [14] Y. J. Wu and T. D. Abhayapala,  “Spatial multizone soundfield reproduction: theory and design,”  IEEE Trans. Audio, Speech, Lang. Process., vol. 19, no. 6, pp. 1711–1720, 2010.
  * [15] J.-W. Choi and Y.-H. Kim,  “Generation of an acoustically bright zone with an illuminated region using multiple sources,”  J. Acoust. Soc. Amer., vol. 111, no. 4, pp. 1695–1700, 2002.
  * [16] T. Abe, S. Koyama, N. Ueno, and H. Saruwatari,  “Amplitude matching for multizone sound field control,”  IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 31, pp. 656–669, 2023.
  * [17] M. A. Poletti, T. D. Abhayapala, and P. Samarasinghe,  “Interior and exterior sound field control using two dimensional higher-order variable-directivity sources,,”  .
  * [18] N. Ueno, S. Koyama, and H. Saruwatari,  “Sound field reproduction with exterior radiation cancellation using analytical weighting of harmonic coefficients,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), Calgary, Canada, Apr. 2018, pp. 466–470.
  * [19] K. Arikawa, S. Koyama, and H. Saruwatari,  “Kernel-interpolation-based spatial active noise control with exterior radiation suppression,”  in Proc. Int. Congr. Acoust. (ICA), Gyeongju, Republic of Korea, Oct. 2022.
  * [20] T. Kojima, K. Arikawa, S. Koyama, and H. Saruwatari,  “Multichannel active noise control with exterior radiation suppression based on Riemannian optimization,”  in Proc. European Signal Process. Conf. (EUSIPCO), Helsinki, Finland, Sep. 2023.
  * [21] P. A. Martin,  Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles.,  Cambridge University Press, New York, NY, 2006.


#  Query-Age-Optimal Scheduling under Sampling and Transmission Constraints

Abolfazl Zakeri, Mohammad Moltafet, Markus Leinonen, and Marian Codreanu  A.
Zakeri and M. Leinonen are with CWC-RT, University of Oulu, Finland, e-mail:
{abolfazl.zakeri<EMAIL_ADDRESS>M. Moltafet is with Department of
Electrical and Computer Engineering University of California Santa Cruz,
email<EMAIL_ADDRESS>M. Codreanu is with Department of Science and
Technology, Linköping University, Sweden, e-mail<EMAIL_ADDRESS>This
research has been financially supported by the Infotech Oulu, the Academy of
Finland (grant 323698), and 6G Flagship program (grant 346208). The work of M.
Leinonen has also been financially supported in part by the Academy of Finland
(grant 340171).

###### Abstract

This letter provides query-age-optimal joint sampling and transmission
scheduling policies for a heterogeneous status update system, consisting of a
stochastic arrival and a generate-at-will source, with an unreliable channel.
Our main goal is to minimize the average query age of information (QAoI)
subject to average sampling, average transmission, and per-slot transmission
constraints. To this end, an optimization problem is formulated and solved by
casting it into a linear program. We also provide a low-complexity near-
optimal policy using the notion of weakly-coupled constrained Markov decision
processes. The numerical results show up to $32\%$ performance improvement by
the proposed policies compared with a benchmark policy.

##  I Introduction

The age of information (AoI) has been proposed to characterize the information
freshness in status update systems [1]. The AoI is defined as the time elapsed
since the latest received status update packet was generated [1, 2].
Minimizing the AoI is mainly associated with sampling and scheduling
optimization. However, in some cases, sampling cannot be controlled [3], e.g.,
a sensor takes a sample whenever it harvests enough energy. At the same time,
controllable and uncontrollable111Controllable source refers to the generate-
at-will source that generates updates upon request, while uncontrollable
source refers to a source that generates updates randomly and cannot be
prompted to produce them. sources can coexist and utilize shared resources.
Besides, sampling and transmission of a status update incur different costs
involving distinct non-transferable resources, leading to individual trade-
offs between each resource usage and freshness. These motivate us to address a
freshness optimization problem in a heterogeneous status update system,
consisting of controllable and uncontrollable sources, subject to individual
sampling and transmission constraints.
We consider two sources, one generate-at-will and one stochastic arrival, with
a transmitter communicating with a remote monitor over an error-prone channel
(see Fig. 1). The transmitter is in charge to sample and transmit updates to
the monitor. The freshness at the monitor is captured by a generalized form of
the AoI, the query AoI [4] (QAoI).

We formulate an optimization problem aiming to minimize the sum (time) average
QAoI at the monitor subject to average sampling, average transmission, and
per-slot transmission constraints. We then develop optimal (joint sampling and
scheduling) polices by solving the problem via casting it as a constrained
Markov decision process (CMDP) problem and subsequently, as a linear program
(LP). To strike a balance between complexity and performance, a near-optimal
low-complexity policy is also developed. The following summarizes the main
contributions of the letter:
$\bullet$ We provide a query-age-optimal policy for a heterogeneous status
update system under average sampling and transmission constraints. To this
end, we use the LP approach to solve a multi-constraint CMDP problem.
$\bullet$ We devise a near-optimal low-complexity policy. To do so, we cast
and solve a weakly coupled CMDP problem [5] and propose a dynamic truncation
algorithm.
$\bullet$ Finally, simulation evaluations are conducted that verify the
effectiveness of the proposed policies.

Related work: Recently, the information freshness has gained much attention
from the sampling and/or scheduling optimization perspectives by using the AoI
metric in, e.g., [6, 7, 8, 9, 10, 11, 12, 13], and using the query/request-
based AoI metric in, e.g., [14, 4, 15]. Assuming the generate-at-will model,
the sampling times in [6, 8, 7, 9, 11, 10, 12, 4, 15], and the sampling rate
in [6] were optimized. The work [9] considered a limit on the sampling
frequency. The works [7, 4] considered the transmission cost, which limits the
average number of transmissions [7]. Nevertheless, in some applications,
sampling of the monitored process and transmission of an update both incur
significant, yet different costs [12]. However, only a few works [6, 8, 12,
13] studied joint sampling and scheduling control with sampling and
transmission costs. These models incorporated a weighted sum of the costs. In
contrast, we incorporate the costs separately and consider two individual
constraints, one on the average sampling cost and one on the average
transmission cost. This is motivated by the fact that generally, the sampling
and transmission costs might need to be imposed on the design by distinct and
non-transferable budgets, i.e., the budgets on different types of resources
(e.g., CPU and spectrum), so that savings in one budget do not create a
surplus to be used in the other one. Consequently, casting the underlying
optimization problem as a CMDP inevitably leads to a multi-constraint CMDP
problem. Solving such a problem is in general more challenging than solving
the most studied single-constraint CMDP problems in, e.g., [7, 10, 6, 12].

In summary, the main distinctive contribution of this letter is to consider
distinct and non-transferable sampling cost and transmission cost constraints,
whereas the prior works have considered a single constraint to account for the
both costs. To that extent, our provided solutions can be applied to multi-
constraint CMDP problems, whereas most prior works are limited to a single-
constraint CMDP problem.

##  II System Model and Problem Formulation

We consider a heterogeneous status update system consisting of two independent
sources: one random arrival (Source 1) and one generate-at-will (Source 2), a
buffer-aided transmitter, and a monitor, as depicted in Fig. 1. A discrete-
time system with unit time slots $t\in\\{0,1,\dots\\}$ is considered. The
sources are indexed by $i\in\\{1,2\\}$. The status update packets of Source 1
are generated according to the Bernoulli distribution with parameter $\mu$.
To maintain the freshest information currently available at the transmitter,
the most recently arrived update of Source 1 and the most recently sampled
update of Source 2 are stored in the transmitter’s buffers of size one packet
per source. Notice that considering a buffer of size one packet for each
source is sufficient in our system, because storing and transmitting the
outdated packets does not improve the AoI.

######  Remark 1.

Two sources are considered for the clarity of the presentation. However, the
formulated problem and solutions can be extended to an arbitrary number of
sources. The extension can be obtained by introducing a set of random arrival
sources $\mathcal{L}$ and a set of generate-at-will sources $\mathcal{K}$,
then replacing $i=1$ by $l,\,\forall\,l\in\mathcal{L}$ and $i=2$ by
$k,\,\forall\,k\in\mathcal{K}$ in the formulas. We will examine the impact of
the number of sources on the system performance in Sec. V (see Fig. 4).

Figure 1:  The considered heterogeneous status update system.

To capture the randomness of the wireless link between the transmitter and the
monitor, we assume the link is error-prone, i.e., the reception of updates by
the monitor are successful with a constant (over slots) probability
${p\in(0,1]}$. Moreover, we assume that at each slot, at most one packet can
be transmitted over the link, which stands for the per-slot transmission
constraint in the system. Besides, we assume that each transmission occupies
one slot, and perfect feedback (i.e., instantaneous and error-free) is
available for the link.
Decision Variables: We define two decision variables: one for the scheduling
at the transmitter and the other for the sampling decision for Source 2. Let
$\alpha_{i}(t)\in\\{0,1\\}$, $i\in\\{1,2\\}$, denote the transmission decision
of source $i$ at the transmitter in slot $t$, where $\alpha_{i}(t)=1$ means
that the transmitter sends the packet of source $i$, and $\alpha_{i}(t)=0$
otherwise. Moreover, let $\beta(t)\in\\{0,1\\}$ denote the sampling decision
of Source 2: $\beta(t)=1$ indicates that the transmitter takes a new sample of
Source 2 in slot $t$, and $\beta(t)=0$ otherwise.
Age of Information: First, we make a common assumption (see e.g., [12, 8] and
references therein) that the AoI values are upper-bounded by a sufficiently
large value $N$. Besides tractability, this accounts for the fact that when
the status information becomes excessively stale by reaching $N$, the time-
critical end application would not be affected if counting further. Let
$\theta_{i}(t)$ be the AoI of source $i$ at the transmitter in slot $t$, and
$\delta_{i}(t)$ be the AoI of source $i$ at the monitor in slot $t$. Let
$q_{i}(t)\in\\{0,1\\},\,i\in\\{1,2\\},$ be an indicator for source $i$, where
$q_{i}(t)=1$ means that the packet of source $i$ is transmitted and
successfully received in slot $t$, and ${q_{i}(t)=0}$ otherwise.

Assuming that the packet arrivals occur at the beginning of slots, the
evolution of the ages for $t=1,2,\dots$, with the initial values
$\theta_{i}(0)=0$ and $\delta_{i}(0)=0$, are given as

$\begin{array}[]{ll}&\theta_{1}(t+1)=\left\\{\begin{array}[]{ll}0,&\text{if a
packet arrives in $t+1$},\\\
\left[\theta_{1}(t)+1\right]^{N},&\text{otherwise},\end{array}\right.\\\
&\theta_{2}(t+1)=\left\\{\begin{array}[]{ll}0,&\text{if}\leavevmode\nobreak\
\leavevmode\nobreak\ \beta(t+1)=1,\\\
\left[\theta_{2}(t)+1\right]^{N},&\text{otherwise}.\end{array}\right.\\\
&\delta_{1}(t+1)=\left\\{\begin{array}[]{ll}\left[\theta_{1}(t)+1\right]^{N},&\leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ q_{1}(t)=1,\\\
\left[\delta_{1}(t)+1\right]^{N},&\leavevmode\nobreak\
\text{otherwise},\end{array}\right.\\\
&\delta_{2}(t+1)=\left\\{\begin{array}[]{ll}1,&\text{if}\leavevmode\nobreak\
q_{2}(t)=1,\beta(t)=1,\\\
\left[\theta_{2}(t)+1\right]^{N},&\text{if}\leavevmode\nobreak\
q_{2}(t)=1,\beta(t)\neq 1,\\\
\left[\delta_{2}(t)+1\right]^{N},&\text{otherwise},\end{array}\right.\end{array}$
(1)

where ${\left[.\right]^{N}}\triangleq\min(.\,,N)$.
We consider that the AoI at the monitor becomes important (only) in the
instants when the updates are actually needed/used by the time-critical end
application. We capture this by the QAoI [4] by defining a per-source query
flag denoted by ${r_{i}(t)\in\\{0,1\\},\,i\in\\{1,2\\}},$ where $r_{i}(t)=1$
means that the monitor queries the information of source $i$ in slot $t$, and
$r_{i}(t)=0$ otherwise. Accordingly, the QAoI for source $i$ in slot $t$ is
defined as $r_{i}(t)\delta_{i}(t)$. We assume that the query flags are
generated according to the binary Markov chains, as depicted in Fig. 1, where
$\rho_{i}$ and $\bar{\rho}_{i}$ are, respectively, the self-transition
probabilities of State $1$ and State $0$.

We denote the sum (discounted) average QAoI at the monitor by $\bar{\delta}$,
defined by

$\begin{array}[]{ll}&\displaystyle\bar{\delta}\triangleq\textstyle\bar{\lambda}\sum_{t=1}^{\infty}{\lambda^{t-1}}\mathbb{E}\left\\{\sum_{i}r_{i}(t)\delta_{i}(t)\right\\},\end{array}$
(2)

where $\lambda\in(0,1)$ is the discount factor, $\bar{\lambda}=(1-\lambda)$ is
the normalization factor [16, Remark 2.1], and $\mathbb{E}\\{\cdot\\}$ is the
expectation with respect to the random channel, the packet arrival process,
and the (possibly random) decision variables $\alpha_{i}(t)$ and $\beta(t)$
made in reaction to the past and current AoI and query values. We also define
the (discounted) average number of transmissions, denoted by $\bar{\alpha}$,
and the (discounted) average number of sampling actions, denoted by
$\bar{\beta}$, as

$\displaystyle\displaystyle\bar{\alpha}\triangleq\bar{\lambda}\textstyle\sum_{t=1}^{\infty}\lambda^{t-1}\mathbb{E}\left\\{\sum_{i}\alpha_{i}(t)\right\\},\leavevmode\nobreak\
\displaystyle\bar{\beta}\triangleq\bar{\lambda}\textstyle\sum_{t=1}^{\infty}\lambda^{t-1}\mathbb{E}\left\\{\beta(t)\right\\}.$

Our main goal is to solve the following problem222 As proposing the CMDP
approach to the main problem (3), the provided solutions can be extended to
the average (expected) case (i.e., $\lambda=1$) if the unichain structure [17,
Sec. 8.3.1] exits. The unichain structure, however, does not exist for our
problem since one can construct a deterministic policy which induces a Markov
chain with two recurrent classes. Thus, for the average case, finding a
solution could be, in general, difficult [18, Sec. 9.2.6].:

$\displaystyle\underset{\\{\\{\alpha_{i}(t)\\}_{i=1}^{2},\beta(t)\\}_{t=1,2,\ldots}}{\mbox{minimize}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle\bar{\delta}\vspace{-2em}$ (3a) subject to
$\displaystyle\bar{\alpha}\leq\Gamma^{\mathrm{tr}},$ (3b)
$\displaystyle\bar{\beta}\leq\Gamma^{\mathrm{sm}},$ (3c)
$\displaystyle\textstyle\sum_{i}\alpha_{i}(t)\leq 1,\leavevmode\nobreak\
\forall\,t,$ (3d) $\displaystyle\beta(t)\leq\alpha_{2}(t),\leavevmode\nobreak\
\forall\,t,\vspace{-1.2 em}$ (3e)

with the decision variables
${\\{\\{\alpha_{i}(t)\\}_{i=1}^{2},\beta(t)\\}_{t=1,2,\ldots}}$, where
${\Gamma^{\mathrm{tr}}\in(0,1]}$ and ${\Gamma^{\mathrm{sm}}\in(0,1]}$ are
limits on the average number of transmissions and sampling actions in the
system. The per-slot constraint (3d) limits that at most one packet can be
transmitted over the link in each slot. The per-slot constraint (3e) ensures
the sampled update in slot $t$ is transmitted at the same slot. Note that
constraint (3e) is not obligatory for problem (3), but it eliminates a sub-
optimal action, i.e., taking a new sample without its concurrent transmission,
which simplifies solving the CMDP problem (introduced next).

Next, we provide an optimal policy to problem (3).

##  III  Optimal Policy

To solve problem (3), we first recast it into a CMDP problem and then
(optimally) solve it via solving its equivalent LP [16].

###  III-A CMDP Formulation of Problem (3)

The CMDP is specified by the following elements:
$\bullet$ State: We define the state in slot $t$ by
${\mathbb{s}(t)\triangleq\big{(}r_{1}(t),\theta_{1}(t),\delta_{1}(t),r_{2}(t),\theta_{2}(t),\delta_{2}(t)\big{)}}$.
We denote the state space by $\mathcal{S}$ which is a finite set.
$\bullet$ Action:  Actions determine the transmission and sampling decisions
while incorporating constraints (3d) and (3e). Thus, we define the action
taken in slot $t$ by ${a(t)\in\\{0,1,2_{\mathrm{r}},2_{\mathrm{n}}\\}}$, where
${a(t)=0}$ means that the transmitter stays idle, ${a(t)=1}$ means that the
packet of Source 1 is transmitted/re-transmitted, ${a(t)=2_{\mathrm{r}}}$
means that the packet of Source 2 is re-transmitted, and
${a(t)=2_{\mathrm{n}}}$ means that a new sample of Source 2 is transmitted.
Also, $\mathcal{A}$ is the action space of size $|\mathcal{A}|=4$. Actions are
determined by a policy $\pi$, which is a (possibly randomized) mapping from
$\mathcal{S}$ to $\mathcal{A}$.
$\bullet$ State Transition Probabilities: We denote the state transition
probability from state $\mathbb{s}$ to next state $\mathbb{s}^{\prime}$ under
an action ${a}$ by $\mathcal{P}_{\mathbb{s}\mathbb{s}^{\prime}}(a)$. Since the
evolution of the AoIs in (LABEL:Eq_AoIDy) and the query processes are
independent among the sources, the transition probability can be written as
${\mathcal{P}_{\mathbb{s}\mathbb{s}^{\prime}}(a)=\prod_{i}\Pr\\{r_{i}^{\prime}\,|\,r_{i}\\}\Pr\\{\mathbb{\underline{s}}_{i}^{\prime}\,\big{|}\,\mathbb{\underline{s}}_{i},a\\}}$,
where $\Pr\\{r_{i}^{\prime}\,|\,r_{i}\\}$ is given by (20), and
$\Pr\\{\mathbb{\underline{s}}_{i}^{\prime}\,\big{|}\,\mathbb{\underline{s}}_{i},a\\}$
is given by (4), in which
$\mathbb{\underline{s}}_{i}\triangleq(\theta_{i},\delta_{i})$ contains the
current age values, as a part of the current state, associated with source
$i$,
$\mathbb{\underline{s}}_{i}^{\prime}\triangleq(\theta_{i}^{\prime},\delta^{\prime}_{i})$
contains the next age values, as a part of the next state, associated with
source $i$, and $\bar{\mu}=1-\mu$.

$\displaystyle\Pr\\{\mathbb{\underline{s}}_{1}^{\prime}\,\big{|}\,\mathbb{\underline{s}}_{1},a\\}=$
(4) $\displaystyle\left\\{\begin{array}[]{ll}\mu,&a\neq 1;\leavevmode\nobreak\
\theta^{\prime}_{1}=0,\leavevmode\nobreak\
\delta^{\prime}_{1}=\left[\delta_{1}+1\right]^{N},\\\ \bar{\mu},&a\neq
1;\leavevmode\nobreak\
\theta^{\prime}_{1}=\left[\theta_{1}+1\right]^{N},\,\delta^{\prime}_{1}=\left[\delta_{1}+1\right]^{N},\\\
\mu p,&a=1;\leavevmode\nobreak\ \theta^{\prime}_{1}=0,\leavevmode\nobreak\
\delta^{\prime}_{1}=\left[\theta_{1}+1\right]^{N},\\\
\mu(1-p),&a=1;\leavevmode\nobreak\ \theta^{\prime}_{1}=0,\leavevmode\nobreak\
\delta^{\prime}_{1}=\left[\delta_{1}+1\right]^{N},\\\
\bar{\mu}p,&a=1;\leavevmode\nobreak\
\theta^{\prime}_{1}=\left[\theta_{1}+1\right]^{N},\leavevmode\nobreak\
\delta^{\prime}_{1}=\left[\theta_{1}+1\right]^{N},\\\
\bar{\mu}(1-p),&a=1;\leavevmode\nobreak\
\theta^{\prime}_{1}=\left[\theta_{1}+1\right]^{N},\leavevmode\nobreak\
\delta^{\prime}_{1}=\left[\delta_{1}+1\right]^{N},\\\
0&\text{otherwise}.\end{array}\right.$ (12)
$\displaystyle\Pr\\{\mathbb{\underline{s}}_{2}^{\prime}\leavevmode\nobreak\
\big{|}\leavevmode\nobreak\ \mathbb{\underline{s}}_{2},a\\}=$
$\displaystyle\left\\{\begin{array}[]{ll}1,&a\neq\\{2_{\mathrm{r}},2_{\mathrm{n}}\\};\,\theta^{\prime}_{2}=\left[\theta_{2}+1\right]^{N},\,\delta^{\prime}_{2}=\left[\delta_{2}+1\right]^{N},\\\
p,&a=2_{\mathrm{r}};\leavevmode\nobreak\
\theta^{\prime}_{2}=\left[\theta_{2}+1\right]^{N},\leavevmode\nobreak\
\delta^{\prime}_{2}=\left[\theta_{2}+1\right]^{N},\\\
1-p,&a=2_{\mathrm{r}};\leavevmode\nobreak\
\theta^{\prime}_{2}=\left[\theta_{2}+1\right]^{N},\leavevmode\nobreak\
\delta^{\prime}_{2}=\left[\delta_{2}+1\right]^{N},\\\
p,&a=2_{\mathrm{n}};\leavevmode\nobreak\
\theta^{\prime}_{2}=1,\leavevmode\nobreak\ \delta^{\prime}_{2}=1,\\\
1-p,&a=2_{\mathrm{n}};\leavevmode\nobreak\
\theta^{\prime}_{2}=1,\leavevmode\nobreak\
\delta^{\prime}_{2}=\left[\delta_{2}+1\right]^{N},\\\
0&\text{otherwise}.\end{array}\right.$ (19)
$\Pr\\{r_{i}^{\prime}\,|\,r_{i}\\}=\left\\{\begin{array}[]{ll}\rho_{i},&r_{i}=1;\leavevmode\nobreak\
r^{\prime}_{i}=1,\\\ 1-\rho_{i},&r_{i}=1;\leavevmode\nobreak\
r^{\prime}_{i}=0,\\\ \bar{\rho}_{i},&r_{i}=0;\leavevmode\nobreak\
r^{\prime}_{i}=0,\\\ 1-\bar{\rho}_{i},&r_{i}=0;\leavevmode\nobreak\
r^{\prime}_{i}=1.\end{array}\right.$ (20)

$\bullet$ Cost Functions: The cost functions include: 1) the QAoI cost defined
by $c(\mathbb{s}(t))=\sum_{i}r_{i}(t)\delta_{i}(t)$, 2) the transmission cost
defined by $d^{\mathrm{tr}}(a(t))=\mathds{1}_{\\{a\neq 0\\}}$, and 3) the
sampling cost defined by
$d^{\mathrm{sm}}(a(t))=\mathds{1}_{\\{a=2_{\mathrm{n}}\\}}$, where
$\mathds{1}_{\\{\cdot\\}}$ is an indicator function which equals to $1$ when
the condition in $\\{\cdot\\}$ holds.
By the above definitions, problem (3) can be recast as the following CMDP
problem

$\displaystyle\underset{\pi\in\Pi}{\text{minimize}}\leavevmode\nobreak\ $
$\displaystyle\bar{\lambda}\textstyle\sum_{t=1}^{\infty}\lambda^{t-1}\mathbb{E}^{\pi,\eta}\left\\{c(\mathbb{s}(t))\right\\},$
(21a) subject to
$\displaystyle\bar{\lambda}\textstyle\sum_{t=1}^{\infty}\lambda^{t-1}\mathbb{E}^{\pi,\eta}\left\\{d^{\mathrm{tr}}(a(t))\right\\}\leq\Gamma^{\mathrm{tr}},$
(21b)
$\displaystyle\bar{\lambda}\textstyle\sum_{t=1}^{\infty}\lambda^{t-1}\mathbb{E}^{\pi,\eta}\left\\{d^{\mathrm{sm}}(a(t))\right\\}\leq\Gamma^{\mathrm{sm}},$
(21c)

where $\Pi$ is the set of all stationary randomized policies, and
$\mathbb{E}^{\pi,\eta}$ denotes the expectation when following a policy $\pi$
for a given initial distribution (over the state space) $\eta$.
Next, we provide an optimal policy to problem (21).

###  III-B Linear Programming of the CMDP Problem (21)

Here, we transform problem (21) into an LP. For a notional simplicity, we
denote state $\mathbb{s}$ by $s$. Let
${x(s,a)\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}}$,
${\forall\,s\in\mathcal{S},\leavevmode\nobreak\ \forall\,a\in\mathcal{A}}$, be
defined as [16, Ch. 3]

$\displaystyle\begin{array}[]{ll}x(s,a)\triangleq\bar{\lambda}\sum_{t=1}^{\infty}\lambda^{t-1}\Pr^{\eta}{\\{s(t)=s,\,a(t)=a\\}},\end{array}$
(23)

where $\Pr^{\eta}{\\{s(t)=s,\,a(t)=a\\}}$ denotes the probability of taking
action $a$ at state $s$ in slot $t$ given the initial distribution $\eta$. It
is noteworthy that $x(s,a),\leavevmode\nobreak\
\forall\,s\in\mathcal{S},\leavevmode\nobreak\ \forall\,a\in\mathcal{A},$ can
be interpreted as the long run discounted time that the system is in state $s$
and action $a$ is chosen. Then, the CMDP problem (21) can be transformed into
the following LP [16, Ch. 3]:

$\displaystyle\underset{\mathbb{X}}{\text{minimize}}\leavevmode\nobreak\
\leavevmode\nobreak\ $
$\displaystyle\textstyle\sum_{s\in\mathcal{S}}\sum_{a\in\mathcal{A}}x(s,a)c(s)$
(24a) subject to
$\displaystyle\textstyle\sum_{s\in\mathcal{S}}\sum_{a\in\mathcal{A}}x(s,a)d^{\mathrm{tr}}(a)\leq\Gamma^{\mathrm{tr}},$
(24b)
$\displaystyle\textstyle\sum_{s\in\mathcal{S}}\sum_{a\in\mathcal{A}}x(s,a)d^{\mathrm{sm}}(a)\leq\Gamma^{\mathrm{sm}},$
(24c)
$\displaystyle\textstyle\sum_{s^{\prime}\in\mathcal{S}}\sum_{a\in\mathcal{A}}x(s^{\prime},a)\big{(}\mathds{1}_{\\{s=s^{\prime}\\}}-\lambda\mathcal{P}_{s^{\prime}s}(a)\big{)}$
$\displaystyle=\bar{\lambda}\eta(s),\forall\,s\in\mathcal{S},$ (24d)
$\displaystyle x(s,a)\geq 0,\leavevmode\nobreak\
\forall\,s\in\mathcal{S},\leavevmode\nobreak\ \forall\,a\in\mathcal{A},$ (24e)

with variables $\mathbb{X}=[x(s,a)]$, where
$\eta(s),\,\forall\,s\in\mathcal{S}$, is the probability that the initial
state is state $s$.

By [16, Theorem 3.3], the following theorem relates a solution of the LP (24)
to a solution of problem (21).

######  Theorem 1.

First, the optimal value of the CMDP problem (21) equals to the optimal value
of the LP (24). Moreover, let ${\mathbb{X}^{*}=[x^{*}(s,a)]},$ be a solution
of the LP (24). Then, a stationary randomized policy
$\pi\triangleq\\{f(s,a)\\}$ with

$\displaystyle\begin{array}[]{cc}f(s,a)=\left\\{\begin{array}[]{ll}x^{*}(s,a)/\bar{x}_{s},&\text{if}\leavevmode\nobreak\
\bar{x}_{s}>0,\\\
\text{arbitrary},&\text{otherwise},\end{array}\right.\end{array}$ (28)

is an optimal policy to the CMDP problem (21), where
${\bar{x}_{s}=\textstyle\sum_{a^{\prime}\in\mathcal{A}}x^{*}(s,a^{\prime})}$,
and $f(s,a),\leavevmode\nobreak\ \forall\,s\in\mathcal{S},\leavevmode\nobreak\
\forall\,a\in\mathcal{A},$ is the probability that an action $a$ is chosen at
state $s$.

Complexity analysis: The complexity of finding the optimal policy (determined
by (28)) comes from the complexity of the LP (24), which is exponential in the
number of sources. This is because its number of variables and constraints
vary with the state space size of the CMDP problem (21), which grows
exponentially with the number of sources. Thus, the optimal policy becomes
computationally inefficient when applied for large numbers of sources. To
strike a balance between the performance and complexity, a heuristic low-
complexity policy is provided in the next section.

##  IV Low-Complexity Policy

In this section, we develop a low-complexity policy to the main problem (3)
via the following procedure: 1) we relax the per-slot transmission constraint
(3d), 2) transform the relaxed problem to a weakly coupled CMDP problem which
is then solved by its equivalent LP [5], and 3) propose a dynamic truncation
algorithm to satisfy the transmission constraint (3d).
We relax the per-slot constraint (3d) to the (discounted) time average
constraint
$\bar{\lambda}\sum_{t=1}^{\infty}\lambda^{t-1}\mathbb{E}\\{\textstyle\sum_{i}\alpha_{i}(t)\\}\leq
1$. This constraint then becomes inactive due to constraint (3b). Then, the
main problem (3) without constraint (3d) can be cast as a weakly coupled CMDP
problem. This is because the transition probabilities, the action space, and
the immediate cost functions become independent among the sources; now, the
average transmission constraint (3b) is the only constraint coupling the
sources.
To proceed, we define the state, action, transition probabilities, and cost
functions associated with each source $i$. The state of source $i$ in slot $t$
is ${s_{i}(t)\triangleq\big{(}r_{i}(t),\theta_{i}(t),\delta_{i}(t)\big{)}}$,
with state space $\mathcal{S}_{i}$, which is a finite set. The action of
source $i$ in slot $t$ is denoted by $a_{i}(t)\in\mathcal{A}_{i}$ with action
spaces $\mathcal{A}_{1}=\\{0,1\\}$ and
$\mathcal{A}_{2}=\\{0,2_{\mathrm{r}},2_{\mathrm{n}}\\}$. The state transition
probabilities of source $i$, denoted as
${\mathcal{\tilde{P}}_{s_{i}s_{i}^{\prime}}^{i}(a),\leavevmode\nobreak\
\forall\,s_{i},\,s_{i}^{\prime}\in\mathcal{S}_{i}}$, is given by
$\Pr\\{r_{i}^{\prime}\,|\,r_{i}\\}\Pr\\{\mathbb{\underline{s}}_{i}^{\prime}\,\big{|}\,\mathbb{\underline{s}}_{i},a\\}$,
where $\Pr\\{r_{i}^{\prime}\,|\,r_{i}\\}$ and
$\Pr\\{\mathbb{\underline{s}}_{i}^{\prime}\,\big{|}\,\mathbb{\underline{s}}_{i},a\\}$
were defined in (20) and (4). The QAoI cost of source $i$ is
${c_{i}({s}_{i}(t))=r_{i}(t)\delta_{i}(t)}$. The transmission cost of source
$i$ is ${d_{i}^{\mathrm{tr}}(a_{i}(t))=\mathds{1}_{\\{a_{i}(t)\neq 0\\}}}$.
The sampling cost of Source 2 is
${d^{\mathrm{sm}}(a_{2}(t))=\mathds{1}_{\\{a_{2}(t)=2_{n}\\}}}$.
First, for a presentation simplicity, we remove the subscript $i$ from
$s_{i}(t)$ and $a_{i}(t)$. Then, for each source $i\in\\{1,2\\}$, let
$x_{i}(s,a)\in\mathbb{R}^{|\mathcal{S}_{i}||\mathcal{A}_{i}|},\leavevmode\nobreak\
\forall\,s\in\mathcal{S}_{i},\leavevmode\nobreak\
\forall\,a\in\mathcal{A}_{i},$ be defined as

$\displaystyle\begin{array}[]{ll}x_{i}(s,a)&\triangleq\bar{\lambda}\sum_{t=1}^{\infty}\lambda^{t-1}\Pr_{i}^{\eta_{i}}{\\{s(t)=s,\,a(t)=a\\}},\end{array}$

where $\Pr_{i}^{\eta_{i}}{\\{s(t)=s,\,a(t)=a\\}}$ denotes the probability of
taking action $a\in\mathcal{A}_{i}$ at state $s\in\mathcal{S}_{i}$ in slot $t$
given the initial distribution $\eta_{i}$. Then, the corresponding LP of the
weakly coupled CMDP problem is given by

$\displaystyle\underset{\mathbb{X}_{i}}{\mbox{minimize}}\leavevmode\nobreak\
\leavevmode\nobreak\ $
$\displaystyle\textstyle\sum_{i}\textstyle\sum_{s\in\mathcal{S}_{i}}\sum_{a\in\mathcal{A}_{i}}x_{i}(s,a)c_{i}(s)$
(29a) subject to
$\displaystyle\textstyle\sum_{i}\textstyle\sum_{s\in\mathcal{S}_{i}}\sum_{a\in\mathcal{A}_{i}}x_{i}(s,a)d_{i}^{\mathrm{tr}}(a)\leq\Gamma^{\mathrm{tr}},$
(29b)
$\displaystyle\textstyle\sum_{s\in\mathcal{S}_{2}}\sum_{a\in\mathcal{A}_{2}}x_{2}(s,a)d^{\mathrm{sm}}(a)\leq\Gamma^{\mathrm{sm}},$
(29c)
$\displaystyle\textstyle\sum_{s^{\prime}\in\mathcal{S}_{i}}\sum_{a\in\mathcal{A}_{i}}x_{i}(s^{\prime},a)\big{(}\mathds{1}_{\\{s=s^{\prime}\\}}-\lambda\mathcal{\tilde{P}}^{i}_{s^{\prime}s}(a)\big{)}$
$\displaystyle=\bar{\lambda}\eta_{i}(s),\forall\,i,\leavevmode\nobreak\
\forall\,s\in\mathcal{S}_{i},$ (29d) $\displaystyle x_{i}(s,a)\geq
0,\forall\,i,\leavevmode\nobreak\
\forall\,s\in\mathcal{S}_{i},\leavevmode\nobreak\
\forall\,a\in\mathcal{A}_{i},$ (29e)

with variables $\mathbb{X}_{i}=[x_{i}(s,a)],\leavevmode\nobreak\
i\in\\{1,2\\}$.

Let $\mathbb{X}_{i}^{*}=[x_{i}^{*}(s,a)]$ be a solution of the LP (29). By
results of [5], an optimal stationary randomized policy for each source $i$,
$\pi^{*}_{i}\triangleq\\{f_{i}(s,a)\\}$, is given by

$\begin{array}[]{cc}f_{i}(s,a)=\left\\{\begin{array}[]{ll}x_{i}^{*}(s,a)/\bar{x}_{i,s},&\text{if}\leavevmode\nobreak\
\bar{x}_{i,s}>0,\\\
\text{arbitrary},&\text{otherwise},\end{array}\right.\end{array}$ (30)

where
${\bar{x}_{i,s}=\textstyle\sum_{a^{\prime}\in\mathcal{A}_{i}}x_{i}^{*}(s,a^{\prime})}$
and ${f_{i}(s,a)}$, ${\forall\,s\in\mathcal{S}_{i}}$,
${\forall\,a\in\mathcal{A}_{i}}$, is the probability that an action
$a\in\mathcal{A}_{i}$ is chosen at state $s\in\mathcal{S}_{i}$ for source $i$.
Notably, the per-source policies $\\{\pi_{i}^{*}\\}$ construct a lower-bound
policy to the main problem (3); this lower bound is tight as empirically shown
in Fig. 4.
Having constructed the above weakly coupled CMDP and derived its solution, we
now propose a near-optimal low-complexity policy for the main problem (3) that
operates as follows: (i) it determines an action of each source $i$ in each
slot $t$ according to $\pi^{*}_{i}$, and (ii) having the actions determined,
if there is more than one packet to be transmitted, the packet of source
$\underline{i}$ will be transmitted, where
${\underline{i}=\operatorname*{arg\,max}_{i\in\\{1,2\\}}\\{h_{i}(t)\\}}$ in
which $h_{1}(t)$ and $h_{2}(t)$ are given as

$\displaystyle\begin{array}[]{cc}h_{1}(t)\triangleq
r_{1}(t)(\delta_{1}(t)-\theta_{1}(t)),\leavevmode\nobreak\ h_{2}(t)\triangleq
r_{2}(t)\delta_{2}(t).\end{array}$ (32)

We term the above Step (ii) as a dynamic truncation algorithm, which still
couples the sources. The main idea of the truncation algorithm is to transmit
the packet that is more fresh from the monitor’s perspective, as it has
greater potential to decrease the AoI at the monitor. To capture such relative
freshness, we have used the AoI difference between the monitor and the
transmitter for Source 1, and the AoI at the monitor itself for Source 2.

Complexity analysis: The main complexity of finding the low-complexity policy
comes from solving the LP (29), which is (only) linear in the number of
sources. This is because its number of variables and constraints vary with the
state space size of the weakly-coupled CMDP problem, which grows linearly with
the number of sources.  Thus, the proposed policy has low complexity, and as
numerically shown in the next section, obtains near-optimal performance.

##  V Numerical Results

Here, we assess the performance of the proposed policies. We set $N=20$,
$\rho_{1}=\rho_{2}=0.7$, $\bar{\rho}_{1}=\bar{\rho}_{2}=0.4$, and
$\lambda=0.99999$, and initialize the query flags by the steady-state
probabilities of the Markov chains, unless otherwise stated. For a benchmark,
we consider a greedy-based baseline policy with the decision rule: If
$\bar{D}^{\mathrm{tr}}(t)\leq\Gamma^{\mathrm{tr}}$, find
${\underline{i}=\operatorname*{arg\,max}_{i\in\\{1,2\\}}\\{h_{i}(t)\\}}$;
otherwise, $a(t)=0$. If $\underline{i}=1$, then $a(t)=1$; otherwise, do the
following. If $\bar{D}^{\mathrm{sm}}(t)\leq\Gamma^{\mathrm{sm}}$, then
$a(t)=2_{\mathrm{n}}$; otherwise, $a(t)=2_{\mathrm{r}}$. Here,
$\bar{D}^{\mathrm{tr}}(t)$ and $\bar{D}^{\mathrm{sm}}(t)$ denote,
respectively, the discounted average number of transmissions and sampling
actions until slot $t$, and $h_{i}(t)$ was given by (32).
Fig. 2 shows the sum (discounted) average QAoI performance of the proposed
policies and the baseline policy as a function of the transmission budget
$\Gamma^{\mathrm{tr}}$ (i.e., the limit on the average number of
transmissions). First, Fig. 2 (as well as Fig. 3) reveals that the heuristic
policy has near-optimal performance, i.e., it nearly coincides with the
optimal policy. The figure shows the importance of an optimal trade-off design
between freshness and resource usage, as the performance gap between the
baseline policy and the proposed policies is significant when the transmission
budget is small. Moreover, the figure reveals that the QAoI is more sensitive
to the transmission budget than the sampling budget, as expected.
Fig. 3 examines the impact of the channel reliability $p$ on the sum average
QAoI performance of the different policies for different packet arrival rate
$\mu$ for Source 1. First, the figure shows that the performance gap between
the proposed policies and the baseline policy is large, especially when the
channel reliability is small. The figure also shows that the QAoI considerably
decreases as the channel reliability increases because it increases the chance
that the monitor receives updates more frequently. Moreover, it can be seen
that the QAoI would not change much by increasing the arrival rate under a
good channel condition.

Figure 2:  The sum (discounted) average QAoI versus the transmission budget
$\Gamma^{\mathrm{tr}}$ for different channel reliabilities $p$, where
$\mu=0.6$.  Figure 3:  The sum (discounted) average QAoI versus the channel
reliability $p$, where $\Gamma^{\mathrm{tr}}=0.5$ and
$\Gamma^{\mathrm{sm}}=0.3$.  Figure 4:  The sum (discounted) average AoI and
QAoI versus the total number of sources with equal number of random arrival
and generate-at-will sources, where $N=150$, $\Gamma^{\mathrm{tr}}=0.8$,
$\Gamma^{\mathrm{sm}}=0.3$, and $p=0.99999$. Moreover, $\mu_{i}=0.5$ (for all
the random arrival sources) and $\rho_{i}=\bar{\rho}_{i}=0.7$ (for all
sources).

Fig. 4 shows the impact of the number of sources on the sum average QAoI/AoI
performance of different policies. First, the figure shows that the proposed
low-complexity policy indeed preserves near-optimal performance in a general
multi-source setup, as it nearly coincides with the lower-bound policy.
Moreover, the performance gap between the baseline policy and the proposed
policy increases in the number of sources.

##  VI Conclusion

We provided a query-age-optimal joint sampling and scheduling policy for a
heterogeneous status update system under average sampling, average
transmission, and per-slot transmission constraints. To this end, a CMDP
problem with two average constraints was cast and solved via its equivalent
LP. We also provided a low-complexity policy by casting a weakly-coupled CMDP
problem and solving it via its equivalent LP. Numerical results showed the
effectiveness of the proposed policies, revealing that age-optimal sampling
and scheduling is crucial for resource-constrained status update systems,
where a greedy-based policy is inefficient. Moreover, the near-optimal
performance of the proposed low-complexity policy substantiates that the LP
approach provides a viable solution for practical systems having many
sources/sensors.

## References

  * [1] S. Kaul, R. Yates, and M. Gruteser, “Real-time status: How often should one update?,” in Proc. IEEE Int. Conf. on Computer Commun., pp. 2731–2735, Orlando, FL, USA, Mar. 2012.
  * [2] Y. Sun, I. Kadota, R. Talak, and E. Modiano, “Age of information: A new metric for information freshness,” Synthesis Lectures on Communication Networks, vol. 12, no. 2, pp. 1–224, Dec. 2019.
  * [3] A. Maatouk, M. Assaad, and A. Ephremides, “On the age of information in a csma environment,” IEEE/ACM Trans. Netw., vol. 28, no. 2, pp. 818–831, Apr. 2020.
  * [4] F. ChiarIoTti, J. Holm, A. E. Kalør, B. Soret, S. K. Jensen, T. B. Pedersen, and P. Popovski, “Query age of information: Freshness in pull-based communication,” IEEE Trans. Commun., vol. 70, no. 3, pp. 1606–1622, Mar. 2022.
  * [5] M. Gagrani and A. Nayyar, “Weakly coupled constrained markov decision processes in borel spaces,” in Proc. American Control Conf., pp. 2790–2795, Denver, CO, USA, Jul. 2020.
  * [6] Y. Cao, Y. Teng, M. Song, and N. Wang, “Age-driven joint sampling and non-slot based scheduling for industrial internet of things,” arXiv preprint arXiv:2205.04092, May, 2022.
  * [7] E. T. Ceran, D. Gündüz, and A. György, “Average age of information with hybrid ARQ under a resource constraint,” IEEE Trans. Wireless Commun., vol. 18, no. 3, pp. 1900–1913, Mar. 2019.
  * [8] E. Fountoulakis, M. Codreanu, A. Ephremides, and N. Pappas, “Joint sampling and transmission policies for minimizing cost under AoI constraints,” arXiv preprint arXiv:2103.15450,, Mar. 2021.
  * [9] J. Pan, A. M. Bedewy, Y. Sun, and N. B. Shroff, “Optimizing sampling for data freshness: Unreliable transmissions with random two-way delay,” arXiv preprint arXiv:2201.02929, 2022.
  * [10] A. M. Bedewy, Y. Sun, S. Kompella, and N. B. Shroff, “Optimal sampling and scheduling for timely status updates in multi-source networks,” IEEE Trans. Inf. Theory, vol. 67, no. 6, pp. 4019–4034, Jun. 2021.
  * [11] I. Kadota, A. Sinha, E. Uysal-Biyikoglu, R. Singh, and E. Modiano, “Scheduling policies for minimizing age of information in broadcast wireless networks,” IEEE/ACM Trans. Netw., vol. 26, no. 6, pp. 2637–2650, Dec. 2018.
  * [12] B. Zhou and W. Saad, “Joint status sampling and updating for minimizing age of information in the internet of things,” IEEE Trans. Commun., vol. 67, no. 11, pp. 7468–7482, Nov. 2019.
  * [13] S. Wang, M. Chen, Z. Yang, C. Yin, W. Saad, S. Cui, and H. V. Poor, “Distributed reinforcement learning for age of information minimization in real-time IoT systems,” IEEE J. Sel. Topics Signal Process., vol. 16, no. 3, pp. 501–515, Apr. 2022.
  * [14] F. Li, Y. Sang, Z. Liu, B. Li, H. Wu, and B. Ji, “Waiting but not aging: Optimizing information freshness under the pull model,” IEEE/ACM Trans. Netw., vol. 29, no. 1, pp. 465–478, Feb. 2021.
  * [15] M. E. Ildiz, O. T. Yavascan, E. Uysal, and O. T. Kartal, “Pull or wait: How to optimize query age of information,” arXiv preprint arXiv:2111.02309, Nov. 2021.
  * [16] E. Altman, Constrained Markov Decision Processes.  volume 7. CRC Press, 1999.
  * [17] M. L. Puterman, Markov Decision Processes: Discrete Stochastic Dynamic Programming.  The MIT Press, 1994.
  * [18] L. Kallenberg, Markov decision processes.  Lect. Notes. Univ. Leiden, Oct. 2016.


# A Convolutional Spiking Network for Gesture Recognition in Brain-Computer
Interfaces

Yiming Ai  Department of Engineering
King’s College London
London, United Kindom
<EMAIL_ADDRESS>Bipin Rajendran  Department of Engineering
King’s College London
London, United Kindom
<EMAIL_ADDRESS>
###### Abstract

Brain-computer interfaces are being explored for a wide variety of therapeutic
applications. Typically, this involves measuring and analyzing continuous-time
electrical brain activity via techniques such as electrocorticogram (ECoG) or
electroencephalography (EEG) to drive external devices. However, due to the
inherent noise and variability in the measurements, the analysis of these
signals is challenging and requires offline processing with significant
computational resources. In this paper, we propose a simple yet efficient
machine learning-based approach for the exemplary problem of hand gesture
classification based on brain signals. We use a hybrid machine learning
approach that uses a convolutional spiking neural network employing a bio-
inspired event-driven synaptic plasticity rule for unsupervised feature
learning of the measured analog signals encoded in the spike domain. We
demonstrate that this approach generalizes to different subjects with both EEG
and ECoG data and achieves superior accuracy in the range of 92.74-97.07% in
identifying different hand gesture classes and motor imagery tasks.

###### Index Terms:

Spiking Neural Network, Brain-computer interface, Event-driven plasticity,
K-means clustering

##  I Introduction

Context and motivation. Brain-computer interface (BCI) is a cutting-edge
technology that adopts state-of-the-art signal processing methods and machine
learning algorithms to detect underlying patterns in time-varying signals
measured from the brain. Nerve illnesses due to neurological disorders, nerve
pathway blockage, poliomyelitis or strokes, and muscular impairments like
paralysis, amyotrophy, or accidents all cause both physical and mental
injuries to patients. In this scenario, BCI provides a potential future for
these patients to alleviate or even eliminate the disability and improve their
quality of life by enabling them to use motor imagery to operate computers,
robotic prostheses, and even their own organs[1, 2, 3, 4, 5]. One of the most
crucial components of a reliable BCI technology is the interpretation of
measured biological signals to drive prosthetic devices. How to best correlate
motor imagery – the cognitive process involving the imagination of a movement
without actually executing it, and the actual motor intention – remains a
challenging problem [6, 7, 8].

Figure 1: The proposed overall scheme for spatio-temporal signal analysis from
BCI recordings.

Traditional signal processing methods as well as machine learning algorithms
[9, 10] have made remarkable progress in solving BCI-based behavior
recognition tasks [11], although efficient analysis of neural signals that are
time-varying and multi-dimensional remains a challenging task. Combining the
abilities of convolutional neural networks (CNN) for pattern recognition and
recurrent neural networks (RNN) for sequential signal processing, a hybrid
CNN-RNN architecture has been implemented for surface electromyography
(sEMG)-based gesture recognition [12]. As the size of the measured data
increases, aforementioned deep learning algorithms that use a large parameter
space make real-time predictions with high accuracy and reliability across
different subjects particularly challenging.

Spiking neural networks (SNNs) are machine learning models that mimic the
spike-triggered information processing modality of the brain [13]. They have
some natural advantages in the processing of time-varying signals as spiking
neuron models naturally incorporate recurrent dynamics. They also hold promise
for event-driven learning and adaptation as well as continual learning
[14][15]. They are also widely studied for efficient hardware implementations
using both CMOS and post-CMOS device technologies [16].

Main contributions. In this paper, we develop a hybrid machine-learning
approach with high classification accuracy that could be used for real-time
gesture recognition. Thanks to the small number of parameters involved in the
model, this method is also amenable to implementation in battery-operated
embedded processors.

Our approach treats neural signals belonging to different cortical regions
separately while extracting temporal correlation information inherent to
different time scales. We highlight some key features of our approach:

  * •

We first identify electrode channels that belong to the same neighborhood in
the 3D-euclidean space and consider the measured signals from the cluster
jointly for further processing.

  * •

We apply a temporal contrast coding scheme to translate the measured analog
signals to spike streams that consist of positive and negative spikes.

  * •

Signals from the same cluster are then analyzed using a spatiotemporal
convolution filter whose output drives an integrate-and-fire neuron.

  * •

The kernel weights of the convolution filter are modified in an event-driven
manner inspired by the biological spike-timing-dependent-plasticity (STDP)
rule.

  * •

The weights at the end of the presentation of a signal from all the clusters
are then analyzed using a k-nearest neighbor (KNN) classifier to identify the
gesture class.

A key advantage of our scheme is that the kernel weights of the network evolve
in an event-driven manner based on the history of spikes over a very short
duration (3 most recent time instants), thus avoiding the need to record and
store analog signals for long durations for signal processing.

This approach is used to classify gesture signals from the Stanford ECOG
dataset and GigaScience EEG data set. For the subjects with electrode
placement information, we obtain classification accuracy in the ranges of
$92.74-97.07\%$, surpassing the performance of previously reported work in the
literature [17, 18, 19, 20, 21, 22, 23, 24, 25].

##  II Background

The ECoG signals used in this work are acquired from the public Stanford
dataset library[26]; we limit our study to the two subjects from the dataset
for whom the electrode placement information was available. To assess brain
activity, an array of platinum electrodes were placed on various cortical
regions. The subjects were instructed to move a single finger in response to
instructions displayed on a computer next to them. For instance, when the
computer displayed “thumb”, the patients were asked to move the corresponding
finger. The ECoG signals were gathered over a period of time at a sampling
rate of $1000\,$Hz. Finger flexion was recorded during the $2-$second movement
trials, where the computer displays randomly-chosen instructions $2$ seconds
per time and the patients performed self-paced movements $2-5$ times
throughout the trial. There was a $2-$second rest interval between two
instructions, during which the computer displayed complete darkness. For each
finger on each patient, the entire dataset captured finger flexion 30 times.

Additionally, five EEG subjects are selected from the GigaScience database
[27], where the motor imagery-based neural signals are collected in a similar
manner.

To train our models, we split the neural signals from the whole trial into
multiple segments. While each gesture recording lasted $2$ seconds, we have
split the signals in each two-second duration into two samples corresponding
to the same gesture (Figure 2). In total, we obtain $400$ samples for each
subject.

Figure 2: Each sample for training our model is obtained from one second of
the recording of the neural signals.

##  III Methods

We now describe the different steps in our approach for gesture
classification.

Spike encoding. We first encode the continuous time analog brain signal to a
discrete spike train with positive and negative spikes as per the standard
temporal contrast (TC) coding scheme [28]. As will be discussed later, we use
the sign of the spike to trigger positive and negative weight adaptation
during the learning phase.

The TD encoding algorithm is implemented as follows. Given a fixed threshold
$\theta_{th}$, for every time point starting from zero, the difference between
the value of the analog signal at the current time step $t_{k}$ and the
previous time step $t_{k-1}$ is calculated as $\Delta
u_{k}=f(t_{k})-f(t_{k-1})+u_{k-1}$. The absolute value of $\Delta u_{k}$ is
then compared with a threshold $\theta_{th}$. If $\Delta u$ is larger than the
$\theta_{th}$, a spike is issued at the current time instant, and the sign of
the spike is set to be the same as the sign of $\Delta u_{k}$. If $\Delta
u_{k}$ is less than the $\theta_{th}$, no spike will be issued. Also, when a
spike is issued, $\Delta u_{k}$ is reset to 0.

Spatial clustering of electrode channels. We next employ a spatial clustering
approach as a means to determine neuronal recordings that may be closely
correlated and to avoid combining measurements that are from vastly separated
regions. For the ECoG signals, the standard k-means clustering algorithm is
used to partition the electrodes into subgroups based on their location in the
euclidean space. Figure 3 illustrates $n_{c}=5$ clusters obtained for one of
the subjects in the ECoG dataset using the standard Elbow method.

For the data from the EEG dataset, the spatial clustering step using k-means
algorithm can be omitted as the clusters can be obtained automatically using
the spatial distribution information available for the 64 channels (Fig. 4).

Figure 3: Spatial cluster outcome as well as its 3D cortical projection
according to the ‘location’ variable for one of the subjects from the ECoG
dataset. Figure 4: EEG channel configuration, adapted from [29].

Convolutional spiking neural network. The spike signals from all the channels
in each cluster are then processed by a convolutional spiking neural network.
We illustrate our approach in Fig. 5 for an exemplary case of 4 channels in a
cluster. The spike values from the channels {-1,0,1} are unfolded along the
$1000$ timesteps of the temporal dimension for each sample and convolved with
a $3\times 3$ convolution kernel. We use a distinct kernel with the same
initial random values for every three channels in the spatial dimension with a
stride of 1, which means that each kernel will only stride across the temporal
dimension once and record the spike dynamics in three adjacent channels. The
output of the kernel feeds into a standard integrate-and-fire spiking neuron.

Figure 5: Illustration of the convolutional SNN implementation.

Bio-inspired unsupervised learning. In order to distinguish the spiking
dynamics from different samples, a modified event-driven weight adaptation
learning rule is employed to update the convolution kernel. The learning in
our model is triggered when the postsynaptic neuron spikes as a result of the
integration of the spikes through the convolutional filter during the
feedforward computation. The detailed description of the algorithm is as
follows.

At every convolution step $t$, if the output neuron issues a spike, we update
the weights corresponding to all the pixels of the kernel that have an
incident spike. The kernel weight is potentiated if the output spike is
positive according to the following rule depending on the sign of the input
spike:

$\Delta W=\begin{cases}\exp(\tau-t-\tau_{r})&\text{if input spike is
positive,}\\\ \exp(-(\tau-t-\tau_{r})^{2})&\text{if input spike is
negative.}\\\ \end{cases}$ (1)

Here, $\tau$ is the temporal location of the input spike relative to the
current time instant, $t-\tau=0,1,2$ as we use a $3\times 3$ kernel;
$\tau_{r}$ represents a fixed learning rate which is a hyperparameter for our
model. If the output spike is negative, the weights are depressed by the same
magnitude as above. If there are no spikes at the output neuron, the kernel
weights remain unchanged. The membrane potential of the output neuron does not
have a leak and is reset only when there is a spike.

KNN Classifier. After the input spikes for the full duration are presented and
the kernel weights updated in an event-triggered manner as described above,
each sample is represented by a concatenated kernel weight vector of dimension
$[3\times 3\times(n_{e}-2n_{c})]\times 1$. In the classification phase, these
kernel features are fed into a standard KNN classifier along with the labels
of the gesture. During the testing phase, the kernel feature vectors are
similarly calculated using the above-described weight update rule, which is
then fed to the KNN classifier to obtain a decision class.

##  IV Results

Datasets obtained from the Standford Human Cortical Activity [26] as well as
GigaScience database Library [27] are used for evaluating the proposed
algorithmic framework. We first report the effect of the number of clusters on
the classification accuracy (Fig. 6). In our study, the best classification
accuracies were obtained when the two subjects in the ECOG dataset were
partitioned into $n_{c}=5$ clusters and the five subjects from the EEG dataset
into $n_{c}=9$ clusters. The number of electrodes or channels varies between
$6$ and $24$ in the ECOG subjects, and between $3$ and $9$ in EEG subjects.

For the two subjects selected from the ECoG dataset, we report a
classification accuracy of $97.07\%$ (subject ‘bp’), and $95.85\%$ (subject
‘ca’). For the EEG data, with minor adjustments to the hyperparameters, our
methodology results in classification accuracy ranging between $92.74\%$ to
$96.51\%$ for $5$ different subjects in this dataset. The confusion matrix for
the 6 gestures is reported in Fig. 7. The values of hyperparameters used in
our simulations are presented in table I.

To benchmark our results, we report results from the literature that has
applied various signal processing techniques for similar BCI classification
tasks (Table II). To the best of our knowledge, the algorithmic framework we
have presented achieves superior accuracy compared to all the prior art.

Figure 6:  Classification accuracy as a function of the number of clusters to
which the electrodes are grouped. The best performance is obtained with 5
clusters for the ECOG data and 9 clusters for the EEG data. Figure 7: Combined
confusion matrix for the two subjects from the ECoG dataset (left) and 5
subjects from the EEG dataset (right). ‘ML’ is left hand movement, ‘MR’ is
right hand movement, ‘MIL’ is motor imagery left hand, ‘MIR’ is motor imagery
right hand. TABLE I: List of Hyperparameters used in our simulation.

Parameter ECoG EEG Encoding Threshold ($\theta_{th}$) 0.22 0.18 K-means
Cluster Number ($n_{c}$) 5 9 Kernel Dimension 3$\times$3 3$\times$3 Temporal
Dimension Stride 3 3 Channel Dimension Stride 1 1 Unsupervised Learning Rate
($\tau_{r}$) 5 6 Convolution Threshold 0.1 0.12 KNN’s K Value 5 5

TABLE II: Comparison of our work with state-of-art BCI classification.

Reference Task Probe Location Method Accuracy (%) or Cross-correlation(cc)
[17] Hand movement EEG BI-SNN 0.55-0.74 (cc) [18] Hand movement EEG CNN and
WNN 77.9%-90.1% [19] Hand movement EEG MSGMTL 73.10% [20] Hand movement EEG
expert CNN 73.78%-91.17% [21] Hand movement EEG OPTICAL+ 69.59% [22] Hand
movement EEG CNNs 68.19% [23] Finger Flexion Primary Motor Cortex ECoG
Temporal lightGBM 77.0% [24] Finger movement Primary Motor Cortex ECoG linear
SVM 60.6%-70.7% [25] Finger movement Primary Motor Cortex ECoG ResOT-$L_{2}$
61.3% Our Work Finger Flexion Primary Motor Cortex ECoG Convolutional SNN-
based Spatiotemporal Signal Analysis 97.07% and 95.85% for two subjects Our
Work Motor Imagery Hand Movement EEG Convolutional SNN- based Spatiotemporal
Signal Analysis 92.74%-96.51% for five subjects

##  V Conclusion

Time-varying signals obtained from BCI techniques such as ECOG and EEG encode
rich spatiotemporal patterns capturing motor imagery, though their analysis is
challenging due to subject-to-subject variations and inherent noise in
measurements. In this paper, we proposed an efficient algorithmic framework
for the classification of finger gestures based on two exemplary datasets,
that combine spatial clustering of the measurement channels and temporal
analysis of the spike-encoded signals based on a convolutional spiking neural
network. The kernel weights of the network evolve in an event-driven manner
based on the history of spikes over a short duration, which can then be used
to identify the gesture class with superior accuracy compared to methods in
the prior art. Hence, our method can be used for in situ and real-time
classification of motor imagery based on locally available spatial and
temporal information for BCI applications.

## Acknowledgment

This research was supported in part by an Open Fellowship of the EPSRC
(EP/X011356/1).

## References

  * [1] Eric C Leuthardt, Gerwin Schalk, Jonathan R Wolpaw, Jeffrey G Ojemann, and Daniel W Moran.  A brain–computer interface using electrocorticographic signals in humans.  Journal of neural engineering, 1(2):63, 2004.
  * [2] Luca Citi, Riccardo Poli, Caterina Cinel, and Francisco Sepulveda.  P300-based BCI mouse with genetically-optimized analogue control.  IEEE transactions on neural systems and rehabilitation engineering, 16(1):51–61, 2008.
  * [3] Sung-Phil Kim, John D Simeral, Leigh R Hochberg, John P Donoghue, and Michael J Black.  Neural control of computer cursor velocity by decoding motor cortical spiking activity in humans with tetraplegia.  Journal of neural engineering, 5(4):455, 2008.
  * [4] Eleanor A Curran and Maria J Stokes.  Learning to control brain activity: A review of the production and control of EEG components for driving brain–computer interface (BCI) systems.  Brain and cognition, 51(3):326–336, 2003.
  * [5] Dennis J McFarland and Jonathan R Wolpaw.  Brain-computer interface operation of robotic and prosthetic devices.  Computer, 41(10):52–56, 2008.
  * [6] Nanying Liang and Laurent Bougrain.  Decoding finger flexion from band-specific ECoG signals in humans.  Frontiers in neuroscience, 6:91, 2012.
  * [7] Rami Alazrai, Hisham Alwanni, and Mohammad I Daoud.  EEG-based BCI system for decoding finger movements within the same hand.  Neuroscience letters, 698:113–120, 2019.
  * [8] Jan Kubanek, Kai J Miller, Jeffrey G Ojemann, Jonathan R Wolpaw, and Gerwin Schalk.  Decoding flexion of individual fingers using electrocorticographic signals in humans.  Journal of neural engineering, 6(6):066001, 2009.
  * [9] Xiang Zhang, Lina Yao, Xianzhi Wang, Jessica Monaghan, David Mcalpine, and Yu Zhang.  A survey on deep learning based brain computer interface: Recent advances and new frontiers.  J Neural Eng., 66, 2021.
  * [10] Zhongke Gao, Tao Yuan, Xinjun Zhou, Chao Ma, Kai Ma, and Pan Hui.  A deep learning method for improving the classification accuracy of SSMVEP-based BCI.  IEEE Transactions on Circuits and Systems II: Express Briefs, 67(12):3447–3451, 2020.
  * [11] Michael Tangermann, Klaus-Robert Müller, Ad Aertsen, Niels Birbaumer, Christoph Braun, Clemens Brunner, Robert Leeb, Carsten Mehring, Kai Miller, Gernot Mueller-Putz, Guido Nolte, Gert Pfurtscheller, Hubert Preissl, Gerwin Schalk, Alois Schlögl, Carmen Vidaurre, Stephan Waldert, and Benjamin Blankertz.  Review of the BCI competition iv.  Frontiers in Neuroscience, 6, 2012.
  * [12] Yu Hu, Yongkang Wong, Wentao Wei, Yu Du, Mohan Kankanhalli, and Weidong Geng.  A novel attention-based hybrid CNN-RNN architecture for sEMG-based gesture recognition.  PloS one, 13(10):e0206049, 2018.
  * [13] Wolfgang Maass.  Networks of spiking neurons: the third generation of neural network models.  Neural networks, 10(9):1659–1671, 1997.
  * [14] Bleema Rosenfeld, Osvaldo Simeone, and Bipin Rajendran.  Spiking generative adversarial networks with a neural network discriminator: Local training, bayesian models, and continual meta-learning.  IEEE Transactions on Computers, 71(11):2778–2791, 2022.
  * [15] Jiechen Chen, Nicolas Skatchkovsky, and Osvaldo Simeone.  Neuromorphic wireless cognition: Event-driven semantic communications for remote inference.  IEEE Transactions on Cognitive Communications and Networking, 9(2):252–265, 2023.
  * [16] Bipin Rajendran, Abu Sebastian, Michael Schmuker, Narayan Srinivasa, and Evangelos Eleftheriou.  Low-power neuromorphic hardware for signal processing applications: A review of architectural and system-level design approaches.  IEEE Signal Processing Magazine, 36(6):97–110, 2019.
  * [17] Kaushalya Kumarasinghe, Nikola Kasabov, and Denise Taylor.  Brain-inspired spiking neural networks for decoding and understanding muscle activity and kinematics from electroencephalography signals during hand movements.  Scientific reports, 11(1):2486, 2021.
  * [18] Zhiwen Zhang, Feng Duan, Jordi Sole-Casals, Josep Dinares-Ferran, Andrzej Cichocki, Zhenglu Yang, and Zhe Sun.  A novel deep learning approach with data augmentation to classify motor imagery signals.  IEEE Access, 7:15945–15954, 2019.
  * [19] Xianxiong Zhang, Qingshan She, Tongcai Tan, Yunyuan Gao, Yuliang Ma, and Jianhai Zhang.  Multi-source geometric metric transfer learning for EEG classification.  Biomedical Signal Processing and Control, 81:104435, 2023.
  * [20] Brenda E Olivas-Padilla and Mario I Chacon-Murguia.  Classification of multiple motor imagery using deep convolutional neural networks and spatial filters.  Applied Soft Computing, 75:461–472, 2019.
  * [21] Shiu Kumar, Ronesh Sharma, and Alok Sharma.  OPTICAL+: a frequency-based deep learning scheme for recognizing brain wave signals.  Peerj Computer Science, 7:e375, 2021.
  * [22] James R Stieger, Stephen A Engel, Daniel Suma, and Bin He.  Benefits of deep learning classification of continuous noninvasive brain–computer interface control.  Journal of neural engineering, 18(4):046082, 2021.
  * [23] Lin Yao, Bingzhao Zhu, and Mahsa Shoaran.  Fast and accurate decoding of finger movements from ECoG through riemannian features and modern machine learning techniques.  Journal of Neural Engineering, 19(1):016037, 2022.
  * [24] Hyemin S Lee, Leonhard Schreiner, Seong-Hyeon Jo, Sebastian Sieghartsleitner, Michael Jordan, Harald Pretl, Christoph Guger, and Hyung-Soon Park.  Individual finger movement decoding using a novel ultra-high-density electroencephalography-based brain-computer interface system.  Frontiers in Neuroscience, 16, 2022.
  * [25] Bingzhao Zhu, Masoud Farivar, and Mahsa Shoaran.  Resot: Resource-efficient oblique trees for neural signal classification.  IEEE Transactions on Biomedical Circuits and Systems, 14(4):692–704, 2020.
  * [26] Kai J Miller, Dora Hermes, Christopher J Honey, Adam O Hebb, Nick F Ramsey, Robert T Knight, Jeffrey G Ojemann, and Eberhard E Fetz.  Human motor cortical activity is selectively phase-entrained on underlying rhythms.  PLOS Computational Biology, 8(9): e1002655, 2012.
  * [27] Hohyun Cho, Minkyu Ahn, Sangtae Ahn, Moonyoung Kwon, and Sung Chan Jun.  EEG datasets for motor imagery brain–computer interface.  GigaScience, 6(7), 2017.
  * [28] Sander M Bohte, Han La Poutré, and Joost N Kok.  Unsupervised clustering with spiking neurons by sparse temporal coding and multilayer RBF networks.  IEEE Transactions on neural networks, 13(2):426–435, 2002.
  * [29] Lu Qiu and Wenya Nan.  Brain network constancy and participant recognition: An integrated approach to big data and complex network analysis.  Frontiers in Psychology, 11:1003, 2020.