11institutetext: 1Lenovo

# 1st Place Solution for CVPR2023 BURST Long Tail and Open World Challenges

Kaer Huang1

###### Abstract

Currently, Video Instance Segmentation (VIS) aims at segmenting and
categorizing objects in videos from a closed set of training categories that
contain only a few dozen of categories, lacking the ability to handle diverse
objects in real-world videos. As TAO and BURST datasets release, we have the
opportunity to research VIS in long-tailed and open-world scenarios.
Traditional VIS methods are evaluated on benchmarks limited to a small number
of common classes, But practical applications require trackers that go beyond
these common classes, detecting and tracking rare and even never-before-seen
objects. Inspired by the latest MOT paper for the long tail task (Tracking
Every Thing in the Wild, Siyuan Li etl), we use the same idea for VIS long-
tailed and open-world setting where classification is the bottleneck. We add a
mask branch on the TETer-SwinL model for segmentation and tuned the
classification loss weight to 1/100 of the original setting to ignore more
effort in optimizing classification loss. and then using the cluster method
(Class Exemplar Matching) to identify classification labels. We use BUSRT mask
labels to generate ground box annotation for training and inference. for the
BURST long tail challenge, we train our model on a combination of LVISv0.5 and
the COCO dataset using repeat factor sampling. First, train the detector with
segmentation and CEM on LVISv0.5 + COCO dataset. And then, train the instance
appearance similarity head on the TAO dataset. at last, our method (LeTracker)
gets 14.9 HOTAall in the BURST test set, ranking 1st in the benchmark. for the
open-world challenges, we only use 64 classes (Intersection classes of BURST
Train subset and COCO dataset, without LVIS dataset) annotations data
training, and testing on BURST test set data and get 61.4 OWTAall, ranking 1st
in the benchmark. Our code will be released to facilitate future research.

###### Keywords:

VIS, Long Tail, Open World

##  1 Introduction

Video Instance Segmentation (VIS) aims to recognize, localize, segment, and
track objects in a given video sequence. It is a cornerstone of dynamic scene
analysis and vital for many real-world applications such as autonomous
driving, augmented reality, video surveillance, and short video analysis.
Despite impressive efforts, Video Instance Segmentation is fundamentally
constrained to segment and classify objects from a small closed set of
training categories, thus limiting the capacity to generalize to large scale
and diverse real-world[2]. Most VIS tasks are evaluated on benchmarks limited
to a small number of common classes. Practical applications require trackers
that go beyond these common classes, detecting and tracking rare and even
never-before-seen objects. As TAO and BURST datasets release, we have the
opportunity to push VIS research to long-tail or open-world tasks.

Figure 1: TETer

##  2 Related Work

Video Instance Segmentation (VIS) is tracking all instances belonging to a
fixed category set and obtaining their independent masks in every frame. As a
combination of identification, tracking, and segmentation, its output consists
of the corresponding predictions for different instances, each result
including a category label, a confidence score, and a binary mask sequence of
this instance. Most of the early research focuses on frame-based online
methods. MaskTrack R-CNN proposes an extra tracking branch adding to the
existing Mask R-CNN framework, setting a simple but strong baseline for VIS. A
series of subsequent works follow this tracking-by-detection paradigm, coupled
with explicit instance associations. For example, IDOL learns discriminative
instance representations through comparative learning, achieving good
performance results. Although with good real-time performance, these frame-
based methods lack the ability to capture long-term temporary dependency.
However, this can be realized naturally in another clip-based manner, which
divides frames of videos into several clips. Multi-frame modeling within the
clip is performed, followed by instance matching of the inter-clip. These
offline clip-based methods achieve better performance because they capture
richer spatiotemporal features. For example, VisTR uses a transformer inside
the clip for parallel end-to-end sequence prediction. Despite the performance
gains achieved by the above methods, they are less effective in scenarios with
congestion, occlusion, high-speed motion, and similar instances. To address
this issue, MDQE and InstMove propose to use spatiotemporal priors to
initialize object queries and to explicitly predict position and motion
deformations respectively, providing significant improvements in the
segmentation of challenging scenes.

Long Tail Tracking is a new research area that requires building trackers that
work for rare objects, that may only contain a few examples in the training
set. Contemporary MOT are designed for closed-set scenarios where all objects
appear frequently in the training and testing data distributions. Hence, Dave
et al proposed a new benchmark, TAO, that focuses on studying MOT in the long
tail of the object category distribution. On this benchmark, AOA, QDTrack, and
TET achieve impressive performance. However, those works are still limited to
pre-defined object categories and thus do not scale to the diversity of real-
world settings.

Open World VIS is a new research area that requires building trackers that can
generalize to never-before-seen objects. OWTB proposed Open World Tracking to
segment and track all the objects in videos. OWTB achieves state-of-the-art
performance on several multi-object tracking datasets. Recent work on UVO
focuses on class-agnostic object segmentation and tracking in videos,
neglecting the object classification performance during evaluation. BURST[14]
mainly follows the evaluation protocols of class-agnostic Multiple Object
Tracking, while only the 80 training categories in MS-COCO are measured for
category-wise evaluation. So methods adopted in BURST, such as STCN, have no
ability to classify objects from novel categories.

Multi Object Tracking (MOT) is a very general algorithm and has been studied
for many years. The mainstream methods follow the tracking-by-detection
paradigm [8, 10, 12, 15]. With the development of deep learning in recent
years, the performance of the detection model is improved rapidly. Currently,
most of the work relies on YOLOX [16, 18]. Our method selected a stronger
performance network CBNetV2 [9] which is used to verify the potential of the
detector in our hypothesis. Another important component of MOT is an
association strategy. Popular association methods include motion-based (IoU
matching, Kalman filter) [1], appearance-based (ReID embedding) [14],
transformer-based [17], or the combination of them [19, 11]. Our methods
remove all motion information and use only a high-performance appearance
model.

Multi Object Tracking and Segmentation (MOTS) is highly related to MOT by
changing the form of boxes to fine-grained mask representation [16]. Many MOTS
methods are developed upon MOT trackers [6, 13, 5]. Our ideas are similar to
theirs. A mask header was added on the basis of the MOT network in our MOTS
solution.

Self-Supervised Learning has made significant progress in representation
learning in recent years. Contrastive learning, one of the self-supervised
learning methods such as MoCo[4], SimCLR[2], BYOL[3], etc, has performance
that is getting closer to results of supervised learning methods in ImageNet
dataset. We leveraged Momentum Contrastive Learning (MoCo-v2)[4] to train a
new appearance embedding model without using tracking annotations. The
technique not only meets the requirements of SSMOT and SSMOTS[5] but also
improves the performance of the appearance model.

##  3 Method

###  3.1 Overall Architecture

Our method gets inspiration from the Long Tail MOT paper, which takes the
cluster method to learn the classification information of the target object.
As shown in Table 1, We add a mask branch to TETer for the VIS task. We also
tuned the classification loss weight to 1/100 of the original setting to
ignore more effort in optimizing classification loss. Then we train and infer
this model on the BURST dataset. Next, we will detail explain every part of
the model in the next section.

###  3.2 Detection and Segmentation

Our method has a general class-agnostic detection with segmentation, As shown
in Figure 1. Detail implementation based on FasterRCNN with SwinL as the
backbone. The difference with normal FasterRCNN detection is the
classification branch network which changes from class num (N) to binary
classes (Yes/No) which means whether it includes objects. Detailed
implementation is shown in Figure 2.

###  3.3 Clustering Class Exemplar Label

Like TETer[7], our method uses a clustering solution as a class learning
module which does not depend on the detection classification branch. As shown
in Figure 3.

###  3.4 TETer

As shown in Figure 4. The Main Architect of our solution is the same as TETer
except we adapt it to the VIS task.

##  4 Results

###  4.1 Dataset

BURST provides 2,914 videos with pixel-precise labels for 16,089 unique object
tracks (600,000 per-frame masks) spanning 482 object classes. For the Long
Tail challenge which allows the use of the LVIS dataset. So we train detection
and CEM module parameters by merging COCO trainset with LVISv0.5, and then,
train the instance appearance similarity head on the BURST train dataset. For
Open World Challenge, LVIS data is not allowed on this challenge. So we just
only Use classes (64) that intersect COCO and BURST trainset.

###  4.2 Metrics

We evaluate all tasks using Higher Order Tracking Accuracy (HOTA) because it
strikes a good balance between measuring frame-level detection and temporal
association accuracy. For the open-world task, a slightly modified, recall-
based variant of HOTA called Open World Tracking Accuracy (OWTA) is used.

###  4.3 Main Results

Our detailed results are shown in Table 2 and Table 3. Compare with other
solutions, we get both champions in Long Tail and Open World Challenges. from
the analysis comparison results between common and uncommon classes, our
solution has a big improvement with uncommon class tracking.

Table 1: Long Tail Result Tracker | HOTAall | DETAall | AssAall | HOTAcom | DETAcom | AssAcom | HOTAunc | DETAunc | AssAunc
---|---|---|---|---|---|---|---|---|---
STCN Tracker | 4.5 | 5.4 | 4.6 | 17.1 | 19.6 | 16.7 | 2.0 | 2.6 | 2.2
Box Tracker | 5.7 | 5.4 | 6.8 | 20.1 | 19.6 | 23.2 | 2.9 | 2.6 | 3.6
LeTracker(Our) | 14.9 | 12.3 | 20.1 | 36.1 | 30.7 | 47.8 | 10.8 | 8.7 | 14.6
Table 2: Open World Result Tracker | OWTAall | DETReall | AssAall | OWTAcom | DETRecom | AssAcom | OWTAunc | DETReunc | AssAunc
---|---|---|---|---|---|---|---|---|---
STCN Tracker2 | 57.5 | 61.6 | 54.1 | 62.9 | 71.5 | 55.7 | 23.9 | 21.0 | 28.6
OWTB | 56.4 | 70.7 | 45.5 | 59.9 | 76.6 | 47.4 | 38.3 | 45.7 | 33.6
Box Tracker | 55.9 | 61.5 | 51.1 | 61.0 | 71.4 | 52.5 | 24.6 | 21.1 | 30.0
LeTracker(Our) | 61.4 | 74 | 51.5 | 64.3 | 78.7 | 52.9 | 47.1 | 54.6 | 42.0

##  5 Conclusition

Traditional VIS methods are evaluated on benchmarks limited to a small number
of common classes and do not have long tail issues, But practical applications
require trackers that go beyond these common classes, detecting and tracking
rare and even never-before-seen objects. Inspired by the latest MOT paper for
the long tail task (Tracking Every Thing in the Wild, Siyuan Li et al), we use
the same idea for VIS long-tailed and open-world setting where classification
is the bottleneck. We add a mask branch on the TETer-SwinL model for
segmentation and tuned the classification loss weight to 1/100 of the original
setting to ignore more effort in optimizing classification loss. and then
using the cluster method (Class Exemplar Matching) to identify classification
labels. We use BUSRT mask labels to generate ground box annotation for
training and inference. for the BURST long tail challenge, we train our model
on a combination of LVISv0.5 and the COCO dataset using repeat factor
sampling. First, train the detector with segmentation and CEM on LVISv0.5 +
COCO dataset. And then, train the instance appearance similarity head on the
TAO train dataset. at last, our method (LeTracker) gets 14.9 HOTAall in the
BURST test set, ranking 1st in the benchmark. for the open-world challenges,
we only use 64 classes (Intersection classes of BURST Train subset and COCO
dataset, without LVIS dataset) annotations data training and testing on BURST
test set data and get 61.4 OWTAall, ranking 1st in the benchmark.

## References

  * [1] Bewley, A., Ge, Z., Ott, L., Ramos, F., Upcroft, B.: Simple online and realtime tracking. In: 2016 IEEE international conference on image processing (ICIP). pp. 3464–3468. IEEE (2016)
  * [2] Chen, T., Kornblith, S., Norouzi, M., Hinton, G.: A simple framework for contrastive learning of visual representations. In: International conference on machine learning. pp. 1597–1607. PMLR (2020)
  * [3] Grill, J.B., Strub, F., Altché, F., Tallec, C., Richemond, P., Buchatskaya, E., Doersch, C., Avila Pires, B., Guo, Z., Gheshlaghi Azar, M., et al.: Bootstrap your own latent-a new approach to self-supervised learning. Advances in neural information processing systems 33, 21271–21284 (2020)
  * [4] He, K., Fan, H., Wu, Y., Xie, S., Girshick, R.: Momentum contrast for unsupervised visual representation learning. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 9729–9738 (2020)
  * [5] Huang, K., Lertniphonphan, K., Chen, F., Li, J., Wang, Z.: Multi-object tracking by self-supervised learning appearance model. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 3162–3168 (2023)
  * [6] Ke, L., Li, X., Danelljan, M., Tai, Y.W., Tang, C.K., Yu, F.: Prototypical cross-attention networks for multiple object tracking and segmentation. Advances in Neural Information Processing Systems 34, 1192–1203 (2021)
  * [7] Li, S., Danelljan, M., Ding, H., Huang, T.E., Yu, F.: Tracking every thing in the wild. In: Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXII. pp. 498–515. Springer (2022)
  * [8] Liang, C., Zhang, Z., Zhou, X., Li, B., Zhu, S., Hu, W.: Rethinking the competition between detection and reid in multiobject tracking. IEEE Transactions on Image Processing 31, 3182–3196 (2022)
  * [9] Liang, T., Chu, X., Liu, Y., Wang, Y., Tang, Z., Chu, W., Chen, J., Ling, H.: Cbnetv2: A composite backbone network architecture for object detection. arXiv preprint arXiv:2107.00420 (2021)
  * [10] Lu, Z., Rathod, V., Votel, R., Huang, J.: Retinatrack: Online single stage joint detection and tracking. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 14668–14678 (2020)
  * [11] Pang, J., Qiu, L., Li, X., Chen, H., Li, Q., Darrell, T., Yu, F.: Quasi-dense similarity learning for multiple object tracking. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 164–173 (2021)
  * [12] Sun, P., Cao, J., Jiang, Y., Zhang, R., Xie, E., Yuan, Z., Wang, C., Luo, P.: Transtrack: Multiple object tracking with transformer. arXiv preprint arXiv:2012.15460 (2020)
  * [13] Voigtlaender, P., Krause, M., Osep, A., Luiten, J., Sekar, B.B.G., Geiger, A., Leibe, B.: Mots: Multi-object tracking and segmentation. In: Proceedings of the ieee/cvf conference on computer vision and pattern recognition. pp. 7942–7951 (2019)
  * [14] Wang, Z., Zhao, H., Li, Y.L., Wang, S., Torr, P., Bertinetto, L.: Do different tracking tasks require different appearance models? Advances in Neural Information Processing Systems 34, 726–738 (2021)
  * [15] Wu, J., Cao, J., Song, L., Wang, Y., Yang, M., Yuan, J.: Track to detect and segment: An online multi-object tracker. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 12352–12361 (2021)
  * [16] Yan, B., Jiang, Y., Sun, P., Wang, D., Yuan, Z., Luo, P., Lu, H.: Towards grand unification of object tracking. arXiv preprint arXiv:2207.07078 (2022)
  * [17] Zeng, F., Dong, B., Wang, T., Zhang, X., Wei, Y.: Motr: End-to-end multiple-object tracking with transformer. arXiv preprint arXiv:2105.03247 (2021)
  * [18] Zhang, Y., Sun, P., Jiang, Y., Yu, D., Yuan, Z., Luo, P., Liu, W., Wang, X.: Bytetrack: Multi-object tracking by associating every detection box. arXiv preprint arXiv:2110.06864 (2021)
  * [19] Zhang, Y., Wang, C., Wang, X., Zeng, W., Liu, W.: Fairmot: On the fairness of detection and re-identification in multiple object tracking. International Journal of Computer Vision 129(11), 3069–3087 (2021)


# The jump effect of a general eccentric cylinder rolling on a ramp

E. Aldo Arroyo$~\\!{}^{(a)}$ and M. Aparicio Alcalde$~\\!{}^{(b)}$,
(a)Centro de Ciências Naturais e Humanas, Universidade Federal do ABC,
Santo André, 09210-170 São Paulo, SP, Brazil.
(b)Instituto de Ciências Exatas e Tecnológicas, Universidade Federal de
Viçosa,
38810-000, Rio Paranaíba, MG, Brazil aldo.arroyo<EMAIL_ADDRESS>
###### Abstract

Interesting phenomena occur when an eccentric rigid body rolls on an inclined
or horizontal plane. For example, a variety of motions between rolling and
sliding are exhibited until suddenly a jump occurs. We provide a detailed
theoretical description of the jump effect for a general eccentric cylinder.
Before the jump, when the cylinder moves along the ramp, we can assume a pure
rolling motion. However, it turns out that when the cylinder reaches its
jumping position, both the normal and static frictional forces approach zero.
Thus, it seems that there will no longer be sufficient force to maintain
rolling without slip. In order to have a jump without slipping, we prove that
the parameters that characterize the dynamic behavior of the cylinder must
belong to some restricted region.

##  1 Introduction

The physical system consists of a rigid body, such as hoops, wheels, disks,
and spheres, whose center of mass is located at a distance $d\neq 0$ from the
geometric center and rolls on a horizontal or inclined plane with friction.
This system has interesting and unexpected dynamic behavior, which has
attracted the attention of the community. In Fig. 1, a general configuration
of the physical system is exhibited, where we can identify an inclined ramp
making an angle $\alpha$ with the horizontal, and a cylinder of radius $R$
with its center of mass located at a distance $d$ from the geometric center.
The cylinder rolls over the ramp, and its motion is tracked by the angle
$\theta$.

Figure 1: Schematic configuration of the physical system showing a cylinder
rolling down a ramp of angle $\alpha$. The position of the geometric center
$C$, and the center of mass $CM$ of the cylinder with respect to point $O$ are
given by the vectors $\vec{l}+\vec{b}$ and $\vec{l}+\vec{b}+\vec{d}$,
respectively. The vector $\vec{F}$ represents the force acting at the contact
point between the cylinder and the ramp, and $M\vec{g}$ is the cylinder’s
weight.

In the case of a horizontal plane (i.e., $\alpha=0^{\circ}$), there is
literature on the dynamics of a particular cylinder. This cylinder consists of
a very thin, massless cylindrical shell with a point mass stuck to its
surface. The motion of this cylinder has been studied and debated by several
authors [1, 2]. When the point mass attached to the cylinder is initially at
the highest point of the cylinder, and the system is released from rest (i.e.,
the initial conditions are such that $\theta_{0}=0$ and $\dot{\theta}_{0}=0$),
Tokieda [2] showed that the cylinder must jump at $\theta=90^{\circ}$.
However, subsequent articles have shown that the assumption made by the
author, that the jump happens immediately after the pure rolling motion, is
incorrect. For instance, it has been shown that before the jump, there must be
sliding [3, 4, 5]. Despite this fact, Pritchett [4] obtained numerical and
experimental results (for a hula hoop with a stuck point mass) showing that
the jump happens around $\theta=90^{\circ}$, even when considering slippage.
Moreover, to test the jump at this angle of $90^{\circ}$, Theron [6] claimed
that the elasticity of the hula hoop must be considered. As observed, despite
the physical system being simple, the motion of the cylinder turns out to be,
in general, complex. Further analysis has shown that a variety of these
motions can include self-induced jumping motions [7, 8, 9, 10, 11], as well as
multiple transitions back and forth from rolling to slipping [12, 13, 14, 15].

In the case of an inclined plane (i.e., $\alpha\neq 0$), we have found a few
models that consider this case. One of them studied an eccentric wheel [6],
and the other studied an eccentric disk [16]. These works hypothesize that the
wheel (or disk) must slip before jumping. Other references [17, 18, 19, 20]
have studied the equations of motion for this cylinder by means of the
Lagrangian formulation or a cumbersome torque analysis. However, these works
did not determine the position where the cylinder jumps.

In all the references we have found (with the exception of [21]), it was
assumed that the necessary condition for the jump is when the normal force
$F_{y}$ (acting perpendicularly from the ramp over the cylinder) vanishes at
the instant of the hop. Let us remark that a different jump condition will be
used in this work, and this condition turns out to be equivalent to that of
[21]. As argued in reference [21], the equation that the angle $\theta$ and
the angular velocity $\dot{\theta}$ must satisfy at the position where the
cylinder jumps are given by

$\displaystyle-g\cos\alpha+d\,\dot{\theta}^{2}\cos\theta=0.$ (1)

Assuming pure rolling motion, which implies conservation of energy, and using
the initial conditions $\theta_{0}=0$ and $\dot{\theta}_{0}=0$, it was
possible to find an angle $\theta_{J}$ that is a solution of Eq. (1).
Consequently, the position $l_{J}=R\theta_{J}$ along the inclined plane where
the cylinder jumps can be determined. However, there is a subtlety in the pure
rolling assumption. For a given value of the coefficient of static friction
$\mu_{s}$ between the cylinder and inclined plane, before the cylinder reaches
the position determined by the angle $\theta_{J}$, it may happen that the
static friction exceeds its maximum value. This would imply that the cylinder
slips before jumping.

Therefore, aside from the jump condition issue, we remark that there is an
important open question: Are there cases where an eccentric rigid body rolling
on an inclined (or horizontal) plane jumps just after pure rolling motion? As
we are going to see, the answer to this question is not trivial and depends on
the values of $\mu_{s}$, the angle $\alpha$, and the moment of inertia of the
rigid body considered. A general body that encompasses all types of eccentric
bodies having a cylindrical shape will be referred to as a general eccentric
cylinder. It turns out that the dynamic characteristics of this cylinder can
be parameterized using two parameters, $\chi$ and $k_{m}$. By exploring this
generality that we consider in our work, we prove that there are (quite
common) situations in which the jump occurs immediately after sliding roll,
and there are (less common) situations where the jump occurs immediately after
pure rolling. We also prove that when the jump occurs immediately after
sliding roll, we have $F_{y}=0$, and in the cases where the jump occurs
immediately after pure rolling, we have $F_{y}>0$. Therefore, it would be
correct to consider $F_{y}=0$ as a necessary condition for the jump whenever
there is sliding rolling immediately before the jump, as occurs in many
studies in the literature, but it would be incorrect to state the same when
the jump occurs immediately after pure rolling or in general cases.

This paper is organized as follows. In section 2, we will describe the
physical system to be analyzed and derive the corresponding equations of
motion. Then we will propose a model for the general eccentric cylinder. In
section 3, we determine the region where the parameters that characterize the
dynamical behavior of the general eccentric cylinder must belong in order to
have a jump without slip, for any value of $\alpha$ and $\mu_{s}$. In section
4, we will provide a summary and suggest further directions for exploration.

##  2 Description of the system and derivation of the equations of motion

In this section, we will define the physical system to be studied and derive
the corresponding equations of motion. The mechanical system consists of a
general eccentric cylinder rolling down an inclined ramp, as shown in Fig. 1.
The term eccentric means that the center of mass is located at a distance of
$d$ from the geometric center. By general, we mean that the mass distribution
within the cylinder is arbitrary, with the only restriction being that this
distribution is invariant under translations along the $z$-axis. Throughout
the motion, we will assume that the principal axis of the cylinder passing
through point $C$ always remains parallel to the $z$-axis. Note that when the
cylinder’s height is very small compared to the radius, the cylinder
represents a disk, and if the mass distribution is concentrated on the edge,
the cylinder becomes a hoop or wheel. In this sense, we can say that this
general eccentric cylinder can encompass all types of eccentric bodies that
have cylindrical shape.

In Fig. 1, we also show the coordinate system $xy$, where the origin is set at
the point $O$ located at the top of the ramp. The $x$-axis and $y$-axis are
parallel and perpendicular to the ramp, respectively. The position of the
center of mass with respect to the geometric center of the cylinder is given
by the vector

$\displaystyle\vec{d}=d\sin\theta\,\hat{i}+d\cos\theta\,\hat{j}\,,$ (2)

where $\theta$ is the angle between the vector $\vec{d}$ and the $y$-axis.

Using the coordinate system $xy$, which is an inertial reference frame, we can
write the position of the center of mass as follows

$\displaystyle\vec{r}_{CM}=\vec{l}+\vec{b}+\vec{d}=(l+d\sin\theta)\hat{i}+(b+d\cos\theta)\hat{j}\,.$
(3)

In the case where the cylinder remains in contact with the ramp, the vector
$\vec{b}$ is given by $\vec{b}=R\,\hat{j}$, where $R$ is the radius of the
cylinder. While in the case where the cylinder loses contact with the ramp and
flies, we have $|\vec{b}|>R$.

The two equations that are used to determine the dynamics of this rigid body
are:

  1. (i)

Newton’s second law for the motion of the center of mass:

$\displaystyle\vec{F}_{T}=M\ddot{\vec{r}}_{CM}\,,$ (4)

where the subscript $T$ means that we are considering the total force acting
on the cylinder; and

  2. (ii)

Newton’s second law for rotations:

$\displaystyle\vec{\tau}_{T}=I_{CM}\,\ddot{\vec{\theta}}\,,$ (5)

where $\ddot{\vec{\theta}}$ is defined as
$\ddot{\vec{\theta}}=-\ddot{\theta}\,\hat{k}$, and the total torque
$\vec{\tau}_{T}$ and moment of inertia $I_{CM}$ are taken around the center of
mass.

From the configuration of the system shown in Fig. 1, we can write Eq. (5) as
follows:

$\displaystyle-(\vec{b}+\vec{d})\times\vec{F}=-I_{CM}\ddot{\theta}\,\hat{k}\,.$
(6)

Since the total force acting on the cylinder is given by
$\vec{F}_{T}=\vec{F}+M\vec{g}$, using Eqs. (3) and (4), we obtain

$\displaystyle\vec{F}=M(\ddot{\vec{l}}+\ddot{\vec{b}}+\ddot{\vec{d}}-\vec{g})\,.$
(7)

Let us write the force $\vec{F}$ in terms of its components in the directions
of the $x$ and $y$-axis

$\displaystyle\vec{F}=F_{x}\hat{i}+F_{y}\hat{j}\,,$ (8)

note that these components $F_{x}$ and $F_{y}$ are the friction and normal
force, respectively. Now substituting Eqs. (3) and (8) into Eq. (7), we get

$\displaystyle F_{x}=M\ddot{l}+Md(\sin\theta)\ddot{~{}}-Mg\sin\alpha\,,$
$\displaystyle F_{y}=M\ddot{b}+Md(\cos\theta)\ddot{~{}}+Mg\cos\alpha\,.$ (9)

Using Eqs. (2), (8) and writing $\vec{b}=b\hat{j}$, from Eq. (6) we obtain

$\displaystyle I_{CM}\ddot{\theta}=d\sin\theta\,F_{y}-(b+d\cos\theta)F_{x}\,.$
(10)

The equations (9) and (10) will be useful in analyzing three possible types of
motion for the cylinder: pure rolling, rolling with slipping, and flight
motion. For each type of motion, additional relations between the forces and
positions must be established. In the next two subsections, we will study the
equations in the case of pure rolling and flight motion.

###  2.1 Equations in the case of pure rolling motion

Since the cylinder is in contact with the ramp, we have $b=R$, which means
that $b$ is constant and, therefore, $\ddot{b}=0$. Moreover, due to the pure
rolling condition, we also have $l=l_{0}+R(\theta-\theta_{0})$, and
consequently $\ddot{l}=R\,\ddot{\theta}$. The moments of inertia $I_{C}$ and
$I_{CM}$ with respect to the geometric center $C$ and the center of mass $CM$
of the cylinder are related by the equation $I_{C}=I_{CM}+Md^{2}$.
Substituting these equations into Eqs. (9) and (10), we can derive the
following nonlinear second-order differential equation

$\displaystyle\ddot{\theta}\left(\frac{I_{C}}{MR^{2}}+1+2\chi\cos\theta\right)-\chi\dot{\theta}^{2}\sin\theta-\frac{g}{R}(\chi\sin(\alpha+\theta)+\sin\alpha)=0\,,$
(11)

where we have defined the following dimensionless parameter $\chi$ as follows

$\displaystyle\chi=\frac{d}{R}\,.$ (12)

###  2.2 Equations in the case of flight motion

In the case of flight motion, the cylinder has no contact with the ramp,
therefore the contact force is null, i.e., $F_{x}=F_{y}=0$. From Eq. (10) it
is straightforward to show that the angular velocity $\dot{\theta}$ is
constant. And from Eq. (9) we obtain:

$\displaystyle\ddot{l}=d\,\dot{\theta}^{2}\sin\theta+g\sin\alpha\,,$
$\displaystyle\ddot{b}=d\,\dot{\theta}^{2}\cos\theta-g\cos\alpha\,.$ (13)

These equations are consistent with the common knowledge about the free-fall
motion of a rigid body, where the center of mass performs a parabolic
trajectory and the angular velocity is constant.

###  2.3 Transition from rolling to flight motion

In order to understand the condition for the transition from rolling to flight
motion, we will perform the following analysis. Before the jump, the cylinder
stays in contact with the ramp, so $b(t)=R$, $\dot{b}(t)=\ddot{b}(t)=0$, and
the values of $\theta(t)$, $\dot{\theta}(t)$, and $\ddot{\theta}(t)$ are
related by means of Eq. (10). At the moment of the jump, the values of $b(t)$,
$\theta(t)$, $\dot{b}(t)$, and $\dot{\theta}(t)$ change continuously. The
continuity of $\dot{b}(t)$ and $\dot{\theta}(t)$ is supported by the fact that
there are no additional external forces acting on the cylinder at the instant
of the jump that would change the cylinder’s momentum.

After the jump, the cylinder performs a flight motion where the relations in
Eq. (13) are valid. Due to the continuities mentioned above, the values of
$b=R$ and $\dot{b}=0$ are the initial conditions for $b(t)$ in the flight
motion. Therefore, when the value of $\ddot{b}$ given by Eq. (13) starts to be
greater than zero, it implies that the values of $\dot{b}$ and $b$ start to
increase. It is interesting to note that this increase in the values of
$\dot{b}$ and $b$ happens because we assumed that the cylinder is not attached
to the ramp. Consequently, at the point where the transition from rolling to
flight motion occurs, we must have $\ddot{b}=0$ in Eq. (13). This condition
means 111This jump condition was obtained in Ref. [21] using an alternative
approach.:

$\displaystyle d\,\dot{\theta}^{2}\cos\theta-g\cos\alpha=0\,.$ (14)

Before the jump, in order to understand the behavior of $\ddot{b}$ as defined
in Eq. (13), we need to track the value of $\dot{\theta}$. The motion of the
cylinder starts with $\dot{\theta}=0$, and due to the action of gravity,
$\dot{\theta}$ increases. Therefore, at the beginning,
$\ddot{b}=-g\cos\alpha<0$, and subsequently, $\ddot{b}$ increases because of
the contribution of $\dot{\theta}^{2}$ (with some oscillation due to the
factor $\cos\theta$). This implies that at some future instant, when
$\ddot{b}$ approaches zero, i.e., when Eq. (14) is satisfied, we reach the
moment when the cylinder jumps.

Let us comment that in other works [5, 6, 16], the jump condition has been
given by $F_{y}=0$. Namely, the jump happens at the point where the normal
force vanishes. However, these works consider the hypothesis that the jump is
not possible from pure rolling, which means that the cylinder must slide
before jumping. Although the main focus of our work is not on the slip case,
using equations (9) and (10), together with the relation
$F_{x}=\sigma\mu_{k}F_{y}$, where $\mu_{k}$ is the coefficient of kinetic
friction, and $\sigma=-1$ in case $\dot{l}>R\dot{\theta}$ (i.e., skidding
motion) and $\sigma=+1$ in case $\dot{l}<R\dot{\theta}$ (i.e., spinning
motion), we can show that the normal force $F_{y}$ is given by

$\displaystyle
F_{y}=\frac{M\,I_{CM}\left(d\,\dot{\theta}^{2}\cos\theta-g\cos\alpha\right)}{Md\left(\sigma\,\mu_{k}\sin\theta(d\cos\theta+R)-d\sin^{2}\theta\right)-I_{CM}}\,.$
(15)

From Eq. (15), we observe that the jump condition, as given by Eq. (14),
clearly implies that the normal force $F_{y}$ vanishes. It should be noted
that the expression for the normal force, as given by Eq. (15), is only valid
when there is slippage. If the jump occurs from pure rolling motion, employing
Eq. (14), we will show that the normal force does not necessarily vanish.

###  2.4 Scale invariance of the dynamics and a model for the general
eccentric cylinder

Considering a general mass distribution inside the cylinder, with the only
restriction being that this distribution is invariant under translations along
the cylinder’s principal axis, such that the center of mass does not coincide
with the cylinder’s geometric center, in this subsection we analyze the
equations of motion of this general eccentric cylinder. It turns out that the
equations of motion have scale invariance. By using this invariance, we can
find common characteristics between two different cylinders in a way that the
dynamics of both cylinders are equivalent. These common characteristics
between two different cylinders can be parameterized using two independent
parameters. We also propose a simplified model of the eccentric cylinder where
the two independent parameters have a simple geometric and mass distribution
interpretation. We prove that this simplified model is dynamically equivalent
to any general cylinder considered in our study.

Let us consider a cartesian coordinate system fixed to the cylinder, such that
the $z$-axis coincides with the principal axis, and the coordinate origin is
located at the cylinder’s geometric center $C$. The center of mass $CM$ and
the cylinder’s moment of inertia with respect to the $z$-axis are computed as
$\vec{r}_{CM}=\frac{1}{M}\int dxdydz\,\vec{r}\rho(x,y)$ and $I_{C}=\int
dxdydz\,(x^{2}+y^{2})\rho(x,y)$, respectively, where $\rho(x,y)$ is the mass
density bounded by the cylindrical surface. Since the mass distribution is
invariant under translations along the $z$-axis, the density $\rho(x,y)$ does
not depend on $z$. These specifications provide two general implications:

  1. 1.

The $CM$ can be located at any distance $d$ from the cylinder’s geometric
center $C$, with this distance restricted to the interval $d\in[0,R]$. For
example, the particular case where $d=R$ occurs when the entire cylinder’s
mass is distributed along a line on the cylindrical lateral surface that is
parallel to the $z$-axis.

  2. 2.

For fixed values of $M$, $R$, and $d$, the moment of inertia $I_{CM}$ (or
equivalently $I_{C}$, thanks to the relation $I_{C}=I_{CM}+Md^{2}$), has
minimum and maximum values, where these bounding values are:

    1. (a)

The minimum: this value corresponds to the case where the whole cylinder’s
mass is distributed along a line (parallel to the $z$-axis) that passes
through the $CM$, therefore in this case we have: $I_{CM}=0$ or
$I_{C}=Md^{2}$.

    2. (b)

The maximum: this value corresponds to the case where the whole cylinder’s
mass is distributed on the cylinder’s lateral surface, so in this case we
have: $I_{CM}=MR^{2}-Md^{2}$ or $I_{C}=MR^{2}$.

Given a value of $I_{C}$ such that $Md^{2}\leq I_{C}\leq MR^{2}$, there are a
variety of possibilities for the mass distribution $\rho(x,y)$ that yield the
same value of $I_{C}$. Since for a fixed values of $M$, $R$, $d$ and $g$, the
equations of motion, given by Eqs. (9) and (10), depend only on $I_{C}$, the
cylinder’s dynamical behavior does not depend on specific details of the mass
distribution $\rho(x,y)$.

From Eqs. (9) and (10), we can write the following equation

$\displaystyle\tilde{I}_{CM}\frac{d^{2}\tilde{\theta}}{d\tilde{t}^{2}}=\chi\sin\tilde{\theta}\,\left(\frac{d^{2}\tilde{b}}{d\tilde{t}^{2}}+\chi\frac{d^{2}\cos\tilde{\theta}}{d\tilde{t}^{2}}+\cos\alpha\right)-(\tilde{b}+\chi\cos\tilde{\theta})\,\left(\frac{d^{2}\tilde{l}}{d\tilde{t}^{2}}+\chi\frac{d^{2}\sin\tilde{\theta}}{d\tilde{t}^{2}}-\sin\alpha\right)\,,$
(16)

where we have defined the adimensional quantities
$\tilde{I}_{CM}=\frac{I_{CM}}{MR^{2}}$, $\tilde{t}=\sqrt{\frac{g}{R}}t$,
$\tilde{\theta}=\theta$, $\tilde{b}=b/R$, $\tilde{l}=l/R$ and $\chi$ is given
by Eq. (12). The functions $\tilde{\theta}$, $\tilde{b}$, and $\tilde{l}$ are
related to the solutions $\theta$, $b$, and $l$ of the equations of motion.
More explicitly, these relations are given by
$\theta(t)=\tilde{\theta}(\sqrt{\frac{g}{R}}t)$,
$b(t)=R\tilde{b}(\sqrt{\frac{g}{R}}t)$, and
$l(t)=R\tilde{l}(\sqrt{\frac{g}{R}}t)$. At this point, we are ready to analyze
the dependence of the cylinder’s dynamics with respect to $M$ and $R$.
Consider two cylinders, one with mass and radius given by the set $(M,R)$ and
the other by $(M^{\prime},R^{\prime})$. If $\chi$ and $\tilde{I}_{CM}$ are the
same for both cylinders, then due to Eq. (16), they share the same solutions
$\tilde{\theta}$, $\tilde{b}$, and $\tilde{l}$. Therefore, the set of
functions $(\theta^{\prime},b^{\prime},l^{\prime})$ and $(\theta,b,l)$ are
related by $\theta^{\prime}=\theta$, $b^{\prime}=\frac{R^{\prime}}{R}b$,
$l^{\prime}=\frac{R^{\prime}}{R}l$. Namely, $\theta$ is equal for both sets,
while $b$ and $l$ do not depend on $M$ and scale by a factor proportional to
$R$. Similarly, since $t=\sqrt{\frac{R}{g}}\tilde{t}$, the time of the events
that happen throughout the motion (like the time when the jump happens) scales
by a factor proportional to the square root of $R$.

From the above observations, and defining the scale factor
$\lambda=\frac{R^{\prime}}{R}$, we conclude that the equations of motion are
invariant under the scale transformation $\theta^{\prime}=\theta$,
$b^{\prime}=\lambda\,b$, $l^{\prime}=\lambda\,l$,
$t^{\prime}=\sqrt{\lambda}\,t$ provided that $\chi$ and $\tilde{I}_{CM}$ are
the same for the two different sets of values $(M,R)$ and
$(M^{\prime},R^{\prime})$. Since $\chi=\chi^{\prime}$, we have that
$d^{\prime}=\lambda\,d$. Furthermore, from
$\tilde{I}_{CM}=\tilde{I}^{\prime}_{CM}$, namely
$\frac{I_{CM}}{MR^{2}}=\frac{I^{\prime}_{CM}}{M^{\prime}R^{\prime 2}}$, or
$\frac{I_{C}}{MR^{2}}=\frac{I^{\prime}_{C}}{M^{\prime}R^{\prime 2}}$, and
using $I_{C}=\int dxdydz\,(x^{2}+y^{2})\rho(x,y)$, we get the relation

$\displaystyle\int
dxdydz\,\left(\left(\frac{x}{R}\right)^{2}+\left(\frac{y}{R}\right)^{2}\right)\frac{\rho(x,y)}{M}=\int
dx^{\prime}dy^{\prime}dz^{\prime}\,\left(\left(\frac{x^{\prime}}{R^{\prime}}\right)^{2}+\left(\frac{y^{\prime}}{R^{\prime}}\right)^{2}\right)\frac{\rho^{\prime}(x^{\prime},y^{\prime})}{M^{\prime}}\,,$

which means that under the transformation
$(x^{\prime},y^{\prime},z^{\prime})=(\lambda x,\lambda y,\lambda z)$, and the
change of mass from $M$ to $M^{\prime}$, we have that
$\frac{\rho(x,y)}{M}dxdydz=\frac{\rho^{\prime}(x^{\prime},y^{\prime})}{M^{\prime}}dx^{\prime}dy^{\prime}dz^{\prime}$,
or $\frac{dm}{M}=\frac{dm^{\prime}}{M^{\prime}}$, that is the mass densities
(divided by the total mass) of the cylinders are related by
$\frac{\rho(x,y)}{M}=\lambda^{3}\,\frac{\rho^{\prime}(x^{\prime},y^{\prime})}{M^{\prime}}$.

A summary of the previous analysis is as follows: given two different
cylinders, where to calculate their centers of mass and moments of inertia, we
use cartesian coordinate systems fixed to them. Let
$(x^{\prime},y^{\prime},z^{\prime})$ the coordinate system fixed to one
cylinder, and $(x,y,z)$ to the other one. If these coordinate systems are
related by the scale transformation
$(x^{\prime},y^{\prime},z^{\prime})=(\lambda x,\lambda y,\lambda z)$, and
given the mass densities $\rho(x,y)$ and
$\rho^{\prime}(x^{\prime},y^{\prime})$ such that
$\frac{\rho(x,y)}{M}=\lambda^{3}\,\frac{\rho^{\prime}(x^{\prime},y^{\prime})}{M^{\prime}}$,
this condition guarantees that $\chi=\chi^{\prime}$ and
$\tilde{I}_{CM}=\tilde{I}^{\prime}_{CM}$, then the dynamic behavior of the
cylinders are equivalent. Since the cylinder’s dynamics is basically ruled by
the two parameters $\chi$ and $\tilde{I}_{CM}$, we will propose a particular
construction of a cylinder which will be characterized by other two
parameters, where one of them can be chosen as being the parameter $\chi$, and
the second one will be related to the parameter $\tilde{I}_{CM}$.

In the Fig. 2, we show a cross section of our cylinder model, which consists
of a thin cylindrical shell of mass $m_{C}$ with uniform mass distribution,
attached to a mass line $m_{P}$, parallel to the cylinder’s principal axis.
Besides the theoretical importance of this cylinder model that will be used in
the presentation of our main results, this model could also be used for an
experimental study of the jump effect.

Let us argue that this cylinder model encompasses all possible cases of
general eccentric cylinders characterized by the parameters $\chi$ and
$\tilde{I}_{CM}$. Note that the mass of the cylinder model is given by
$M=m_{C}+m_{P}$. Aside the parameter $\chi=d/R$, where $R$ is the radius of
the thin cylindrical shell, we introduce the parameter $k_{m}=m_{C}/M$.

Figure 2: Cross section of a thin cylindrical shell of mass $m_{C}$ plus a
mass line $m_{P}$ parallel to the cylinder’s principal axis.

For example, some particular cases are:

  * •

One where $CM$ coincides with $C$, so we have $d=0$ (namely $\chi=0$), this
case happens when the mass line $m_{P}$ passes through $C$. The parameter
$k_{m}$ controls the different possibilities for the values of the moment of
inertia.

  * •

One where $CM$ is over the border of the cylinder, it means $d=R$, this case
happens when $\chi=1$ and $k_{m}=0$.

Through the relation that defines the $CM$, we can show that the distance $s$
between the $CM$ and the mass line $m_{P}$, is given by
$s=\frac{k_{m}}{1-k_{m}}\chi R$. Then the moment of inertia $I_{C}$ is
computed as follows

$\displaystyle
I_{C}=m_{C}R^{2}+m_{P}(d+s)^{2}=m_{C}R^{2}+\frac{M^{2}}{m_{P}}d^{2}\,,$
$\displaystyle\Rightarrow\frac{I_{C}}{MR^{2}}=k_{m}+\frac{\chi^{2}}{1-k_{m}}\,.$
(17)

Since $Md^{2}\leq I_{C}\leq MR^{2}$, we have that:
$\chi^{2}\leq\frac{I_{C}}{MR^{2}}\leq 1$, and therefore $0\leq\chi\leq 1$.
Then from Eq. (17), we have that $0\leq k_{m}\leq 1-\chi$. Moreover, according
to Eq. (17) for the different values of $\chi$ and $k_{m}$ in the former
intervals, we cover all the possible values for $\frac{I_{C}}{MR^{2}}$. From
Eq. (17), we have that
$\tilde{I}_{CM}=\frac{I_{CM}}{MR^{2}}=k_{m}+\frac{k_{m}\chi^{2}}{1-k_{m}}$,
where this parameter $\tilde{I}_{CM}$ appears in Eq. (16), which explicitly
exhibits the independence of the dynamics in relation to the mass $M$ and the
radius $R$ when $\chi$ and $k_{m}$ are fixed.

##  3 Conditions for a slip-free transition from pure rolling to flight motion

As mentioned in the previous section, the cylinder whole motion is basically
composed of three types of particular motions: pure rolling, rolling plus
slipping and flight motion. From results of numerical simulation (work in
progress [22]), we have observed that depending on the initial conditions, the
values of the parameters $\chi$, $k_{m}$ and the ramp inclination $\alpha$,
the cylinder performs a variety of interesting motions. The most common
situation corresponds to the case where the cylinder initially performs a pure
rolling motion, then the motion is alternated between pure rolling and rolling
plus slipping, until at some point the cylinder loses contact with the ramp
and jumps.

In this section, we will study the following sequence of motions for the
cylinder and the conditions required to have such a sequence.

$\displaystyle\theta=\theta_{0}\;\;$ $\displaystyle\Rightarrow\text{initial
position,}$ $\displaystyle\theta_{0}<\theta<\theta_{J}\;\;$
$\displaystyle\Rightarrow\text{pure rolling motion,}$
$\displaystyle\theta=\theta_{J}\;\;$ $\displaystyle\Rightarrow\text{jump
position,}$ $\displaystyle\theta>\theta_{J}\;\;$
$\displaystyle\Rightarrow\text{flight motion.}$

Note that the transition from pure rolling to flight motion happens at the
point where $\theta=\theta_{J}$. The subscript $J$ means that we are
considering the value of the quantities at the instant of the jump. From here
to the rest of the paper, we will use the following acronym JARM to mean: jump
after a pure rolling motion.

As discussed in section 2, at the point where the cylinder jumps, the angle
$\theta$ and the angular velocity $\dot{\theta}$ must satisfy Eq. (14). We
have denoted by $\theta_{J}$ the angle that is a solution of Eq. (14).
Therefore from this equation, we can write

$\displaystyle\dot{\theta}_{J}^{2}=\frac{g\cos\alpha}{d\cos\theta_{J}}.$ (18)

Substituting Eq. (18) into the equation of motion Eq. (11), we obtain

$\displaystyle\ddot{\theta}_{J}=\frac{g(1+\chi\cos\theta_{J})\sec\theta_{J}\sin(\alpha+\theta_{J})}{R\left(k_{m}+\frac{\chi^{2}}{1-k_{m}}+1+2\chi\cos\theta_{J}\right)}.$
(19)

In order to have a JARM, since by definition the cylinder does not slip, the
force of friction between the cylinder and the ramp should be static, and

$\displaystyle\Big{|}\frac{F_{x}}{F_{y}}\Big{|}\leq\mu_{s}\,,$ (20)

where $F_{x}$ and $F_{y}$ are the friction and normal force, respectively.
Remember that the definition of these forces as given in Eqs. (9) are valid
for general rolling motions (pure rolling or rolling plus slipping). In order
to restrict our analysis to the pure rolling motion, we should substitute the
equations $\ddot{b}=0$ and $\ddot{l}=R\,\ddot{\theta}$, into Eqs. (9); after
such substitution, we obtain

$\displaystyle F_{x}$
$\displaystyle=MR\big{(}-\frac{g}{R}\sin\alpha-\chi\dot{\theta}^{2}\sin\theta+(1+\chi\cos\theta)\ddot{\theta}\,\big{)}\,,$
$\displaystyle F_{y}$
$\displaystyle=MR\big{(}\frac{g}{R}\cos\alpha-\chi\dot{\theta}^{2}\cos\theta-\chi\ddot{\theta}\sin\theta\big{)}\,.$
(21)

It is interesting to note that at the beginning of the cylinder’s motion, the
kinetic energy is low (or zero if the cylinder is left from rest), and the
$CM$ of the cylinder is in a higher position so that this initial
configuration allows the cylinder to roll down. Clearly at the beginning,
$F_{y}$ is positive because it is dominated by the term $Mg\cos\alpha$,
subsequently the values of $\dot{\theta}$ and $\ddot{\theta}$ grow and so the
value of the normal force approaches to zero and at some point becomes
negative. Since the cylinder is not attached to the ramp, negative values of
the normal force are not physical allowed in our study.

Before the normal force becomes negative (when it is positive and approaching
zero), the inequality given by Eq. (20) is satisfied. However, at some point,
this inequality will no longer hold, which means that the cylinder will begin
to slip. Numerical inspection [22] has revealed that the point where
inequality (20) is no longer valid is located in proximity to two other
points: the point where the normal force becomes zero and the point where the
jump condition is satisfied (as given by Eq. (14)). Based on these
observations, and in order to have a JARM in the case of pure rolling motion,
we will analyze the following condition for the normal force:

$\displaystyle F_{y}>0.$ (22)

Since this condition is less restrictive than the condition given by Eq. (20),
it is clear that we will need to complement Eq. (22) with Eq. (20).

###  3.1 Conditions for a JARM independent of the initial conditions

In order to perform a general analysis of the values and signs of the friction
and normal force at the jump point, namely when $\theta=\theta_{J}$, we
substitute Eqs. (18) and (19) into Eqs. (3), so that we obtain

$\displaystyle F_{x,J}$
$\displaystyle=-\frac{gM\left(k_{m}+\frac{\chi^{2}}{1-k_{m}}-\chi^{2}\cos^{2}\theta_{J}\right)\sec\theta_{J}\sin(\alpha+\theta_{J})}{k_{m}+\frac{\chi^{2}}{1-k_{m}}+1+2\chi\cos\theta_{J}}\,,$
(23) $\displaystyle F_{y,J}$
$\displaystyle=-\frac{gM\chi(1+\chi\cos\theta_{J})\tan\theta_{J}\sin(\alpha+\theta_{J})}{k_{m}+\frac{\chi^{2}}{1-k_{m}}+1+2\chi\cos\theta_{J}}\,.$
(24)

Regarding the result of the normal force $F_{y,J}$ from equation (24), we note
that $F_{y,J}$ is not necessarily equal to zero. In the case where
$F_{y,J}>0$, and since for $\theta>\theta_{J}$ we have $F_{y}=0$, it is clear
that the normal force changes discontinuously from a non-vanishing to zero
value at the point where the cylinder jumps.

Note that some terms in Eqs. (23) and (24) are positive, these are:
$k_{m}+\frac{\chi^{2}}{1-k_{m}}+1+2\chi\cos\theta_{J}>0$; and
$1+\chi\cos\theta_{J}>0$. Also from Eq. (18) we have: $\cos\theta_{J}>0$; and
thus $k_{m}+\frac{\chi^{2}}{1-k_{m}}-\chi^{2}\cos^{2}\theta_{J}>0$. After
consideration of the positivity condition of these terms, employing Eqs. (23)
and (24), we get

$\displaystyle sign(F_{x,J})$ $\displaystyle=-sign(\sin(\alpha+\theta_{J})),$
(25) $\displaystyle sign(F_{y,J})$
$\displaystyle=-sign(\sin\theta_{J}\sin(\alpha+\theta_{J})).$ (26)

A subtle issue that can be observed from equation (26) is that the normal
force at the jump point could even be negative. Physically, this would mean
that before the cylinder jumps, the normal force could have been zero.
Therefore, the cylinder may have slipped before jumping, implying that the
pure rolling assumption is no longer valid. In what follows, we will address
this issue in more detail.

As mentioned before, in order to have a JARM, the condition that the normal
force should be positive needs to be complemented with Eq. (20). Therefore,
let us start our analysis by searching for possible allowed values of
$\theta_{J}$ through Eq. (26), so that the condition given by Eq. (22) is
satisfied. In that sense, it is not difficult to show that

$\displaystyle 2\pi n-\alpha<\theta_{J}<2\pi n,\;\;\;n=1,2,3,\dots$ (27)

From Fig. 3a, we can see that the angle $\alpha+\theta$ is measured between
$\vec{d}$ and the vertical line. Therefore from this geometrical configuration
of the angles, we can easily prove that the whole shadowed regions (gray and
green) corresponds to the regions where $\cos\theta>0$. Now from Eq. (26), it
is not difficult to see that in the gray region $F_{y,J}<0$, and in the green
region $F_{y,J}>0$.

Therefore, the only domain where a JARM could happen (due to the condition in
Eq. (22)), is when the $CM$ is inside the green region, on the other regions a
JARM is precluded. Notice that the green region agrees with the interval for
$\theta_{J}$ given by Eq. (27), where the integer number $n$ is interpreted as
the counting of the full turns completed by the cylinder. Finally, in the
green region we can check that $F_{x,J}<0$ (this result is obtained from Eq.
(25)), which means that the static friction force over the cylinder points
along and upward the ramp.

Figure 3: (a) Representation of the angles $\alpha$ and $\theta$ on the cross
section of the cylinder. (b) Plot of the function
$\mu(\chi,k_{m},\theta_{J})$, for $\chi=0.5$ and $k_{m}=0.4$.

More restrictions for the possible allowed values of $\theta_{J}$ can be
obtained from the condition given by Eq. (20). So, substituting Eqs. (23) and
(24) into Eq. (20) we obtain

$\displaystyle\mu(\chi,k_{m},\theta_{J})\equiv\frac{\frac{-k_{m}^{2}+k_{m}+\chi^{2}}{(k_{m}-1)\chi^{2}}+\cos^{2}\theta_{J}}{\sin\theta_{J}\,\left(\frac{1}{\chi}+\cos\theta_{J}\right)}\leq\mu_{s}\,.$
(28)

In the Fig. 3b, a typical plot of the function $\mu(\chi,k_{m},\theta_{J})$ is
shown. By fixing the value of $\mu_{s}$, and using the inequality given in Eq.
(28), we can calculate the region to which $\theta_{J}$ belongs. It is worth
noting that the domain of $\theta_{J}$ as shown in Fig. 3b is bigger than the
one defined by Eq. (27), however it is important to keep in mind that the
lower bound for $\theta_{J}$ must be greater than or equal to $2\pi n-\alpha$.
We also note that due to the usual periodic property of the trigonometric
functions that appear in $\mu(\chi,k_{m},\theta_{J})$, the plot of this
function is the same for any integer value $n$. After some analysis of the
Fig. 3b, we conclude that to have a JARM:

  1. 1.

There is a minimum value $\mu_{s-min}$ for the coefficient of static friction.

  2. 2.

For each $\mu_{s}>\mu_{s-min}$, we have that $max\\{2\pi
n-\alpha,\theta_{J-min}\\}<\theta_{J}<\theta_{J-max}$; where $\theta_{J-min}$
and $\theta_{J-max}$ are identified as follows, when $\mu_{s}>\mu_{s-min}$ and
$\mu_{s}\approx\mu_{s-min}$, there are two solutions to the equation
$\mu(\chi,k_{m},\theta_{J})=\mu_{s}$, these solutions are precisely
$\theta_{J-min}$ and $\theta_{J-max}$, while for values of $\mu_{s}$ such that
$\mu_{s}\gg\mu_{s-min}$, there is only one solution that is $\theta_{J-max}$.

  3. 3.

Since there is a maximum value for $\theta_{J}$, which was denoted by
$\theta_{J-max}$, then there is a minimum value for $\alpha$, denoted by
$\alpha_{min}$ which satisfies $\alpha_{min}=2\pi n-\theta_{J-max}$, this is
true because $\theta_{J}$ satisfies the inequality given in Eq. (27).

As a specific application of the general results presented above, let us
consider the case where $\alpha=0$, which corresponds to the horizontal plane.
From Eq. (27), we can see that the only valid value for the angle $\theta_{J}$
in this case is $\theta_{J}=2\pi n$. However, substituting this value into the
left-hand side of Eq. (28) results in a divergence, indicating that to
maintain roll without slipping, the coefficient of static friction $\mu_{s}$
must approach infinity. Since such a value is physically impossible, we can
conclude that the no-slip condition must be violated before the cylinder
jumps, regardless of the values of $\chi$ and $k_{m}$. In the following
discussion, we will consider the case where $\alpha\neq 0$.

###  3.2 Conditions for a JARM using initial conditions

Our study in the previous subsection is interesting because it proves that
there are restrictions for the possible allowed values of $\alpha$, $\mu_{s}$
and $\theta_{J}$. Also note that to derive our previous results, we have not
used any particular initial condition.

In this subsection, we consider the following quite standard initial
conditions: the cylinder is left from rest $\dot{\theta}_{0}=0$ at the top of
the ramp, with $\theta_{0}=0$. As we are going to see, the setting of these
initial conditions will further restrict the existence of a JARM.

Using the conservation of energy, which is valid for pure rolling motion,
together with the initial conditions $\dot{\theta}_{0}=0$ and $\theta_{0}=0$,
from the jump condition given by Eq. (14), we can show that the angle
$\theta_{J}$ at the jump point is the root of the function:

$\displaystyle
J(\theta)=-1+\frac{\big{(}\chi\cos\alpha-\chi\cos(\alpha+\theta)+\theta\sin\alpha\big{)}\cos\theta}{(2\pi\gamma+\cos\theta-1)\cos\alpha},$
(29)

where the parameter $\gamma$ is defined by

$\displaystyle\gamma=\frac{k_{m}^{2}+2k_{m}\chi-\chi^{2}-2\chi-1}{4\pi(k_{m}-1)\chi}.$
(30)

Let us provide conditions for the existence of roots $\theta_{J}$ of the
function $J(\theta)$ with the restriction given by equation (27). It turns out
that to guarantee the existence of a root $\theta_{J}\in[2\pi n-\alpha,2\pi
n]$ of the function $J(\theta)$, we must have that

$\displaystyle J(2\pi n-\alpha)$
$\displaystyle=\Big{(}\frac{\chi\cos\alpha+(2\pi
n-\alpha)\sin\alpha-\chi}{\cos\alpha+2\pi\gamma-1}-1\Big{)}<0,\;\;\;\text{and},$
(31) $\displaystyle J(2\pi n)$
$\displaystyle=\Big{(}\frac{n\tan\alpha}{\gamma}-1\Big{)}>0.$ (32)

Note that through the equation: $J(\theta_{J})=0$, it should be possible to
express $\theta_{J}$ as a function of the parameters $\alpha$, $\chi$, and
$k_{m}$, namely

$\displaystyle\theta_{J}=\theta_{n,J}(\alpha,\chi,k_{m}),$ (33)

where the subscript $n$ explicitly indicates that $\theta_{J}$ belongs to the
interval $[2\pi n-\alpha,2\pi n]$. Using this equation (33), we can write the
inequality (28) as follows

$\displaystyle\frac{\frac{-k_{m}^{2}+k_{m}+\chi^{2}}{(k_{m}-1)\chi^{2}}+\cos^{2}\big{(}\theta_{n,J}(\alpha,\chi,k_{m})\big{)}}{\sin\big{(}\theta_{n,J}(\alpha,\chi,k_{m})\big{)}\,\Big{[}\frac{1}{\chi}+\cos\big{(}\theta_{n,J}(\alpha,\chi,k_{m})\big{)}\Big{]}}\leq\mu_{s}\,.$
(34)

We could think that for given values of $n$, $\alpha$, $\chi$, $k_{m}$ and
$\mu_{s}$ such that inequalities (31), (32) and (34) are true, it would be
enough to guarantee that the cylinder will jump without first having slipped.
However, note that the above analysis has been performed at the point where
the cylinder jumps, namely the inequalities (31), (32) and (34) do not
necessarily guarantee pure rolling motion for values of the angle $\theta$
such that $\theta<\theta_{J}$. Therefore we will need to impose more
restrictions.

Since extra conditions will come from the analysis of inequalities of the type
given in (22) and (20), we need to write $F_{x}$ and $F_{y}$ for generic
values of $\theta$ such that $\theta<\theta_{J}$. These components of the
force are given by Eqs. (3).

For the initial conditions $\dot{\theta}_{0}=0$ and $\theta_{0}=0$, using the
equation of motion (11), and the conservation of the energy, we can express
the angular velocity $\dot{\theta}$ and the angular acceleration
$\ddot{\theta}$ in terms of the angle $\theta$, so that the components of the
force given in Eqs. (3) can be written as functions that depend explicitly on
$\theta$

$\displaystyle F_{x}(\theta)=$
$\displaystyle-\frac{gM}{4\chi(2\pi\gamma+\cos\theta-1)^{2}}\Big{[}\sin\alpha\big{(}16\pi^{2}\gamma^{2}\chi-16\pi\gamma\chi-4\pi\gamma-\chi^{2}+3\chi+2$
$\displaystyle+2\pi\gamma\chi^{2}+2(-1+(-3+6\pi\gamma)\chi)\cos\theta+\chi\cos(2\theta)-2\theta\sin\theta-4\theta\chi\sin\theta$
$\displaystyle+8\pi\gamma\theta\chi\sin\theta+\theta\chi\sin(2\theta)\big{)}+\chi(2\cos\alpha(-1-2\chi+4\pi\gamma\chi+\chi\cos\theta)\sin\theta$
$\displaystyle-(-2+4\pi\gamma+\chi)\sin(\alpha+\theta)-\chi((-3+6\pi\gamma)\sin(\alpha+2\theta)+\sin(\alpha+3\theta)))\Big{]},$
(35) $\displaystyle F_{y}(\theta)=$
$\displaystyle\;\;\frac{gM}{4(2\pi\gamma+\cos\theta-1)^{2}}\Big{[}\cos\alpha\big{(}6-16\pi\gamma+16\pi^{2}\gamma^{2}-4\chi+2\pi\gamma\chi$
$\displaystyle+(-8-8\pi\gamma(-2+\chi)+7\chi)\cos\theta+(2-4\chi+6\pi\gamma\chi)\cos(2\theta)+\chi\cos(3\theta)\big{)}$
$\displaystyle-\sin\alpha(3\theta+4(-1+2\pi\gamma)\theta\cos\theta+\theta\cos(2\theta)-2\sin\theta+4\pi\gamma\sin\theta$
$\displaystyle+3\chi\sin\theta+\sin(2\theta)-3\chi\sin(2\theta)+6\pi\gamma\chi\sin(2\theta)+\chi\sin(3\theta))\Big{]}.$
(36)

In order to guarantee a pure rolling motion throughout the entire path from
$\theta=0$ to the point where the cylinder jumps $\theta=\theta_{J}$, for any
value of $\theta$ such that $0<\theta<\theta_{J}$, the following inequalities
must be satisfied

$\displaystyle F_{y}(\theta)$ $\displaystyle>0,\;\;\;\text{and},$ (37)
$\displaystyle\frac{|F_{x}(\theta)|}{F_{y}(\theta)}$
$\displaystyle\leq\mu_{s}.$ (38)

Let us summarize the main result of this subsection. We have shown that using
the initial conditions $\theta_{0}=0$ and $\dot{\theta}_{0}=0$, to have a JARM
the parameters $n$, $\alpha$, $\chi$, $k_{m}$ and $\mu_{s}$ need to be chosen
such that the inequalities (31), (32), (34), (37) and (38) are satisfied.
Therefore, this result imposes nontrivial restrictions on the possible allowed
values of the parameters that appear in the equations of the problem.

###  3.3 Regions in the parameter space ($\alpha$,$\chi$,$k_{m}$) for fixed
values of $n$ and $\mu_{s}$

In order to have a JARM for the initial conditions $\theta_{0}=0$ and
$\dot{\theta}_{0}=0$, in this subsection, we are going to show the region
where the parameters must belong. Since essentially we have five parameters
$n$, $\mu_{s}$, $\alpha$, $\chi$, and $k_{m}$ to visualize the region defined
by the inequalities (31), (32), (34), (37) and (38), we will need to fix at
least two parameters so that the remaining three parameters can be visualized
in a three-dimensional space.

By fixing the value of the parameters $n$ and $\mu_{s}$, in Fig. 4 we present
the regions in the parameter space $(\alpha,\chi,k_{m})$ where the
inequalities mentioned in the previous paragraph are fulfilled, namely if we
choose any set of parameters $\alpha$, $\chi$, and $k_{m}$ that belong to
these regions, we guarantee the occurrence of a JARM. Comparing Fig. 4(b) and
Fig. 4(a), we see that the regions with $\mu_{s}=1$ are bigger than the
regions with $\mu_{s}=0.7$. This result makes sense since for a larger static
coefficient of friction, we expect that the cylinder has more chances to
maintain pure rolling motion.

For a fixed value of $\mu_{s}$, we can also compare the regions obtained with
different values of $n$. For instance, from Fig. 4(a), for the value of
$\mu_{s}=0.7$, we observe that as the values of $n$ increase, the
corresponding regions are getting smaller. In general, for any value of
$\mu_{s}$, this pattern was observed. The physical interpretation of this
result is as follows, first let us remember that $n$ represents the number of
turns performed by the cylinder, so for a given value of $\mu_{s}$, since
$n=1$ is the lowest possible value for $n$, the greatest chance of having a
JARM occurs before the cylinder completes a full turn, and the chance
decreases every time we increase the value of $n$.

Regarding the inclination of the ramp given by the angle $\alpha$ (where
$\pi/2$ is its maximum value), we observed that, inside the regions shown in
Fig. 4, large values of $\alpha$ are in correspondence with small values of
$\chi$. This type of correspondence makes sense since for a high inclination
of the ramp, in order to avoid a slipping, the normal force should not
oscillate to much, and that happens when the value of $\chi$ is small.

When the inclination of the ramp gets closer to the value of $\pi/2$, we
notice that the regions defined by the values of $\chi$ and $k_{m}$ becomes
smaller, where in the limit case $\alpha\rightarrow\pi/2$ the parameters
$\chi$ and $k_{m}$ vanish. This is evidenced by the fact that when $\alpha$ is
close to $\pi/2$, the region where a JARM occurs has the shape of a wedge
where its vertex is given by the point $\alpha=\pi/2$, $\chi=0$ and $k_{m}=0$.
In order to interpret this result, let us remember that inside the region
where a JARM happens, large values of $\alpha$ are in correspondence with
small values of $\chi$, namely when the $CM$ is close to the geometric center
$C$ of the cylinder, there are less oscillations of the normal force, which
implies in less chances to have a slip. We conclude that high values of
$\alpha$ are allowed provided that the parameters $\chi$ and $k_{m}$ are small
enough. Therefore, an unexpected result of having a JARM is obtained in the
limit case where $\alpha\rightarrow\pi/2$, $\chi\rightarrow 0$ and
$k_{m}\rightarrow 0$.

Figure 4: Regions where a JARM happens in the space of parameters $\alpha$,
$\chi$ and $k_{m}$, for fixed values of $\mu_{s}$, $n$ and the initial
conditions $\theta_{0}=0$ and $\dot{\theta}_{0}=0$. (a) Here we have set
$\mu_{s}=0.7$ and $n=1,2,3$ correspond to the green, purple and orange solids
respectively. (b) This plot is the same as the one presented in (a) but with
$\mu_{s}=1$.

##  4 Summary and discussion

For a given value of the coefficient of static friction $\mu_{s}$, we have
shown that a general eccentric cylinder performs a jump starting from pure
rolling motion, provided that the angle $\alpha$, and the parameters $\chi$
and $k_{m}$ that characterize the cylinder, belong to a restricted region. If
these parameters do not belong to the aforementioned region, the cylinder has
to perform another type of motion, such as slipping with rolling, before the
jump. In a future paper [22], we will analyze these other varieties of motion
using our general cylinder.

Another important issue that can be explored is related to the initial
conditions. We have presented a general discussion about the existence of JARM
that is independent of the initial conditions. As a result of this analysis,
we show that the value of the angle $\theta_{J}$ is restricted to some
interval. To fix the value of $\theta_{J}$, some particular initial conditions
are needed. Therefore, we have used the somewhat standard initial conditions,
$\theta_{0}=0$ and $\dot{\theta}_{0}=0$. It will be an interesting and non-
trivial problem to analyze how the regions shown in Fig. 4 change when other
initial conditions are used. For example, we can see if these regions increase
or decrease, which would physically imply a greater or lesser chance of having
JARM. What would be the optimal initial conditions that allow a greater chance
of having JARM?

Finally, since only the slipping motion before the jump has been observed in
experiments and theoretical studies carried out to date, it has been a common
conclusion that the no-slip conditions must be violated before the jump.
Indeed, when $\alpha=0$, on general grounds, we have definitely shown that
this last conclusion is true. However, in the case where $\alpha\neq 0$, there
is a chance to have JARM. Therefore, to empirically test the occurrence of
JARM, it would be important to set up an experiment that takes into account
appropriate values for $\mu_{s}$ and the parameters $\chi$ and $k_{m}$.

## Acknowledgements

We would like to thank Dominique Sugny, Gabriela and Pavao Mardešić for useful
discussions.

## References

  * [1] J. E. Littlewood, A mathematician’s miscellany.  Methuen, 1953.
  * [2] T. F. Tokieda, “The hopping hoop,” The American Mathematical Monthly, vol. 104, no. 2, pp. 152–154, 1997.
  * [3] J. P. Butler, “Hopping hoops don’t hop,” The American Mathematical Monthly, vol. 106, no. 6, pp. 565–568, 1999.
  * [4] T. Pritchett, “The hopping hoop revisited,” The American Mathematical Monthly, vol. 106, no. 7, pp. 609–617, 1999.
  * [5] L. Yanzhu and X. Yun, “Qualitative analysis of a rolling hoop with mass unbalance,” Acta Mechanica Sinica, vol. 20, no. 6, pp. 672–675, 2004.
  * [6] W. F. D. Theron, “The rolling motion of an eccentrically loaded wheel,” American Journal of Physics, vol. 68, no. 9, pp. 812–820, 2000.
  * [7] A. Taylor and M. Fehrs, “The dynamics of an eccentrically loaded hoop,” American Journal of Physics, vol. 78, no. 5, pp. 496–498, 2010.
  * [8] A. Bronars and O. M. O’Reilly, “Gliding motions of a rigid body: the curious dynamics of littlewood’s rolling hoop,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 475, no. 2231, p. 20190440, 2019.
  * [9] M. Batista, “Self-induced jumping of a rigid body of revolution on a smooth horizontal surface,” International Journal of Non-Linear Mechanics, vol. 43, no. 1, pp. 26–35, 2008.
  * [10] A. P. Ivanov, “On detachment conditions in the problem on the motion of a rigid body on a rough plane,” Regular and Chaotic Dynamics, vol. 13, no. 4, pp. 355–368, 2008.
  * [11] Y. Shimomura, M. Branicki, and H. Moffatt, “Dynamics of an axisymmetric body spinning on a horizontal surface. ii. self-induced jumping,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 461, no. 2058, pp. 1753–1774, 2005.
  * [12] P. Kessler and O. M. O’Reilly, “The ringing of euler’s disk,” Regular and Chaotic dynamics, vol. 7, no. 1, pp. 49–60, 2002.
  * [13] R. I. Leine, “Experimental and theoretical investigation of the energy dissipation of a rolling disk during its final stage of motion,” Archive of Applied Mechanics, vol. 79, no. 11, pp. 1063–1082, 2009.
  * [14] W. Theron and M. Maritz, “The amazing variety of motions of a loaded hoop,” Mathematical and Computer Modelling, vol. 47, no. 9, pp. 1077–1088, 2008\.
  * [15] P. Onorato, M. Malgieri, P. Mascheretti, and A. D. Ambrosis, “The surprising rolling spool: librational motion and failure of the pure rolling condition,” European Journal of Physics, vol. 36, p. 038002, mar 2015.
  * [16] S. A. Moore, D. Culver, and B. P. Mann, “The eccentric disk and its eccentric behavior,” European Journal of Physics, vol. 42, p. 065012, oct 2021.
  * [17] A. Carnevali and R. May, “Rolling motion of non-axisymmetric cylinders,” American Journal of Physics, vol. 73, no. 10, pp. 909–913, 2005.
  * [18] L. Turner and A. M. Turner, “Asymmetric rolling bodies and the phantom torque,” American Journal of Physics, vol. 78, no. 9, pp. 905–908, 2010\.
  * [19] J. H. Jensen, “Rules for rolling as a rotation about the instantaneous point of contact,” European Journal of Physics, vol. 32, pp. 389–397, jan 2011\.
  * [20] B. Y.-K. Hu, “Rolling of asymmetric discs on an inclined plane,” European Journal of Physics, vol. 32, pp. L51–L54, oct 2011.
  * [21] R. W. Gómez, J. J. Hernández-Gómez, and V. Marquina, “A jumping cylinder on an inclined plane,” European Journal of Physics, vol. 33, pp. 1359–1365, jul 2012.
  * [22] E. A. Arroyo and et. al., “work in progress.,”


# Perceptual Quality Enhancement of Sound Field Synthesis Based on Combination
of Pressure and Amplitude Matching

###### Abstract

A sound field synthesis method enhancing perceptual quality is proposed. Sound
field synthesis using multiple loudspeakers enables spatial audio reproduction
with a broad listening area; however, synthesis errors at high frequencies
called spatial aliasing artifacts are unavoidable. To minimize these
artifacts, we propose a method based on the combination of pressure and
amplitude matching. On the basis of the human’s auditory properties,
synthesizing the amplitude distribution will be sufficient for horizontal
sound localization. Furthermore, a flat amplitude response should be
synthesized as much as possible to avoid coloration. Therefore, we apply
amplitude matching, which is a method to synthesize the desired amplitude
distribution with arbitrary phase distribution, for high frequencies and
conventional pressure matching for low frequencies. Experimental results of
numerical simulations and listening tests using a practical system indicated
that the perceptual quality of the sound field synthesized by the proposed
method was improved from that synthesized by pressure matching.

Index Terms—  sound field synthesis, sound field reproduction, pressure
matching, amplitude matching, spatial audio

##  1 Introduction

Sound field synthesis/reproduction is a technique to recreate a spatial sound
using multiple loudspeakers (or secondary sources). One of its prospective
applications is spatial audio for virtual/augmented reality, which enables
spatial audio reproduction with a broader listening area than in the case of
conventional spatial audio techniques, such as multichannel surround sound and
binaural synthesis.

Sound field synthesis methods based on the minimization of the squared error
between synthesized and desired sound fields, such as pressure matching and
mode matching [1, 2, 3, 4, 5, 6], have practical advantages compared with
methods based on analytical representations derived from boundary integral
equations such as wave field synthesis and higher-order ambisonics [7, 8, 9,
4, 10, 11, 12, 13], since the array geometry of the secondary sources can be
arbitrary. In particular, pressure matching is widely used because of its
simple implementation.

A well-known issue of the sound field synthesis methods is spatial aliasing
artifacts. That is, depending on the secondary source placement, the synthesis
accuracy can significantly decrease at high frequencies, which can lead to the
degradation of sound localization for listeners and distortion of timbre, or
coloration, of source signals. Thus, the perceptual quality of the synthesized
sound field can considerably deteriorate.

To improve the perceptual quality, we propose a method combining amplitude
matching, which was originally proposed for multizone sound field control [14,
15], with pressure matching. Amplitude matching is a method to synthesize the
desired amplitude (or magnitude) distribution, leaving the phase distribution
arbitrary, whereas pressure matching aims to synthesize both amplitude and
phase distributions, i.e., pressure distribution. We apply amplitude matching
to mitigate the spatial aliasing artifacts by reducing the parameters to be
controlled at high frequencies, keeping the range of the listening area broad.
On the basis of the human’s auditory properties, the interaural level
difference (ILD) is known to be a dominant cue for horizontal sound
localization at high frequencies, compared with the interaural time difference
(ITD) [16, 17]. Therefore, synthesizing the amplitude distribution will be
sufficient for sound localization. Furthermore, by prioritizing the synthesis
of the desired amplitude distribution, a flat amplitude response is reproduced
as much as possible, and coloration effects can be alleviated. We formulate a
new cost function combining amplitude matching for high frequencies and
conventional pressure matching for low frequencies, which can be solved in a
similar manner to amplitude matching. We evaluate the proposed method through
numerical experiments in the frequency domain and listening experiments in a
real environment.

##  2 Sound Field Synthesis Problem

Let $\Omega\subset\mathbb{R}^{3}$ be a target region for synthesizing the
desired sound field. As shown in Fig. 1, $L$ secondary sources are placed at
$\bm{r}_{l}\in\mathbb{R}^{3}\backslash\Omega$ ($l\in\\{1,\ldots,L\\}$). The
driving signal of the $l$th secondary source and its transfer function to the
position $\bm{r}\in\Omega$ at the angular frequency $\omega$ are denoted as
$d_{l}(\omega)$ and $g_{l}(\bm{r},\omega)$, respectively. Then, the
synthesized pressure distribution $u_{\mathrm{syn}}(\bm{r},\omega)$
($\bm{r}\in\Omega$) is represented as

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

where $\bm{g}(\bm{r},\omega)\in\mathbb{C}^{L}$ and
$\bm{d}(\omega)\in\mathbb{C}^{L}$ are the vectors consisting of
$g_{l}(\bm{r},\omega)$ and $d_{l}(\omega)$, respectively. Hereafter, $\omega$
is omitted for notational simplicity.

The objective is to synthesize the desired sound field
$u_{\mathrm{des}}(\bm{r})$ over the target region $\Omega$. We define the
optimization problem to obtain $\bm{d}$ as follows.

$\displaystyle\mathop{\rm
minimize}\limits_{\bm{d}\in\mathbb{C}^{L}}Q(\bm{d}):=\int_{\bm{r}\in\Omega}\left|\bm{g}(\bm{r})^{\mathsf{T}}\bm{d}-u_{\mathrm{des}}(\bm{r})\right|^{2}\mathrm{d}\bm{r}$
(2)

Since this problem is difficult to solve owing to the regional integration,
several methods to approximately solve it, for example, (weighted) pressure
and mode matching [1, 2, 3, 4, 5, 6], have been proposed.

Figure 1: Sound field synthesis over the target region $\Omega$ with multiple
loudspeakers

###  2.1 Pressure Matching

Pressure matching is one of the widely used sound field synthesis methods to
approximately solve (2). First, the region $\Omega$ is discretized into $N$
control points whose positions are denoted as $\\{\bm{r}_{n}\\}_{n=1}^{N}$. We
assume that the control points are arranged densely enough over $\Omega$.
Then, the cost function for pressure matching is formulated as the
minimization problem of the squared error between the synthesized and desired
pressures at the control points as

$\displaystyle\mathop{\rm minimize}\limits_{\bm{d}\in\mathbb{C}^{L}}\
\|\bm{G}\bm{d}-\bm{u}_{\mathrm{des}}\|_{2}^{2}+\beta\|\bm{d}\|_{2}^{2},$ (3)

where $\bm{G}\in\mathbb{C}^{N\times L}$ is the matrix consisting of transfer
functions $\\{g_{l}(\bm{r}_{n})\\}_{n=1}^{N}$,
$\bm{u}_{\mathrm{des}}\in\mathbb{C}^{N}$ is the vector consisting of the
desired pressures $\\{u_{\mathrm{des}}(\bm{r}_{n})\\}_{n=1}^{N}$, and $\beta$
is the regularization parameter. This least-squares problem (3) has a closed-
form solution as follows:

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

Pressure matching is extended as weighted pressure matching [18, 6] by
combining with sound field interpolation techniques [19, 20].

Another strategy to approximate $Q(\bm{d})$ is to represent the sound field by
spherical wavefunction expansion [21, 22], which is referred to as (weighted)
mode matching [4, 5]. It is demonstrated that weighted pressure matching is a
special case of weighted mode matching in [6].

###  2.2 Spatial Aliasing Artifacts

On the basis of the single-layer potential [23, 21], the desired sound field
can be perfectly synthesized if secondary sources are continuously distributed
point sources on a surface of the target region $\Omega$. However, owing to
the discrete placement of the secondary sources, spatial aliasing artifacts
are unavoidable for sound field synthesis methods. The synthesis accuracy can
decrease particularly at high frequencies, which can lead to the degradation
of sound localization and the coloration of source signals. The properties of
spatial aliasing in analytical sound field synthesis methods have been
extensively investigated [24], but spatial aliasing can also occur in
numerical methods, depending on the secondary source placement [25, 26]. One
of the strategies to improve the synthesis accuracy at high frequencies is to
make the target region smaller [27, 5]. The challenge here is to improve the
perceptual quality of the synthesized sound field by mitigating spatial
aliasing artifacts in a broad listening area.

##  3 Proposed Method Based on Combination of Pressure and Amplitude Matching

Even when it is difficult to synthesize the desired pressure distribution,
i.e., amplitude and phase, for high frequencies, it can be considered that
synthesizing only the amplitude distribution can be achieved by leaving the
phase distribution arbitrary. On the basis of the human’s auditory properties,
ILD is known to be a dominant cue for horizontal sound localization above
$1500~{}\mathrm{Hz}$, and the dependence on ITD is markedly reduced from
$1000$–$1500~{}\mathrm{Hz}$ [16, 17], which indicates that synthesizing the
amplitude distribution is sufficient for sound localization above
$1500~{}\mathrm{Hz}$. This perceptual property is also used in a method for
binaural rendering [28]. Furthermore, coloration effects can be alleviated by
reproducing the flat amplitude response as much as possible. Therefore, we
propose a method combining amplitude matching [14, 15], which is aimed to
synthesize the desired amplitude at the control points, with pressure
matching. By applying amplitude matching for high frequencies, the parameters
to be controlled can be reduced from pressure matching, keeping the range of
the target region, i.e., the number of control points. Thus, we can improve
the perceptual quality of the synthesized sound field by reproducing a more
accurate amplitude distribution over $\Omega$.

###  3.1 Proposed Algorithm

We define the optimization problem of the proposed method as a composite of
pressure and amplitude matching as

$\displaystyle\mathop{\rm
minimize}\limits_{\bm{d}\in\mathbb{C}^{L}}J(\bm{d}):=$
$\displaystyle\left(1-\gamma\right)\left\|\bm{G}\bm{d}-\bm{u}_{\mathrm{des}}\right\|_{2}^{2}$
$\displaystyle\ \
+\gamma\left\||\bm{G}\bm{d}|-|\bm{u}_{\mathrm{des}}|\right\|_{2}^{2}+\beta\|\bm{d}\|_{2}^{2},$
(5)

where $|\cdot|$ represents the element-wise absolute value of vectors, and
$\gamma\in\mathbb{R}([0,1])$ is the parameter that determines the balance
between pressure and amplitude matching. When $\gamma=0$, (5) corresponds to
pressure matching, and when $\gamma=1$, it corresponds to amplitude matching;
therefore, $\gamma$ should be set close to $0$ for low frequencies and $1$ for
high frequencies. For example, $\gamma$ can be defined as a sigmoid function
of $\omega$ with the transition angular frequency $\omega_{\mathrm{T}}$ and
parameter $\sigma$ as

$\displaystyle\gamma(\omega)=\frac{1}{1+\mathrm{e}^{-\frac{\sigma}{2\pi}(\omega-\omega_{\mathrm{T}})}}.$
(6)

Since the cost function $J$ in (5) is neither convex nor differentiable owing
to the squared error term of the amplitude matching, (5) does not have a
closed-form solution. However, the alternating direction method of multipliers
(ADMM) can be applied in a similar manner to the algorithm for amplitude
matching [15]. First, we introduce the auxiliary variables of amplitude
$\bm{a}\in\mathbb{R}_{\geq\bm{0}}^{N}$ and phase
$\bm{\theta}\in\mathbb{R}^{N}$ such that
$\bm{G}\bm{d}=\bm{a}\odot\mathrm{e}^{{\mathrm{j}}\bm{\theta}}$, where $\odot$
represents the Hadamard product. Then, (5) is reformulated as

$\displaystyle\mathop{\rm minimize}\limits_{\bm{d},\bm{a},\bm{\theta}}\
(1-\gamma)\left\|\bm{a}\odot\mathrm{e}^{{\mathrm{j}}\bm{\theta}}-\bm{u}_{\mathrm{des}}\right\|_{2}^{2}$
$\displaystyle\hskip
60.0pt+\gamma\left\|\bm{a}-|\bm{u}_{\mathrm{des}}|\right\|_{2}^{2}+\beta\|\bm{d}\|_{2}^{2}$
$\displaystyle\mathrm{subject\ to}\
\bm{G}\bm{d}=\bm{a}\odot\mathrm{e}^{{\mathrm{j}}\bm{\theta}}.$ (7)

The augmented Lagrangian function $L_{\rho}$ for (7) is defined as

$\displaystyle L_{\rho}(\bm{a},\bm{\theta},\bm{d},\bm{\lambda})$
$\displaystyle=(1-\gamma)\left\|\bm{a}\odot\mathrm{e}^{\mathrm{j}\bm{\theta}}-\bm{u}_{\mathrm{des}}\right\|_{2}^{2}+\gamma\left\|\bm{a}-|\bm{u}_{\mathrm{des}}|\right\|_{2}^{2}+\beta\|\bm{d}\|_{2}^{2}$
$\displaystyle\hskip
20.0pt+\Re\left[\bm{\lambda}^{\mathsf{T}}(\bm{G}\bm{d}-\bm{a}\odot\mathrm{e}^{\mathrm{j}\bm{\theta}})\right]+\frac{\rho}{2}\|\bm{G}\bm{d}-\bm{a}\odot\mathrm{e}^{\mathrm{j}\bm{\theta}}\|_{2}^{2},$
(8)

where $\bm{\lambda}\in\mathbb{C}^{N}$ is the Lagrange multiplier, $\Re[\cdot]$
represents the real part of a complex value, and $\rho>0$ is the penalty
parameter.

Each variable is alternately updated on the basis of ADMM as

$\displaystyle\left(\bm{a}^{(i+1)},\bm{\theta}^{(i+1)}\right)=\mathop{\rm
arg~{}min}\limits_{\bm{a},\bm{\theta}}\
L_{\rho}\left(\bm{a},\bm{\theta},\bm{d}^{(i)},\bm{\lambda}^{(i)}\right)$ (9)
$\displaystyle\bm{d}^{(i+1)}=\mathop{\rm arg~{}min}\limits_{\bm{d}}\
L_{\rho}\left(\bm{a}^{(i+1)},\bm{\theta}^{(i+1)},\bm{d},\bm{\lambda}^{(i)}\right)$
(10)
$\displaystyle\bm{\lambda}^{(i+1)}=\bm{\lambda}^{(i)}+\rho\left(\bm{G}\bm{d}^{(i+1)}-\bm{a}^{(i+1)}\odot\mathrm{e}^{\mathrm{j}\bm{\theta}^{(i+1)}}\right),$
(11)

where $i$ denotes the iteration index. (9) is minimized independently for
$\bm{\theta}$ and $\bm{a}$ as

$\displaystyle\bm{\theta}^{(i+1)}=\arg\left((1-\gamma)\bm{u}_{\mathrm{des}}+\frac{\rho}{2}\left(\bm{G}\bm{d}^{(i)}+\frac{\bm{\lambda}^{(i)}}{\rho}\right)\right),$
(12)
$\displaystyle\bm{a}^{(i+1)}=\frac{\left|2(1-\gamma)\bm{u}_{\mathrm{des}}+\rho\left(\bm{G}\bm{d}^{(i)}+\frac{\bm{\lambda}^{(i)}}{\rho}\right)\right|+2\gamma|\bm{u}_{\mathrm{des}}|}{\rho+2}.$
(13)

The update rule for $\bm{d}$ is obtained by solving (10) as

$\displaystyle\bm{d}^{(i+1)}$ $\displaystyle\ \
=\left(\bm{G}^{\mathsf{H}}\bm{G}+\frac{2\beta}{\rho}\bm{I}\right)^{-1}\bm{G}^{\mathsf{H}}\left(\bm{a}^{(i+1)}\odot\mathrm{e}^{{\mathrm{j}}\bm{\theta}^{(i+1)}}-\frac{\bm{\lambda}^{(i)}}{\rho}\right).$
(14)

By iteratively updating $\bm{\theta}^{(i)}$, $\bm{a}^{(i)}$, $\bm{d}^{(i)}$,
and $\bm{\lambda}^{(i)}$ by using (12), (13), (14), and (11), respectively,
starting with initial values, we can obtain the optimal driving signal
$\bm{d}$.

###  3.2 Time-Domain Filter Design

By using the proposed algorithm described in Sect. 3.1, we can obtain the
driving signal in the frequency domain. In practice, a finite impulse response
(FIR) filter is obtained by computing the inverse Fourier transform of
$\bm{d}$ for target frequency bins. However, if the driving signal is
independently determined for each frequency, it can have discontinuities
between frequency bins particularly for $\gamma=1$ because the phase of
$\bm{d}$ is arbitrary. These discontinuities can lead to an unnecessarily
large FIR filter length.

To overcome this issue, the differential-norm penalty, which is also used in
amplitude matching [15], can be applied. By introducing the subscript of the
index of the frequency bin $f\in\\{1,\ldots,F\\}$, we define the differential-
norm penalty for the $f$th frequency bin as

$\displaystyle D(\bm{d}_{f})=\|\bm{d}_{f}-\bm{d}_{f-1}\|_{2}^{2}.$ (15)

The optimization problem for the time-domain filter design is represented as

$\displaystyle\mathop{\rm
minimize}\limits_{\\{\bm{d}_{f}\\}_{f=1}^{F}}\sum_{f=1}^{F}\left[(1-\gamma)\left\|\bm{G}_{f}\bm{d}_{f}-\bm{u}_{{\mathrm{des}},f}\right\|_{2}^{2}\right.$
$\displaystyle\hskip
10.0pt\left.+\gamma\left\||\bm{G}_{f}\bm{d}_{f}|-|\bm{u}_{{\mathrm{des}},f}|\right\|_{2}^{2}\right]+\alpha\sum_{f=2}^{F}D(\bm{d}_{f})+\beta\sum_{f=1}^{F}\|\bm{d}_{f}\|_{2}^{2},$
(16)

where $\alpha$ is the weight for the differential-norm penalty term. ADMM is
similarly applied to solve (16), but its detailed derivation is omitted. The
derivation of ADMM for amplitude matching and a technique to reduce
computational complexity can be found in [15].

##  4 Experiments

We conducted experiments to evaluate the proposed method (Proposed) compared
with pressure matching (PM). First, numerical experimental results are shown
to evaluate the ILD of a listener and the amplitude response of the
synthesized sound field. Second, the results of listening experiments using a
practical system are presented.

###  4.1 Numerical Experiments

Figure 2: Experimental setup. Blue circles and yellow dots represent sources
and control points, respectively.

Numerical experiments in the frequency domain were conducted under the three-
dimensional free-field assumption. The target region $\Omega$ was a cuboid of
$1.0~{}\mathrm{m}\times 1.0~{}\mathrm{m}\times 0.04~{}\mathrm{m}$. As shown in
Fig. 2, 16 omnidirectional loudspeakers were placed along the borders of the
squares of $2.0~{}\mathrm{m}\times 2.0~{}\mathrm{m}$ at the heights of $z=\pm
0.1~{}\mathrm{m}$, and $24\times 24\times 2$ control points were regularly
placed inside $\Omega$. Therefore, the total number of loudspeakers and
control points were $L=32$ and $N=1152$, respectively. The desired sound field
was a spherical wave from the point source at $(2.0,0.0,0.0)~{}\mathrm{m}$. In
Proposed, $\gamma$ was set using (6) with $\omega_{\mathrm{T}}/2\pi=2000$ and
$\sigma=0.01$ so that the phase distribution becomes arbitrary above
$2000~{}\mathrm{Hz}$. The regularization parameter $\beta$ for Proposed and PM
was set as $\|\bm{G}^{\mathsf{H}}\bm{G}\|_{2}^{2}\times 10^{-3}$. The penalty
parameter $\rho$ in (8) was $1.0$.

We evaluated the ILDs of the synthesized sound field when a listener’s head
was in $\Omega$. The binaural signals at $\omega$ for the position
$\bm{r}_{\mathrm{H}}\in\Omega$ and azimuth direction
$\phi_{\mathrm{H}}\in[0,2\pi)$ of the listener’s head were denoted as
$b_{\mathrm{L}}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}},\omega)$ and
$b_{\mathrm{R}}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}},\omega)$ for the left
and right ears, respectively. $b_{\mathrm{L}}$ and $b_{\mathrm{R}}$ were
obtained by calculating the transfer function between the loudspeakers and the
listener’s ears using Mesh2HRTF [29, 30]. The ILD for $\bm{r}_{\mathrm{H}}$
and $\phi_{\mathrm{H}}$ was calculated in the frequency domain as

$\displaystyle\mathrm{ILD}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}})$
$\displaystyle:=10\log_{10}\frac{\sum_{\omega}|b_{\mathrm{L}}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}},\omega)|^{2}}{\sum_{\omega}|b_{\mathrm{R}}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}},\omega)|^{2}}.$
(17)

The evaluation measure was the normalized error ($\mathrm{NE}$) between the
synthesized and desired ILDs ($\mathrm{ILD}_{\mathrm{syn}}$ and
$\mathrm{ILD}_{\mathrm{true}}$, respectively) at $\bm{r}_{\mathrm{H}}$ defined
as

$\displaystyle\mathrm{NE}(\bm{r}_{\mathrm{H}})=\frac{\sum_{\phi_{\mathrm{H}}}|\mathrm{ILD}_{\mathrm{syn}}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}})-\mathrm{ILD}_{\mathrm{true}}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}})|}{\sum_{\phi_{\mathrm{H}}}|\mathrm{ILD}_{\mathrm{true}}(\bm{r}_{\mathrm{H}},\phi_{\mathrm{H}})|},$
(18)

where the summation for $\phi_{\mathrm{H}}$ was calculated for
$\\{0,\pi/32,\pi/16,\ldots,31\pi/32\\}~{}\mathrm{rad}$. The distributions of
$\mathrm{NE}(\bm{r}_{\mathrm{H}})$ on the $x$–$y$-plane at $z=0$ for Proposed
and PM are shown in Fig. 3. The evaluated positions were $5\times 5$ points on
the square of $1.0~{}\mathrm{m}\times 1.0~{}\mathrm{m}$. The $\mathrm{NE}$ of
Proposed was smaller than that of $\mathrm{PM}$ over the region. Note that the
ITDs were accurately synthesized in both methods below $2000~{}\mathrm{Hz}$.

Next, the amplitude response of the synthesized sound field was investigated.
In Fig. 4, the amplitude responses at the center of $\Omega$ of the desired
sound field and synthesized sound field of Proposed and PM are plotted. Owing
to the large variations in the amplitude response of PM above
$2000~{}\mathrm{Hz}$, the timbre of the source signal can be highly distorted.
In contrast, the almost flat amplitude response was achieved by Proposed at
the center of $\Omega$.

(a) PM (b) Proposed

Figure 3: Distribution of $\mathrm{NE}$ of ILDs. Figure 4: Amplitude responses
of desired and synthesized sound fields at the center of $\Omega$.

###  4.2 Listening Experiments

Listening experiments were conducted to evaluate the perceptual quality of the
synthesized sound field by using a practical loudspeaker array. The numbers
and positions of loudspeakers and control points were the same as those in
Fig. 2. The reference loudspeaker was placed at $(2.0,0.5,0.0)~{}\mathrm{m}$
as a primary sound source of the desired field. The driving signals of
Proposed and PM were obtained by assuming the loudspeakers as point sources.
The reverberation time $T_{60}$ of the room was around $0.19~{}\mathrm{s}$.

The perceptual quality was evaluated by using multiple stimuli with a hidden
reference and anchor (MUSHRA) [31]. Test participants were asked to rate the
difference between the reference and test signals of $10~{}\mathrm{s}$ on a
scale from 0 to 100. The reference and test signals are summarized as follows:

  * •

Reference: The source signal played back through the reference loudspeaker.

  * •

C1/Hidden anchor: The lowpass-filtered source signal up to
$3.5~{}\mathrm{kHz}$ played back through the reference loudspeaker.

  * •

C2/PM: The sound synthesized by PM and played back through the loudspeaker
array.

  * •

C3/Proposed: The sound synthesized by Proposed and played back through the
loudspeaker array.

  * •

C4/Hidden reference: The same as Reference.

The participants’ head center was approximately positioned at the center of
the target region by adjusting the chair, but they were able to rotate and
move their heads on the chair freely. The participants were able to listen to
each test signal repeatedly. Fourteen male subjects in their 20s and 30s were
included, and those who scored more than 60 on the hidden anchor, which was
one participant in this test, were excluded from the evaluation. Two source
signals, Vocals and Instrumental, taken from track 10 of MUSDB18-HQ [32] were
investigated.

Fig. 5 shows the box-and-whisker plots of the scores of each test signal. The
median score of Proposed was significantly higher than that of PM for both
Vocals and Instrumental. After the validation of the normality of data for C2
and C3 by the Shapiro–Wilk test, Welch’s $t$-test was conducted at a
significance level of $0.05$. The $p$ values for Vocals and Instrumental were
$9.1\times 10^{-4}$ and $2.0\times 10^{-3}$, respectively; therefore, there
were significant differences in mean scores between C2 and C3 in the both
cases.

(a) Vocals (b) Instrumental

Figure 5: Box-and-whisker plots of scores for each test signal. Circles
indicate outliers, and the notches represent 95% confidence intervals of the
median.

##  5 Conclusion

We proposed a sound field synthesis method based on the combination of
pressure and amplitude matching to improve perceptual quality. The cost
function is defined as square errors of pressure distribution for low
frequencies and amplitude distribution for high frequencies to alleviate the
effects of spatial aliasing artifacts. The ADMM-based algorithm to solve this
cost function is also derived. In the numerical experiments and listening
experiments, it was validated that the perceptual quality of the proposed
method can be improved from that of PM. Future work includes extended
listening experiments to evaluate perceptual quality in more detail.

##  6 ACKNOWLEDGMENT

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

## References

  * [1] P. A. Nelson, “Active control of acoustic fields and the reproduction of sound,” _J. Sound Vibr._ , vol. 177, no. 4, pp. 447–477, 1993.
  * [2] O. Kirkeby, P. A. Nelson, F. O. Bustamante, and H. Hamada, “Local sound field reproduction using digital signal processing,” _J. Acoust. Soc. Amer._ , vol. 100, no. 3, pp. 1584–1593, 1996.
  * [3] T. Betlehem and T. D. Abhayapala, “Theory and design of sound field reproduction in reverberant rooms,” _J. Acoust. Soc. Amer._ , vol. 117, no. 4, pp. 2100–2111, Apr. 2005.
  * [4] M. A. Poletti, “Three-dimensional surround sound systems based on spherical harmonics,” _J. Audio Eng. Soc._ , vol. 53, no. 11, pp. 1004–1025, Nov. 2005\.
  * [5] 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.
  * [6] 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.
  * [7] 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\.
  * [8] S. Spors, R. Rabenstein, and J. Ahrens, “The theory of wave field synthesis revisited,” in _Proc. 124th AES Conv._ , Amsterdam, Netherlands, Oct. 2008\.
  * [9] 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.
  * [10] J. Ahrens and S. Spors, “An analytical approach to sound field reproduction using circular and spherical loudspeaker distributions,” _Acta Acustica united with Acustica_ , vol. 94, no. 6, pp. 988–999, Nov. 2008.
  * [11] 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, Jan. 2009.
  * [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._ , vol. 21, no. 4, pp. 685–696, 2013.
  * [13] S. Koyama, K. Furuya, Y. Hiwasaki, Y. Haneda, and Y. Suzuki, “Wave field reconstruction filtering in cylindrical harmonic domain for with-height recording and reproduction,” _IEEE/ACM Trans. Audio, Speech, Lang. Process._ , vol. 22, no. 10, pp. 1546–1557, 2014.
  * [14] S. Koyama, T. Amakasu, N. Ueno, and H. Saruwatari, “Amplitude matching: Majorization-minimization algorithm for sound field control only with amplitude constraint,” in _Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)_ , Jun. 2021, pp. 411–415.
  * [15] 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.
  * [16] J. Blauert, _Spatial Hearing: The Psychophysicas of Human Sound Localization_. Cambridge: The MIT Press, 1997.
  * [17] A. Brughera, L. Dunai, and W. M. Hartmann, “Human interaural time difference thresholds for sine tones: The high-frequency limit,” _J. Acoust. Soc. Amer._ , pp. 2839–2855, 2013.
  * [18] S. Koyama and K. Arikawa, “Weighted pressure matching based on kernel interpolation for sound field reproduction,” in _Proc. Int. Congr. Acoust. (ICA)_ , Gyeongju, Republic of Korea, Oct. 2022.
  * [19] N. Ueno, S. Koyama, and H. Saruwatari, “Sound field recording using distributed microphones based on harmonic analysis of infinite order,” _IEEE Signal Process. Lett._ , vol. 25, no. 1, pp. 135–139, 2018.
  * [20] ——, “Directionally weighted wave field estimation exploiting prior information on source direction,” _IEEE Trans. Signal Process._ , vol. 69, pp. 2383–2395, 2021.
  * [21] E. G. Williams, _Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography_. London, UK: Academic Press, 1999.
  * [22] P. A. Martin, _Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles_. New York: Cambridge University Press, 2006.
  * [23] D. Colton and R. Kress, _Inverse Acoustic and Electromagnetic Scattering Theory_. New York, NY, USA: Springer, 2013\.
  * [24] S. Spors and R. Rabenstein, “Spatial aliasing artifacts produced by linear and circular loudspeaker arrays used for wave field synthesis,” in _Proc. 120th AES Conv._ , May 2006.
  * [25] 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._ , vol. 28, pp. 686–714, 2020.
  * [26] K. Kimura, S. Koyama, N. Ueno, and H. Saruwatari, “Mean-square-error-based secondary source placement in sound field synthesis with prior information on desired field,” in _Proc. IEEE Int. Workshop Appl. Signal Process. Audio Acoust. (WASPAA)_ , Oct. 2021, pp. 281–285.
  * [27] M. Poletti and T. Betlehem, “Creation of a single sound field for multiple listeners,” in _Proc. INTER-NOISE_ , 2014, pp. 2033–2040.
  * [28] C. Schoörkhuber, M. Zaynschirm, and R. H oldrich, “Binaural rendering of ambisonics signals via magnitude least squares,” in _Proc. DAGA_ , 2018, pp. 339–342.
  * [29] H. Ziegelwanger, W. Kreuzer, and P. Majdak, “Mesh2hrtf: Open-source software package for the numerical calculation of head-related transfer functions,” in _Proc. Int. Congr. Sound Vibr. (ICSV)_ , 2015.
  * [30] H. Ziegelwanger, P. Majdak, and W. Kreuzer, “Numerical calculation of listener-specific head-related transfer functions and sound localization: Microphone model and mesh discretization,” _J. Acoust. Soc. Amer._ , vol. 138, no. 1, pp. 208–222, 2015.
  * [31] International Telecommunication Union, Radiocommunication Sector, “Recommendation itu-r bs.1534-3: Method for the subjective assessment of intermediate quality level of audio systems,” 2014.
  * [32] Z. Rafii, A. Liutkus, F.-R. Stöter, S. I. Mimilakis, and R. Bittner, “MUSDB18-HQ - an uncompressed version of MUSDB18,” Dec. 2019. [Online]. Available: https://doi.org/10.5281/zenodo.3338373