# wav2shape: Hearing the shape of a drum machine

###### Abstract

Disentangling and recovering physical attributes, such as shape and material,
from a few waveform examples is a challenging inverse problem in audio signal
processing, with numerous applications in musical acoustics as well as
structural engineering. We propose to address this problem via a combination
of time–frequency analysis and supervised machine learning. We start by
synthesizing a dataset of sounds using the functional transformation method.
Then, we represent each percussive sound in terms of its time-invariant
scattering transform coefficients and formulate the parametric estimation of
the resonator as multidimensional regression with a deep convolutional neural
network. We interpolate scattering coefficients over the surface of the drum
as a surrogate for potentially missing data, and study the response of the
neural network to interpolated samples. Lastly, we resynthesize drum sounds
from scattering coefficients, therefore paving the way towards a deep
generative model of drum sounds whose latent variables are physically
interpretable.

##  1 Introduction

Throughout musical traditions, drums come in all shapes and sizes (see Figure
1). Such diversity in manufacturing results in a wide range of perceptual
attributes: bright, warm, mellow, and so forth. Yet, current approaches to
drum music transcription, which are based on one-versus-all classification,
fail to capture the multiple factors of variability underlying the timbre
perception of percussive sounds [1]. Instead, they regard each item in the
drum kit as a separate category, and rarely account for the effect of playing
technique. Therefore, in the context of music information retrieval (MIR), the
goal of broadening and refining the vocabulary of percussive sound recognition
systems requires to move away from discrete taxonomies.

Figure 1: Drums of various shapes and materials. Left to right: mbejn, 19th
century, Fang people of Gabon; ceramic drum, 1st century, Nasca people of
Peru; darabukka, 19th century, Syria; tympanum of a Pejeng-type drum, Bronze
age, Indonesia (Sumba); pakhavaj, 19th century, North India; ipu hula, 19th
century, Hawai’i; frame drum, 19th century, Native American people of Dakota;
Union army drum, ca. 1864, Pennsylvania. All images are in the public domain
and can be accessed at: %****␣han2020fa.tex␣Line␣75␣****www.metmuseum.org

In a different context, prior literature on musical acoustics has managed to
simulate the response of a drum from the knowledge of its shape and material.
Among studies on physical modeling of musical instruments, functional
transformation method (FTM) [2] and finite difference method (FDM) [3, 4] play
a central role. They rely on partial differential equations (PDE) to describe
the structural and material constraints imposed by the resonator. The
coefficients governing these equations may be varied continuously. Thus, PDE-
based models for drum sound synthesis offer a fine level of expressive control
while guaranteeing physical plausibility and interpretability.

From a musical standpoint, a major appeal behind physical models lies in
auditory perception: all other things being equal, larger drums tend to sound
lower, stiffer drums tend to sound brighter, and so forth. Yet, a major
drawback of PDE-based modeling for drum sound synthesis is that all shape and
material parameters must be known ahead of time. If, on the contrary, these
parameters are unknown, adjusting the synthesizer to match a predefined audio
sample incurs a process of multidimensional trial and error, which is tedious
and unscalable. This is unlike other methods for audio synthesis, such as
digital waveguide [5] or modal synthesis [6].

In this article, we strive towards resolving the tradeoff between control and
flexibility in drum sound synthesis. To this end, we formulate the
identification of percussive sounds as an inverse problem, thus combining
insights from physical modeling and statistical machine learning. Our main
contribution is wav2shape, i.e., a machine listening system which takes a drum
stroke recording as input and retrieves the shape parameters which produced
it. The methodological novelty of wav2shape lies in its hybrid architecture,
combining feature engineering and feature learning: indeed, it composes a 1-D
scattering transform and a deep convolutional network to learn the task of
shape regression in a supervised way. The advantage of choosing scattering
coefficient over conventional audio descriptor such as MFCC and CQT in
characterizing nonstationary sounds has been discussed in previous works [7,
8].

The subtitle of this paper is a deliberate reference to a famous mathematical
paper named “Can One Hear the Shape of a Drum?” [9]; that is, whether any two
isospectral planar domains are necessarily isometric. Since its publication,
this question has been answered affirmatively in the important particular
cases of circular and rectangular domains; but negatively in the general case,
with the construction of nonconvex counterexamples. Despite the evident
connection with our paper, we note that [9] and wav2shape strive towards
slightly different goals. First, while [9] makes no prior assumption on the
symmetries of the membrane, wav2shape focuses on representing rectangular and
circular membranes, which are by far the most common in music. In return,
while [9] is restricted to the recovery of the domain under forced
oscillations, wav2shape also expresses the effects of stiffness and damping,
both frequency-dependent and frequency-independent. These effects are crucial
for modeling the response of the drum membrane to a localized impulse, e.g.
induced by the player’s hand, a stick, or a mallet.

Our main finding is that, after training, wav2shape is able to generalize to
previously unseen shapes. Add an additional experiment, we interpolate the
value of scattering coefficients over the 2-D surface of the drum and verify
that the convnet in wav2shape generalizes to interpolated drum stroke
locations. Lastly, we invert the scattering transform operator, thus laying
the foundations for turning wav2shape into a deep generative model without
explicit knowledge of the partial differential equation (PDE) underlying the
vibration of the membrane.

##  2 Problem Statement

###  2.1 Multidimensional regression of PDE coefficients

The vibration of a drum obeys a partial differential equation (PDE), involving
both resonant and dissipative terms. In the following, we assume the
analytical form, boundary conditions, and initial conditions of this PDE to be
known—as Section 3 will discuss. Conversely, we take its vector of constant
coefficients $\boldsymbol{\theta}$ to be unknown. We represent the state of
the drum by the displacement field $\mathbf{X}_{\boldsymbol{\theta}}$ of its
membrane as a function of space $u\in[0,l]$ and time $t$. We place the origin
of the Cartesian coordinate system at the center of the drum
($u=u_{0}=(l/2,l/2)$) and the onset of the stroke ($t=0$). The goal of
wav2shape is to recover $\boldsymbol{\theta}$ from a single measurement of
$\mathbf{X}_{\boldsymbol{\theta}}$ near the origin.

###  2.2 Need for geometrical invariants

Let
$\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}:t\mapsto\mathbf{X}_{\boldsymbol{\theta}}\big{(}t,u=u_{0}\big{)}$
be the time series describing the displacement of the drum at its center. For
any given $\boldsymbol{\theta}$, the signal
$\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}$ lasts for about one second and
spans about $20\text{\,}\mathrm{kHz}$ in bandwidth. Therefore, once
discretized uniformly and truncated to a finite duration,
$\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}$ has a typical length of $10^{5}$
samples. Furthermore, Euclidean distances in the waveform domain are not
informative for recovering $\boldsymbol{\theta}$: for example, flipping the
polarity of the signal (i.e., from to
$\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}$ to
$-\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}$) produces a large Euclidean
distance, yet leaves $\boldsymbol{\theta}$ unchanged. More generally,
discrepancies in audio acquisition across samples, e.g. involving changes in
gain and DC bias, imply that the evolution of each
$\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}$ is only known up to a global
affine transformation. For this reason, a major challenge underlying the
development of wav2shape is to represent high-dimensional audio signals in a
feature space which satisfies certain geometrical invariants (such as
$\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}\mapsto
a\mbox{\boldmath{$x$}}_{\boldsymbol{\theta}}+b$) while preserving informative
variability in $\boldsymbol{\theta}$.

###  2.3 Need for phase demodulation

In addition to affine changes in the displacement domain,
$\mbox{\boldmath{$x$}}_{\theta}$ is also subject to random fluctuations in the
spatiotemporal domain. This is because, in practice, the origin ($u=l/2$,
$t=0$) of the Cartesian coordinate system is prone to small measurement
errors. Given that $\mathbf{X}_{\theta}$ oscillates rapidly in time and space
near the origin, such measurement errors incur large phase deviations. These
phase deviations affect Euclidean distances between waveforms. On the
contrary, long-range interactions between wave ridges are informative of modal
resonance and damping, regardless of phase. Hence, wav2shape must demodulate
fast oscillations in $\boldsymbol{x}_{\theta}$ in order to stably characterize
shape parameters $\boldsymbol{\theta}$.

###  2.4 Need for numerical stability to deformations

Let us denote by $\boldsymbol{\Phi}$ an instance of the wav2shape model. The
output of $\boldsymbol{\Phi}$ is a vector of constant coefficients to the PDE
governing the vibration of the drum:
$\boldsymbol{\tilde{\theta}}=\boldsymbol{\Phi}(\boldsymbol{x_{\theta}})$. We
evaluate wav2shape in terms of Euclidean distance between vectors describing
true vs. predicted drum shapes:

$L_{\Phi}(\boldsymbol{\theta})=\big{\|}\boldsymbol{\tilde{\theta}}-\boldsymbol{\theta}\big{\|}_{2}=\big{\|}\boldsymbol{\Phi}(\boldsymbol{x_{\theta}})-\boldsymbol{\theta}\big{\|}_{2}.$
(1)

This Euclidean distance is computed in a vector space of relatively high
dimension—in this article, we encode $\boldsymbol{\theta}$ in dimension five.
Thus, the supervised prediction of $\boldsymbol{\theta}$ is exposed to the
curse of dimensionality. In order to learn the wav2shape function
$\boldsymbol{\Phi}$ from limited annotated data, it is necessary to map
waveform samples $\boldsymbol{x}_{\theta}$ to a feature space in which
coordinate-wise variations of $\boldsymbol{\theta}$ are disentangled and
linearized.

In the context of wav2shape, some factors of variability in
$\boldsymbol{\theta}$ (e.g., pitch) are most intuitive in the frequency
domain, while others (e.g., rate of damping) are most intuitive in the time
domain. Therefore, it is advantageous to train a machine learning system to
regress $\boldsymbol{\theta}$ in the time–frequency domain, rather than the
time domain. Section 4 will present how wav2shape combines a scattering
transform and a deep convolutional neural network, as an unsupervised feature
extraction stage and a supervised nonlinear regression stage respectively.

##  3 From shape to wave:
physical synthesis model

###  3.1 Formulation as a fourth-order PDE

Let us recall the wave equation in dimension two:

$\dfrac{\partial^{2}\mathbf{X}}{\partial
t^{2}}-c^{2}\nabla^{2}\mathbf{X}=\mathbf{Y}(t,u),$ (2)

where $c$ is the speed of sound over the drum membrane; the symbol
$\nabla^{2}$ denotes the spatial Laplacian operator; and the scalar field
$\mathbf{Y}$ represents the gesture of the musician. Throughout this article,
we assume the spatiotemporal field $\mathbf{Y}$ to be factorizable into a
temporal component $\boldsymbol{y}_{\mathrm{t}}$ and a spatial component
$\boldsymbol{y}_{\mathrm{u}}$.

Although the formulation above may be sufficient to identify stationary
eigenmodes in $\mathbf{X}$, it does not faithfully characterize the response
of a drum membrane to a percussive excitation $\mathbf{Y}$ [10]. To address
this issue, we consider the stiffness $S$ of the drum membrane as a function
of its Young’s modulus and its moment of inertia. Furthermore, air drag
induces an energy dissipation in $\mathbf{X}$ through a first-order damping
coefficient $d_{1}$. Lastly, near the boundary of the drum, the mechanical
coupling between the membrane and the body of the drum also causes energy
dissipation through a third-order damping coefficient $d_{3}$.

Once the terms $S$ (stiffness), $d_{1}$ (first-order damping), and $d_{3}$
(third-order damping) have been taken into account, the PDE governing the
displacement field $\mathbf{X}$ becomes:

$\displaystyle\left(\dfrac{\partial^{2}\mathbf{X}}{\partial
t^{2}}(t,u)-c^{2}\nabla^{2}\mathbf{X}(t,u)\right)$
$\displaystyle+S^{4}\big{(}\nabla^{4}\mathbf{X}(t,u)\big{)}+\dfrac{\partial}{\partial
t}\Big{(}d_{1}\mathbf{X}(t,u)+d_{3}\nabla^{2}\mathbf{X}(t,u)\Big{)}$
$\displaystyle=\mathbf{Y}(t,u)=\boldsymbol{y}_{\mathrm{t}}(t)\boldsymbol{y}_{\mathrm{u}}(u),$
(3)

where the spatiotemporal field $\nabla^{4}\mathbf{X}$ denotes the “double
Laplacian” of $\mathbf{X}$, i.e., the Laplacian of $\nabla^{2}\mathbf{X}$.

###  3.2 Boundary conditions

For the sake of simplicity and conciseness, we only address the case of a
rectangular membrane, e.g., that of a cajón. The important case of a circular
membrane (see Figure 1) could be derived from Equation 3 with the same tools
as presented hereafter; yet, it would incur a conversion to polar coordinates,
and the resort to Bessel functions. We direct readers to [11] for the
important case of the circular membrane. Note, in this case, that the transfer
function method (TFM) is an alternative denomination for the functional
transformation method (FTM).

We consider the membrane to be a rectangle of width $l_{1}$, length $l_{2}$,
and aspect ratio $\alpha=l_{1}/l_{2}$. Along the edges of this rectangle, we
assume the displacement field to be null: for every $t$, $\mathbf{X}(t,u)=0$
if $u_{1}=0$, $u_{1}=l_{1}$, $u_{2}=0$, or $u_{2}=l_{2}$. This is tantamount
to assuming that the shape of the drum remains fixed throughout the duration
of the percussive stroke.

###  3.3 Functional transformation method (FTM)

The Laplace transform of $\mathbf{X}$ over the time dimension is

$\mathcal{L}\\{\mathbf{X}\\}(s,u)=\int_{0}^{+\infty}\mathbf{X}(t,u)\,\exp(-st)\;\mathrm{d}t,$
(4)

In the Laplace domain, Equation 3 becomes

$\displaystyle S^{4}\left(\nabla^{4}\mathbf{\mathcal{L}\\{X\\}}(s,u)\right)$
$\displaystyle+$
$\displaystyle\left(sd_{3}-c^{2}\right)\nabla^{2}\mathbf{\mathcal{L}\\{X\\}}(s,u)$
$\displaystyle+$
$\displaystyle\left(s^{2}+sd_{1}\right)\mathbf{\mathcal{L}\\{X\\}}(s,u)$
$\displaystyle=\mathcal{L}\\{\boldsymbol{y}_{\mathrm{t}}\\}(s)\boldsymbol{y}_{\mathrm{u}}(u).$
(5)

The interest of the Laplace domain is that, in comparison with Equation 3, the
equation above replaces temporal derivatives with simpler algebraic terms.
Similarly, spatial derivatives may be eliminated by means of the Sturm-
Liouville transformation (SLT), as detailed in [12, 13, 14]. Once in the
Laplace-SLT domain, the solution of the PDE can be recovered in the
spatiotemporal domain by performing an inverse Sturm-Liouville and inverse
Laplace transform consecutively. In this context, drums with a rectangular
membrane are conceptually simpler: indeed, the inverse Sturm-Liouville
transformation boils down to a Fourier series decomposition [12]. Thus, in the
particular case described in Section 3.2, we may skip the SLT altogether and,
instead, decompose the Laplace domain solution as a Fourier series over the
2-D variable $u$.

At any fixed $s\in\mathbb{C}$, the spatial field
$u\mapsto\mathcal{L}\\{\mathbf{X}\\}(s,u)$ is absolutely continuous. We index
each mode in $\mathcal{L}\\{\mathbf{X}\\}$ by the pair
$m=(m_{1},m_{2})\in\mathbb{Z}^{2}$, and denote by
$\mathcal{\widehat{L}}_{m}(\mathbf{X})(s)\in\mathbb{C}$ the associated Fourier
coefficients:

$\displaystyle\mathcal{L}\\{\mathbf{X}\\}(s,u)$
$\displaystyle=\sum_{m\in\mathbb{N}^{2}}\mathcal{\widehat{L}}_{m}\\{\mathbf{X}\\}(s)\sin\left(\dfrac{m_{1}\pi
u_{1}}{l_{1}}\right)\sin\left(\dfrac{m_{2}\pi u_{2}}{l_{2}}\right)$ (6)

Similarly, we decompose $\boldsymbol{y}^{\mathrm{u}}$ into a series of 2-D
Fourier coefficients $\widehat{y}_{m}^{\mathrm{u}}$. Plugging the equation
above into Equation 5 allows a modal identification of the form:

$\displaystyle\mathcal{\widehat{L}}_{m}\\{\mathbf{X}\\}(s)$
$\displaystyle=\mathcal{L}\\{\boldsymbol{h}_{m}\\}(s)\times\mathcal{L}\\{\boldsymbol{y}^{\mathrm{t}}\\}(s)\times\widehat{y}_{m}^{\mathrm{u}}$
$\displaystyle=\frac{\mathcal{L}\\{\boldsymbol{y}^{\mathrm{t}}\\}(s)\times\widehat{y}_{m}^{\mathrm{u}}}{(s-z_{m})(s-\overline{z_{m}})},$
(7)

where the complex numbers $z_{m}$ and their conjugates $\overline{z_{m}}$
denote the poles of resonance of the impulse response $\boldsymbol{h}_{m}$.

After defining the constant $\gamma_{m}=m_{1}^{2}+m_{2}^{2}/\alpha^{2}$, a
straightforward computation leads to

$\mathfrak{R}(z_{m})=\dfrac{d_{3}\gamma_{m}-d_{1}}{2}$ (8)

for the real part, and

$\mathfrak{I}(z_{m})^{2}=\left(S^{2}-\dfrac{d_{3}^{2}}{4}\right)\gamma_{m}^{2}+\left(c^{2}+\dfrac{d_{1}d_{3}}{2}\right)\gamma_{m}-\dfrac{d_{1}^{2}}{4}$
(9)

for the squared imaginary part. Each impulse response $\boldsymbol{h}_{m}$ is
a real-valued sine wave with an exponential decay:

$\boldsymbol{h}_{m}(t)=\exp\big{(}\mathfrak{R}(z_{m})t\big{)}\sin\big{(}\mathfrak{I}(z_{m})t\big{)}.$
(10)

Lastly, an inverse Laplace transform of every term in Equation 6 yields the
following closed-form expression for $\mathbf{X}$:

$\displaystyle\mathbf{X}(t,u)=$
$\displaystyle\sum_{m\in\mathbb{N}^{2}}\big{(}\boldsymbol{y}^{\mathrm{t}}\ast\boldsymbol{h}_{m})(t)$
$\displaystyle\times\widehat{y}_{m}^{\mathrm{u}}\sin\left(\dfrac{m_{1}\pi
u_{1}}{l_{1}}\right)\sin\left(\dfrac{m_{2}\pi u_{2}}{l_{2}}\right),$ (11)

where the asterisk denotes the convolution operator.

###  3.4 Reparametrization

Although the tuplet $(S,c,d_{1},d_{3},\alpha)$ suffices to describe the
physical system in Equation 3, this tuplet remains unwieldy from a computer
music standpoint. Indeed, software plugins for drum sound synthesis usually
have knobs for “pitch” and “duration”; yet, these two perceptual attributes do
not appear clearly in Equation 3. Therefore, we map the tuplet above to a 5-D
space in which pitch and duration may be controlled intuitively.

Given a mode $\boldsymbol{h}_{m}$ (see Equation 10), we denote its carrier
frequency by the imaginary part $\omega_{m}=\mathfrak{I}(z_{m})$ and its
modulation frequency by the negative real part
$\sigma_{m}=-\mathfrak{R}(z_{m})$. The fundamental frequency of
$\boldsymbol{h}_{m}$ is perceived as proportional to $\omega_{m}$ while its
duration is perceived as inversely proportional to $\sigma_{m}$. By
convention, we take the mode of largest spatial extent as the reference for
the fundamental frequency and duration of $\mathbf{X}$. Setting $m=(1,0)$ in
Equation 9 yields the fundamental frequency:

$\omega=\omega_{(1,0)}=\sqrt{\frac{\beta^{4}}{\alpha^{2}}S^{4}+\dfrac{\beta}{\alpha}c^{2}-\dfrac{1}{4}\left(\dfrac{\beta}{\alpha}d_{3}-d_{1}\right)^{2}},$
(12)

where the dimensionless constant $\beta=\alpha+1/\alpha$ is associated to the
aspect ratio $\alpha$ of the rectangular drum membrane (see Section 3.2).
Finally, we define the duration of $\mathbf{X}$ as the inverse of the
modulation frequency of the mode $\boldsymbol{h}_{(1,0)}$. Equation 8 becomes:

$\tau=\dfrac{1}{\sigma_{(1,0)}}=\dfrac{2}{d_{1}-\dfrac{\beta}{\alpha}d_{3}}.$
(13)

Furthermore, we define the frequency-dependent damping of $\mathbf{X}$ as

$p=\dfrac{d_{3}}{\beta d_{3}-\alpha d_{1}}$ (14)

and its dispersion as

$D=\dfrac{1}{\alpha\omega}\sqrt{S^{4}-\dfrac{d_{3}^{2}}{4}}.$ (15)

We describe the “shape” of the drum as the 5-D vector
$\boldsymbol{\theta}=(\omega,\tau,p,D,\alpha)$. Once defined the value of
$\boldsymbol{\theta}$, we iterate over the multiindex
$m=(m_{1},m_{2})\in\mathbb{N}^{2}$, set
$\gamma_{m}=m_{1}^{2}+m_{2}^{2}/\alpha^{2}$, and define the associated
modulation frequency

$\sigma_{m}=\frac{1+p(\gamma_{m}-1)}{\tau}$ (16)

and squared carrier frequency

$\displaystyle\omega_{m}^{2}$ $\displaystyle=D^{2}\omega^{2}\gamma_{m}^{2}$
$\displaystyle+\left(\dfrac{(1-p)^{2}}{\tau^{2}}+\omega^{2}(1-D^{2})\right)\gamma_{m}$
$\displaystyle-\dfrac{(1-p)^{2}}{\tau^{2}}.$ (17)

Then we define the exponentially modulated sinusoid
$\boldsymbol{h}_{m}:t\mapsto\exp(-\sigma_{m}t)\sin(\omega_{m}t)$ as in
Equation 10. The infinite series $(\boldsymbol{h}_{m})$ fully describes the
response of the drum to an arbitrary excitation $\mathbf{Y}$ (see Equation 3).
In practice, we compute impulse responses $(\boldsymbol{h}_{m})$ over a finite
grid of $M^{2}=100$ modes, i.e., ten modes in each dimension.

Observe that the parameter $\tau$ affects only modulation frequencies
$\sigma_{m}$ without affecting carrier frequencies $\omega_{m}$. Conversely,
the parameters $\omega$ and $D$ only affect carrier frequencies $\omega_{m}$
without affecting modulation frequencies $\sigma_{m}$. As regards $p$ and
$\alpha$, they affect both the carrier frequency and the modulation frequency
of every mode.

###  3.5 Real-time implementation as a VST plugin

We implement the physical model above in the C++ language by means of the JUCE
application framework111Link to download source code and executable binaries:
https://github.com/lylyhan/Thesis. As a result, our drum sound synthesizer is
portable on Windows, Mac OS X, and Linux. Furthermore, we package our software
in the VST (Virtual Studio Technology) format. Thus, it can be integrated into
a digital audio workstation (DAW) such as Ableton Live, Adobe Audition,
Audacity, Cubase, Logic Pro, Max/MSP, or REAPER.

Figure 2: Screenshot of the graphical user interface (GUI) of our real-time
drum synthesis plugin.

The plugin contains three physical models: 1D string model, 2D rectangular
drum model, and 3D cuboid model. All models solve the PDE in Equation 3 via
the functional transformation method (FTM), adapted to different dimensions of
the spatial variable $u$.

Figure 2 illustrates the appearance of the plugin, in its “drum” (i.e., 2-D)
variant. The plugin contains four knobs, mapped as follows: “sustain” $\tau$,
“roundness” $p$, inharmonicity $D$, and “squareness” (i.e., aspect ratio)
$\alpha$. Moreover, we map the logarithm of fundamental frequency $\omega$ to
MIDI note numbers. All knobs can be controlled by dragging the pointer on the
screen or via a MIDI controller.

In practice, we set the excitation field $\mathbf{Y}$ to be a 2-D Gaussian of
fixed width. The user can set the center $u_{0}$ of this Gaussian by clicking
inside a depiction of the drum membrane on the screen. Our real-time
implementation discretizes the time variable $t$ according to the chosen
sampling rate, and renders the spatiotemporal field
$\mathbf{X}_{\boldsymbol{\theta}}(t,u)$ by summing up modes in Equation in 10.
Note that, as long as the shape vector $\boldsymbol{\theta}$ does not vary
between drum hits, the values of $z_{m}$ and $\boldsymbol{h}_{m}$ can be
precomputed without knowledge of $u_{0}$.

##  4 From wave to shape:
Machine listening model

Our problem statement (Section 2) stressed the importance of geometrical
invariants, phase demodulation, and numerical stability to deformations in the
context of regressing shape ($\boldsymbol{\theta}$) from wave
($\boldsymbol{x}_{\boldsymbol{\theta}}$). In this section, we present the
“wav2shape” machine listening model and explain how it satisfies these
mathematical properties. This model has a hybrid architecture: it composes a
feature engineering stage (1-D scattering transform) and a feature learning
stage (deep convolutional network) in a supervised way.

###  4.1 Scattering transform

Let $\boldsymbol{\psi}\in\mathbf{L}^{2}(\mathbb{R},\mathbb{C})$ a Hilbert-
analytic filter with null average, unit center frequency, and quality factor
$Q$ equal to one. We define a wavelet filterbank as the family
$\boldsymbol{\psi}_{j}:t\mapsto 2^{-j}\boldsymbol{\psi}(2^{-j}t)$ for integer
$j$. Each wavelet $\boldsymbol{\psi}_{j}$ has a center frequency proportional
to $2^{-j}$ and an effective receptive field proportional to $2^{j}$ in the
time domain.

We define the scalogram of $y$ as the complex modulus of its discrete wavelet
transform (DWT):

$\mathbf{U_{1}}\boldsymbol{x}:(t,j_{1})\longmapsto\left|\int_{-\infty}^{+\infty}\mbox{\boldmath{$x$}}(t^{\prime})\mbox{\boldmath{$\psi_{j_{1}}$}}(t-t^{\prime})\;\mathrm{d}t^{\prime}\right|.$
(18)

Likewise, we define a second layer of nonlinear transformation for $y$ as the
“scalogram of its scalogram”:

$\mathbf{U_{2}}\boldsymbol{x}:(t,j_{1},j_{2})\longmapsto\Big{|}\big{|}\boldsymbol{x}\ast\mbox{\boldmath{$\psi_{j_{1}}$}}\big{|}\ast\mbox{\boldmath{$\psi_{j_{2}}$}}\Big{|}(t),$
(19)

where the asterisk denotes a convolution product.

Every layer in a scattering network composes an invariant linear system
(namely, the complex DWT) and a pointwise operation (the complex modulus).
Thus, by recurrence over the depth variable $n$, every tensor
$\mathbf{U_{n}}\mbox{\boldmath{$y$}}$ is equivariant to the action of delay
operators. This alternation of convolution and modulus transform provides
complementary high-frequency wavelet coefficients [15].

In order to replace this equivariance property by an invariance property, we
integrate each $\mathbf{U_{n}}$ over some predefined time scale $T=2^{J}$,
yielding the invariant scattering transform:

$\displaystyle\mathbf{S_{n}}\mbox{\boldmath{$x$}}:(t,p)\longmapsto\int_{-\infty}^{+\infty}\mathbf{U_{n}}(t^{\prime},p)\mbox{\boldmath{$\phi$}}_{T}(t-t^{\prime})\;\mathrm{d}t^{\prime}$
(20)

where the $n$-tuple $p=(j_{1}\ldots j_{n})$ is known as a scattering path and
the function $\mbox{\boldmath{$\phi$}}_{T}$ is a real-valued low-pass filter
of time scale $T$. The number of layers is referred to as the order of the
scattering transform. Finally, we concatenate invariant scattering transform
coefficients of different orders:

$\displaystyle\mathbf{S}\mbox{\boldmath{$x$}}(t,p)=[\mathbf{S_{0}}\mbox{\boldmath{$x$}}(t),\mathbf{S_{1}}\mbox{\boldmath{$x$}}(t),...,\mathbf{S_{N}}\mbox{\boldmath{$x$}}(t)](p),$
(21)

where the path $p$ is a multiindex tuple containing between zero and $N$
entries. We direct readers to [16] for further mathematical details on the
scattering transform.

The two most important hyperparameters of the scattering transform are its
scale $J$ and its order $N$. A higher scale determines the window size,
reduces the number of time bins, and produces more scattering coefficients for
each time bin. Higher-order scattering coefficients encode and layer energy
extracted from the maximum to a number of shorter time scales, and thus
introduce a “deep”, nonlinear characterization of spectrotemporal modulations.

In this article, we set $J=8$ and $N=2$ unless specified otherwise. We compute
the scattering transform by means of the Kymatio library [17], using PyTorch
as a backend222Official website of Kymatio library: https://kymat.io.

###  4.2 Deep convolutional network: wav2shape

In order to learn a nonlinear mapping between waveform and the set of physical
parameters, we train a convolutional neural network, dubbed wav2shape (“wave
to shape”). Comprising four 1-D convolutional layers and two fully connected
dense layers, wav2shape is configured as follows:

  * •

layer 1: The input feature matrix passes through a batch normalization layer,
then 16 convolutional filters with a receptive field of 8 temporal samples.
The convolution is followed by a rectified linear unit (ReLU) and average
pooling over 4 temporal samples.

  * •

layer 2, 3, and 4: same as layer 1, except that the batch normalization
happens after the convolution. The average pooling filter in layer 4 has a
receptive field of 2 temporal samples, due to constraint in the time
dimension. After that, layer 4 is followed by a “flattening” operation.

  * •

layer 5: 64 hidden units, followed by a ReLU activation function.

  * •

layer 6: 5 hidden units, followed by a linear activation function.

Instead of supplying “raw” scattering coefficients to the first layer of
wav2shape, we apply a logarithmic transformation of the form

$\rho\big{(}{\mathbf{S}}\boldsymbol{x}\big{)}(t,p)=\log\left(1+\dfrac{\mathbf{S}\boldsymbol{x}(t,p)}{\varepsilon}\right),$
(22)

which has the effect of empirically Gaussianizing the statistical distribution
of each coefficient [7].

We set $\varepsilon=10^{-3}$ after verifying informally that this value yields
features which are similar enough for slightly perturbed audio signals with
imperceptible difference, yet still sufficiently distinct across different
drum shapes $\boldsymbol{\theta}$. Smaller values of the hyperparameter
$\varepsilon$ yields more discriminating feature representations, however too
small an $\varepsilon$ might magnify the otherwise imperceptible difference
between audio signals in feature space.

During training, we minimize mean squared error between the ground truth and
predicted $\boldsymbol{\theta}$ using the Adam optimizer. We use a minibatch
size of 64 and train for 30 epochs with 50 steps per epoch, i.e. 96k samples
in total. The validation set accuracy is checkpointed after each epoch to
identify the best performing model.

##  5 Experiments

###  5.1 Dataset

We synthesize a dataset of percussive sounds by discretizing the physical
parameters

$\boldsymbol{\theta}=\\{\omega,\tau,\log p,\log D,\alpha\\}$ (23)

uniformly, thus resulting in a five-dimensional hypercube. Each sound is
computed with the same temporal excitation; that is, a Dirac impulse in the
time domain ($\boldsymbol{y}^{\mathrm{t}}=\delta_{t}$) and a Gaussian in the
spatial domain:

$\hat{y}_{m}^{u}=\mathit{G}(\mu=l/2,\sigma=0.4),$ (24)

peaking at the center of the drum. Each drum sound lasts for $2^{15}$ audio
samples, i.e., about $1.5$ second at a sample rate of $22050$ Hz.

Along each dimension of the five-dimensional hypercube of shape parameters
$\boldsymbol{\theta}$, we curate the validation set by carving out the
“center” $60\%$ range over every dimension. This results in $0.6^{5}\approx
7.8\%$ of the total 100k samples being assigned to validation set.
Furthermore, we assign the surrounding sample space proportionally to training
set and test set. This ensures that the training set, validation set, and test
set do not overlap. There are 82221 samples in training set, 10k samples in
test set, and 7779 samples in validation set.

###  5.2 Shape Regression

Figure 3: Training the convolutional neural network with different choices of
scale $J$, order $N$ and scaling factor $\varepsilon$ as input scattering
coefficients yields varying learning robustness. The upper 2 diagrams show
effects on validation regression loss by selecting $J\in\\{6,8,10,12,14\\}$
and orders $N\in\\{1,2\\}$, where the upper and lower diagrams are using
scaling factor of $\varepsilon\in\\{10^{-3},10^{-1}\\}$ respectively during
feature preprocessing. With $\varepsilon=10^{-3}$, the best-performing model
is $J=8$, order $N=2$. The bottom diagram demonstrates distribution of the
absolute regression error of each individual physical parameter when applying
the best model on test set. Box and whisker edges denote quartiles and deciles
respectively.

The best performing wav2shape model results from a trial-and-error process of
hyperparameter optimization. We perform ten trials of training with different
values of scale $J$ and order $N$ to find the most successful input feature.
Figure 3 (top) summarizes our findings. Note that the input feature dimension
varies with $J$: thus, when the resulting time dimension is small,the number
of average-pooling filters and their receptive field sizes need to be changed
accordingly. Apart from the case of $J=8$ detailed in 4.2, $J=6$ uses four
average-pooling of receptive field 4; $J=10$ uses two of receptive fields 4
and one of receptive field 2; $J=12$ uses one each of receptive field 4 and 2;
and $J=14$ uses only one of receptive field 2.

We evaluate wav2shape in terms of Euclidean distance between prediction and
the normalized ground truth $\theta$. As a point of comparison, the mean
Euclidean distance between two points drawn uniformly at random in a
5-dimensional hypercube of unit hypervolume is around 0.87. In all of the
models, the minimum validation loss is far below this value: this indicates
that all variations of wav2shape generalize beyond the training set. The best
performing model is achieved with $J=8$ second-order scattering coefficients
scaled by $\varepsilon=10^{-3}$ as input, where the lowest minimum validation
loss across ten trials is around 0.0129.

Consistently with previous publications on the scattering transform, we
observe that, for all values of $J$, shape regression with $N=2$ outperforms
$N=1$. Indeed, the double nonlinearity in second-order scattering transform
contributes to the demodulation of nonstationarities, such as those found at
the onset of a drum sound. Meanwhile, larger scale $J$ increases the maximum
time window size, thus encodes the audio signal
$\boldsymbol{x}_{\boldsymbol{\theta}}$ with a lower sample rate yet more
coefficients along the frequency dimension. As $J$ increases, the audio
descriptor is more stable to deformations yet less discriminative to
variations in drum shape.

Figure 3 (bottom) breaks down the regression error of our best performing
model according to different dimensions of the shape parameter
$\boldsymbol{\theta}$. We observe that our model is the most accurate on
parameters $\tau$ and $\omega$ while being the least accurate on parameters
$p$ and $D$. An explanation is that $\tau$ and $\omega$ are the two parameters
which more directly affect poles of the system.

On the other hand, both $D$ and $p$ have asymptotic influences on the poles as
the mode number increases. Specifically $\tau_{m}\approx\tau_{1}/(pm^{2})$ and
$\omega_{m}\approx D\omega_{1}m^{2}$ for large m. These imply that effects of
changing $p$ and $D$ would be more obvious when sound is synthesized with more
modes. In our dataset each sound is summed only up to mode 10 due to time
constraints. Thus this deficiency in higher modal data might have also caused
the result.

###  5.3 Hearing shapes from neighboring sounds

Figure 4: Diagram of the closed-loop system: we use physical sound synthesis,
scattering transform, gradient descent, and the trained deep neural network
wav2shape to traverse between physical parameters, sound and scattering
domain. DNN: Deep Neural Network. GD: Gradient Descent.

To examine the stability of scattering transform, we construct a closed-loop
system that allows us to traverse between sound, physical and scattering
domains. Figure 4 shows a diagram of this closed-loop system.

We begin by selecting a drum shape, i.e. some random combination of physical
parameters $\boldsymbol{\theta}$. Secondly, we interpolate scattering
transform coefficients $\mathbf{SX}(t,u_{1},u_{2})$ on the drum. Specifically,
we compute the scattering transform at its neighboring points:
$\mathbf{SX}(t,u_{1}-\delta,u_{2})$, $\mathbf{SX}(t,u_{1}+\delta,u_{2})$,
$\mathbf{SX}(t,u_{1},u_{2}-\delta)$, $\mathbf{SX}(t,u_{1},u_{2}+\delta)$ and
approximate the scattering coefficients at $(u_{1},u_{2})$ from those of its
neighbors. Thirdly, we regress physical parameters from
$\widetilde{\mathbf{S}}\mbox{\boldmath{$x$}}(u_{1},u_{2})$ via the wav2shape
model, yielding a vector $\boldsymbol{\theta}^{\star}$. Lastly, we measure the
mean squared error between the predicted shape $\boldsymbol{\theta}^{\star}$
and the true shape $\boldsymbol{\theta}$.

The motivation behind this interpolation procedure is two-fold. First, we
examine the ability of the scattering transform to linearize the dependency of
the drum signal $\boldsymbol{x}_{\theta}$ with respect to the location of the
stroke. Secondly, we inquire whether wav2shape, which is trained on signals
measured at the exact center of the drum, remains capable of predicting the
shape from surrounding measurements.

In order to approximate the scattering coefficients at $(u_{1},u_{2})$, we
apply a four-point linear interpolation, i.e., an unweighted average of
neighboring coefficients along the four cardinal directions: North, East,
South, and West. We measure the approximation error of each scattering path
$p$ in terms of its discretized Laplacian

$\displaystyle\nabla^{2}\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u_{1},u_{2},p)$
$\displaystyle=\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u_{1},u_{2},p)$
$\displaystyle-\frac{1}{4}\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u_{1}-\delta,u_{2},p)$
$\displaystyle-\frac{1}{4}\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u_{1}+\delta,u_{2},p)$
$\displaystyle-\frac{1}{4}\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u_{1},u_{2}-\delta,p)$
$\displaystyle-\frac{1}{4}\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u_{1},u_{2}+\delta,p),$
(25)

where the step size $\delta$ is equal to $10\%$ of the side length of the
drum. For a given scattering path $p$ and time instant $t$, the equation above
measures the curvature of the manifold associated to
$u\mapsto\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u,p)$. If this manifold
is approximately flat, the linear interpolation is relatively accurate and the
discretized Laplacian is relatively small.

We summarize the discretized Laplacian above by taking its $\ell^{2}$ norm
over time and across scattering paths:

$\displaystyle\mathbf{H}\mathbf{X_{\boldsymbol{\theta}}}($ $\displaystyle
u_{1},u_{2})=$
$\displaystyle\sqrt{\int_{\mathbb{R}}\sum_{p}\nabla^{2}\mathbf{S}\mathbf{X}_{\boldsymbol{\theta}}(t,u_{1},u_{2},p)^{2}\;\mathrm{d}t,}$
(26)

thus yielding a matrix which is indexed by the spatial variable
$u=(u_{1},u_{2})$.

As an illustration, Figure 5 (left) shows the matrix
$\mathbf{H}\mathbf{X_{\boldsymbol{\theta}}}$ as a heatmap, for a fixed value
of the vector $\boldsymbol{\theta}$. As a point of comparison, we also compute
a Laplacian heatmap for Fourier modulus coefficients (Figure 5, right). We
observe a symmetric pattern over the surface of the drum. The darker regions
of this pattern correspond to the locations on the drum in which the
approximation of scattering coefficients by means of linear approximation is
the least valid. Interestingly, the locations of best fit do not lie near the
center, but between the four axes of symmetry of the drum.

By application of the Parseval theorem, the scattering transform and the
Fourier transform have comparable $\ell^{2}$ norms, i.e., the norm of the
signal $\boldsymbol{x}_{\boldsymbol{\theta}}$ [16]. Therefore, the heatmaps in
Figure 5 can be compared with the same numerical graduations. Over the surface
of the drum, we observe that the Laplacian of the scattering transform has a
smaller $\ell^{2}$ norm than the Laplacian of the Fourier transform modulus.
This difference reflects the better ability of the scattering transform to
linearize the dependency of the signal $\boldsymbol{x}_{\boldsymbol{\theta}}$
with respect to the origin $u$ of the excitation.

Figure 5: Heatmaps of Laplacian of scattering coefficients (left) and Fourier
coefficents (right) on the drum. The upper pair uses Dirac impulses to model
the excitation in space and the lower pair uses a Gaussian-shaped spatial
envelope. Darker colors reflect a larger $\ell^{2}$ norm of the Laplacian.

As an additional experiment, we apply the wav2shape model to interpolated
scattering coefficients. We sample the shape vector $\theta$ from three
distinct distributions: the validation set (7779 samples), the test set (10k
samples with same distribution as training set), and a previously unseen test
set (10k samples) drawn uniformly at random.

Figure 6 summarizes our results. We find that wav2shape is capable of
recovering the shape vector $\boldsymbol{\theta}$ with a relative mean squared
error around $0.15$. In comparison, a random guess would yield a relative mean
squared error of the order of $0.87$ (see Section 5.2). However, the error of
wav2shape on interpolated scattering coefficients is larger than the error on
true scattering coefficients, i.e., 0.0129 on the validation set. Such
discrepancy in shape regression accuracy results from interpolation error,
manifested by the nonzero Laplacian at $u=l/2$ on the drum (see Figure 5).
Future work will investigate methods to improve the ability of wav2shape to
generalize to off-center stroke locations.

Figure 6: Comparison of prediction error when the best performing wav2shape
model is applied onto scattering coefficients that are synthesized versus
interpolated at $u=0$ on the same drum. 60 drums are randomly selected from
three distributions: validation set (unseen by the model), test set (same
distribution as training set) and random range (unseen by the model). The red
line at 0.87 indicates regression loss achieved by a uniform random guess.

###  5.4 Reconstruction

Figure 7: Top plot: spectrogram of a synthetic drum sound. Bottom six plots:
spectrograms of reconstructed sounds from order 1 and order 2 scattering
coefficients, with different valus of the maximum scattering scale:
$J\in\\{3,5,7\\}$. Brighter colors denote greater values of energy in the
time–frequency domain.

In this section, we propose to invert the scattering transform operator; that
is, for a given shape vector $\boldsymbol{\theta}$, to recover the drum sound
$\boldsymbol{x}_{\boldsymbol{\theta}}$ from its scattering coefficients
$\mathbf{S}\boldsymbol{x}_{\boldsymbol{\theta}}$. The long-term motivation
behind this approach is to convert wav2shape, which is a discriminative model,
into a generative model. Indeed, although the synthesis of drum sounds from
simple shapes can be resolved by solving a partial differential equation in
closed form (see Section 3), the case of irregular shapes is mathematically
intractable. One way to circumvent this problem is to adopt a data-driven
perspective on physical simulations, and let a machine learning model decode
the latent vector $\boldsymbol{\theta}$ into its acoustical correlate
$\boldsymbol{x}_{\boldsymbol{\theta}}$.

Audio texture synthesis from scattering coefficients was inaugurated by [18].
We refer to [19] for a practical application of scattering-based audio
synthesis to contemporary music creation. More recently, the combination of
scattering transform and generative adversarial networks (GAN) has shown a
promising avenue of research [20, 21].

Given a drum shape vector $\boldsymbol{\theta}$, our goal is to minimize the
error functional

$E(\boldsymbol{y})=\big{\|}\mathbf{S}\boldsymbol{y}-\mathbf{S}\boldsymbol{x}_{\boldsymbol{\theta}}\big{\|},$
(27)

defined as a Euclidean distance in the domain of scattering coefficients. Let
us denote by $\boldsymbol{\nabla}E(\boldsymbol{y})$ the gradient of $E$
evaluated at $\boldsymbol{y}$. In practice, we compute
$\boldsymbol{\nabla}E(\boldsymbol{y})$ by reverse-mode automatic
differentiation, via the PyTorch backend of the Kymatio library. We refer to
[22] for a description of gradient backpropagation in a scattering network,
with an application to audio texture synthesis.

Because the function $E$ is nonconvex, gradient descent does not converge to
the global minimizer
$\boldsymbol{y}^{*}=\boldsymbol{x}_{\boldsymbol{\theta}}$, but merely to a
local minimizer. Whether this local minimizer is perceptually to the target
signal $\boldsymbol{x}_{\boldsymbol{\theta}}$ depends upon the parameters of
the scattering transform.

Starting from a colored Gaussian noise $\boldsymbol{y}_{0}(t)$ whose power
spectral density matches $\mathbf{S}_{1}\boldsymbol{x}(t,\lambda)$, we refine
it by additive updates of the form

$\boldsymbol{y}_{n+1}(t)=\boldsymbol{y}_{n}(t)+\boldsymbol{u}_{n}(t).$ (28)

In the simplest form of gradient descent, the signal $\boldsymbol{u}_{n}$ is
equal to the gradient $\boldsymbol{\nabla}E(\boldsymbol{y_{n}})$ multiplied by
a constant learning rate term $\mu$.

To accelerate convergence, we adopt an adaptive learning rate policy, known as
the “bold driver” heuristic [23]. We initialize the learning rate at
$\mu_{0}=0.1$. At every iteration $n$, if the error has reduced, we confirm
the update and increase the learning rate by $10\%$; otherwise, we retract the
step and decrease the learning rate by $50\%$. This leads to a sequence of
learning rates: $\mu_{1},\mu_{2},$ etc. which depends upon the iteration.

Moreover, we add a momentum term to the gradient term
$\boldsymbol{\nabla}E(\boldsymbol{y_{n}})$, leading to an update of the form:

$\boldsymbol{u}_{n}(t)=m\times\boldsymbol{u}_{n-1}(t)+\mu_{n}\boldsymbol{\nabla}E(\boldsymbol{y}_{n})(t).$
(29)

Following a previous publication [24, Section 3.3.2], we set the momentum
hyperparameter to $m=0.9$.

After $n\sim 300$ iterations, the reconstruction error $E(\boldsymbol{y}_{n})$
is two orders of magnitude below the initial reconstruction error
$E(\boldsymbol{y}_{0})$. We repeat the same procedure for scattering networks
of depths $N=1$ and $N=2$ and for various values of the scattering scale
parameter $J$.

Figure 7 illustrates our results for one particular setting of the drum shape
vector $\boldsymbol{\theta}$. In the case of first-order scattering, we find
that lower values of $J$ lead to a sharper reconstruction of
$\boldsymbol{x}_{\boldsymbol{\theta}}$, at the expense of stability to
deformations (see Section 2.4). Conversely, larger values of $J$ elicit
audible artifacts in the reconstructed signal $\boldsymbol{y}$. Nevertheless,
increasing the depth of the scattering network from $N=1$ to $N=2$ layers
yields an improvement of reconstruction quality, which is particularly
noticeable with higher values of $J$. Although previous publications had
reported the same effect in the case of audio textures [18], such as applause
or bubbling water, the originality of our finding is that it applies to
isolated percussive events.

##  6 Conclusions

In this article, we have presented wav2shape: a convolutional neural network
which disentangles and retrieves physical parameters from waveforms of
percussive sounds. First, we have presented a 2-D physical model of a
rectangular membrane, based on a fourth-order partial differential equation
(PDE) in time and space. We have solved the PDE in closed form by means of the
functional transformation method (FTM), and included a freely downloadable VST
plugin which synthesizes drum sounds in real time. Then, we have computed
second-order scattering coefficients of these sounds and designed wav2shape as
a convolutional neural network (CNN) operating on the logarithm of these
coefficients. We have trained wav2shape in a supervised fashion in order to
regress the parameters underlying the PDE, such as pitch, sustain, and
inharmonicity.

From an experimental standpoint, we have found that wav2shape is capable of
generalizing beyond its training set and predicting the shape of previously
unseen sounds (Figure 3). The network’s robustness in shape regression
confirmed that the scattering transform has the ability to linearize the
dependency of the signal upon the position of the drum stroke (Figure 5).
Indeed, when applied on linearly interpolated scattering coefficients, the
wav2shape neural network continues to produce an interpretable outcome.
Lastly, we have used reverse-mode automatic differentiation in the Kymatio
library to synthesize drum sounds directly from scattering coefficients,
without explicitly solving a partial differential equation.

Although the results of wav2shape are promising, we acknowledge that it
suffers from some practical limitations, which hamper its usability in
computer music creation. First, physical parameters such as inharmonicity $D$
and aspect ratio $\alpha$ are not recovered as accurately as pitch $\omega$ or
sustain $\tau$. Secondly, wav2shape is only capable of retrieving the shape
vector $\boldsymbol{\theta}$ if the rectangular drum is stroked exactly at its
center: it would be beneficial, albeit challenging, to generalize the approach
to any stroke location $u_{0}$. Thirdly, we have trained wav2shape on a
relatively large training set of over 82k audio samples. The acquisition of
these samples was only made possible by simulating the response of the
membrane. The prospect of extending autonomous systems from such a simulated
environment towards a real environment is a topic of ongoing research in
reinforcement learning, known as sim2real [25]. Yet, the field of deep
learning for musical acoustics predominantly relies on supervised learning
techniques instead of reinforcement learning. In this context, we believe that
future research is needed to strengthen the interoperability between physical
modeling and data-driven modeling of musical sounds.

##  7 Acknowledgment

This work is partially supported by National Science Foundation award 1633259
(BIRDVOX). We wish to thank Jennie Choi for managing the open-access API of
the Metropolitan Museum of Art and Ivan Selesnick for generously providing his
notes on digital sound synthesis using the functional transformation method.
We also thank Scott Fitzgerald, Amy Hurst, and Mark Plumbley for fruitful
discussions.

## References

  * [1] C. Wu, C. Dittmar, C. Southall, R. Vogl, G. Widmer, J. Hockman, M. Müller, and A. Lerch, “A review of automatic drum transcription,” IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 26, no. 9, pp. 1457–1483, 2018.
  * [2] L. Trautmann and R. Rabenstein, Digital sound synthesis by physical modeling using the functional transformation method.  Springer, 2013.
  * [3] L. Hiller and P. Ruiz, “Synthesizing musical sounds by solving the wave equation for vibrating objects: Part 2,” Journal of the Audio Engineering Society, vol. 19, no. 7, pp. 542–551, 1971.
  * [4] A. Chaigne and A. Askenfelt, “Numerical simulations of piano strings. I. A physical model for a struck string using finite difference methods,” The Journal of the Acoustical Society of America, vol. 95, no. 2, pp. 1112–1118, 1994.
  * [5] J. O. Smith, Physical Audio Signal Processing.  W3K Editions, 2010.
  * [6] J.-M. Adrien, The Missing Link: Modal Synthesis, p. 269–298.  Cambridge, MA, USA: MIT Press, 1991.
  * [7] V. Lostanlen, J. Andén, and M. Lagrange, “Extended playing techniques: The next milestone in musical instrument recognition,” in Proceedings of the International Conference on Digital Libraries for Musicology (DLfM), 2018\.
  * [8] V. Lostanlen, A. Cohen-Hadria, and J. P. Bello, “One or Two Components? The Scattering Transform Answers,” in Proceedings of the European Signal Processing Conference (EUSIPCO), 2020.
  * [9] M. Kac, “Can one hear the shape of a drum?,” The American Mathematical Monthly, vol. 73, no. 4P2, pp. 1–23, 1966.
  * [10] R. Rabenstein, “Digital sound synthesis of string instruments with the functional transformation method,” Signal Processing, vol. 83, no. 8, pp. 1673 – 1688, 2003.
  * [11] L. Trautmann, S. Petrausch, and R. Rabenstein, “Physical modeling of drums by transfer function methods,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), vol. 5, pp. 3385–3388 vol.5, 2001.
  * [12] M. Schäfer, P. Frenstatsky, and R. Rabenstein, “A physical string model with adjustable boundary conditions,” in Proceedings of the International Conference on Digital Audio Effects (DAFX), 2016.
  * [13] M. Schafer and R. Rabenstein, “Calculation of the transformation kernels for the functional transformation method,” in Proceedings of the International Workshop on Multidimensional Systems (nDS), pp. 1–6, 2017.
  * [14] R. Rabenstein, T. Koch, and C. Popp, “Tubular bells: A physical and algorithmic model,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 18, no. 4, pp. 881–890, 2010.
  * [15] J. Andén and S. Mallat, “Deep scattering spectrum,” IEEE Transactions on Signal Processing, vol. 62, p. 4114–4128, Aug 2014.
  * [16] S. Mallat, “Group invariant scattering,” Communications on Pure and Applied Mathematics, vol. 65, no. 10, pp. 1331–1398, 2012.
  * [17] M. Andreux, T. Angles, G. Exarchakis, R. Leonarduzzi, G. Rochette, L. Thiry, J. Zarka, S. Mallat, J. Andén, E. Belilovsky, et al., “Kymatio: Scattering transforms in Python,” Journal of Machine Learning Research, vol. 21, no. 60, pp. 1–6, 2020.
  * [18] J. Bruna and S. Mallat, “Audio texture synthesis with scattering moments,” arXiv preprint arXiv:1311.0407, 2013.
  * [19] V. Lostanlen, “On time-frequency scattering and computer music,” in Florian Hecker: Halluzination, Perspektive, Synthese (V. J. Müller, ed.), Berlin: Sternberg Press, 2019.
  * [20] M. Andreux and S. Mallat, “Music generation and transformation with moment matching–scattering inverse networks,” in Proceedings of the International Society on Music Information Retrieval (ISMIR) Conference, pp. 327–333, 2018.
  * [21] T. Angles and S. Mallat, “Generative networks as inverse problems with scattering transforms,” in Proceedings of the International Conference on Learning Representations (ICLR), 2018.
  * [22] V. Lostanlen and F. Hecker, “The Shape of RemiXXXes to Come: Audio Texture Synthesis with Time–frequency Scattering,” in Proceedings of the International Conference on Digital Audio Effects (DAFX), 2019.
  * [23] I. Sutskever, J. Martens, G. Dahl, and G. Hinton, “On the importance of initialization and momentum in deep learning,” in Proceedings of the International Conference on Machine Learning (ICML), pp. 1139–1147, 2013.
  * [24] V. Lostanlen, Convolutional operators in the time-frequency domain.  PhD thesis, École normale supérieure, 2017.
  * [25] Y. Chebotar, A. Handa, V. Makoviychuk, M. Macklin, J. Issac, N. Ratliff, and D. Fox, “Closing the sim-to-real loop: Adapting simulation randomization with real-world experience,” in Proceedings of the International Conference on Robotics and Automation (ICRA), pp. 8973–8979, IEEE, 2019.


# Digital Twin-Assisted Collaborative Transcoding for Better User Satisfaction
in Live Streaming

Xinyu Huang1, Mushu Li2, Wen Wu3, Conghao Zhou1, and Xuemin (Sherman) Shen1
1Department of Electrical & Computer Engineering, University of Waterloo,
Canada
2Department of Electrical, Computer, and Biomedical Engineering, Toronto
Metropolitan University, Canada
3Frontier Research Center, Peng Cheng Laboratory, China
Email: {x357huan, c89zhou<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a digital twin (DT)-assisted cloud-edge
collaborative transcoding scheme to enhance user satisfaction in live
streaming. We first present a DT-assisted transcoding workload estimation
(TWE) model for the cloud-edge collaborative transcoding. Particularly, two
DTs are constructed for emulating the cloud-edge collaborative transcoding
process by analyzing spatial-temporal information of individual videos and
transcoding configurations of transcoding queues, respectively. Two light-
weight Bayesian neural networks are adopted to fit the TWE models in DTs,
respectively. We then formulate a transcoding-path selection problem to
maximize long-term user satisfaction within an average service delay
threshold, taking into account the dynamics of video arrivals and video
requests. The problem is transformed into a standard Markov decision process
by using the Lyapunov optimization and solved by a deep reinforcement learning
algorithm. Simulation results based on the real-world dataset demonstrate that
the proposed scheme can effectively enhance user satisfaction compared with
benchmark schemes.

##  I Introduction

With the prevalence of smart mobile equipment and ubiquitous Internet access,
an increasing number of amateur or professional broadcasters start to utilize
online video platforms, such as Youtube Live, Facebook Live, Twitch TV, etc.,
to produce live streams and interact with users anytime and anywhere.
According to a recent Grand View Research report, the live streaming market is
expected to increase from $70 billion in 2021 to about $224 billion in 2028
[1]. As a critical technology to guarantee continuous and high-definition
video playback, video transcoding encodes video streams pipelined at online
video platforms into multiple video versions in real time [2]. Since a user’s
device has a small buffer size for live streaming playback, the performance of
video transcoding can directly affect user satisfaction. Video transcoding can
be conducted in a cloud server, i.e., cloud-transcoding, and an edge server,
i.e., edge-transcoding [3]. However, video transcoding is a computation-
intensive process, and the transcoding delay dominates the service delay,
almost $70\%$ [4]. Relying solely on cloud-transcoding or edge-transcoding can
incur a large transmission delay and heavy computation overhead, respectively
[5, 6]. Therefore, efficient cloud-edge collaborative transcoding has
attracted considerable research attention.

In the literature, significant efforts have been devoted to improve the cloud-
edge collaborative transcoding performance. Pang et al. proposed to dispatch
newly generated video streams to appropriate transcoding queues (TQs) in cloud
and edge servers based on predicted interaction intensities in different
streaming channels, which can satisfy users’ heterogeneous quality of
experience (QoE) requirements [7]. To reduce collaborative transcoding cost,
Erfanian et al. constructed a transcoding-cost-based multicast tree for each
cloud and edge computing server to determine where and how many TQs should be
deployed [8]. Zhu et al. proposed an auction-based approach for TQ selection
to further reduce transcoding cost within a prescribed transcoding delay
threshold [5]. The above works utilize a general transcoding workload
estimation (TWE) model for different kinds of video streams in different TQs.
However, the TWE model ignores videos’ spatial-temporal information and
servers’ transcoding configurations, thus bringing TQs’ length estimation
errors. Such estimation errors can result in improper transcoding-path
selections for video streams, thereby rendering user satisfaction degradation.
Hence, it is paramount to construct an accurate TWE model.

The digital twin (DT) technology is a potential solution since it is a digital
representation of a physical entity (PE) that can accurately reflect its
status and feature via real-time synchronization between the DT and PE [9].
The DT technology can be utilized to store and analyze users’ data to
construct user-specific and real-time personalized QoE models for tailored
network resource management [10]. We leverage the DT technology to construct
virtual TQs based on estimated transcoding workloads and transcoding-path
selections, which can reflect the dynamics of physical TQs and emulate
transcoding performance of physical TQs in real time.

In this paper, we present a DT-assisted cloud-edge collaborative transcoding
scheme for live streaming services to enhance long-term user satisfaction
given an average service delay threshold. Specifically, we first propose a DT-
assisted collaborative transcoding model to capture the dynamics of physical
TQs in cloud and edge servers. The cloud and edge transcoding DTs are
established at the network controller, consisting of historical transcoding
data, a TWE model and virtual TQs. The TWE model adopts a Bayesian neural
network (BNN) to fit historical transcoding data in order to reduce estimation
errors. Based on the TWE model, the virtual TQs are constructed to reflect the
dynamics of physical TQs. Secondly, we design a tailored deep reinforcement
learning (DRL) algorithm for making transcoding-path selections. The objective
is to maximize long-term user satisfaction within an average service delay
threshold by optimizing transcoding-path selections. Since the formulated
transcoding-path selection problem is a constrained Markov decision process
(CMDP), we leverage the Lyapunov optimization technique to transform it into a
standard MDP. Extensive experiments conducted on the real-world dataset
demonstrate that the proposed DT-assisted collaborative transcoding scheme can
effectively enhance user satisfaction within an service delay threshold
compared with benchmarks. The main contributions of this paper are summarized
as follows:

  * $\bullet$

We propose a DT-assisted collaborative transcoding model, which can reduce the
TQ length estimation error and capture its dynamics.

  * $\bullet$

We develop a cloud-edge collaborative transcoding algorithm based on DRL,
which can realize an online transcoding-path selection.

The remainder of this paper is organized as follows. The DT-assisted
collaborative transcoding scheme is presented in Section II. The DRL-based
transcoding scheduling algorithm is proposed in Section III. Simulation
results are provided in Section IV, followed by the conclusion in Section V.

##  II DT-Assisted Collaborative Transcoding Scheme

Figure 1: DT-assisted collaborative transcoding scenario.

###  II-A System Model

As shown in Fig. 1, we consider a DT-assisted collaborative transcoding
scenario, which mainly consists of a cloud server, an edge server, and DTs.

Users’ devices generate video requests to the network controller for
requesting live streams to be played soon with different video qualities every
0.5 s. The cloud and edge servers collaboratively transcode requested live
streams and deliver them to users’ devices. Specifically, the cloud server is
deployed with three kinds of TQs111Our proposed scheme can also be applied to
the scenario where the number of TQs in the cloud or edge server is more than
three. for video transcoding based on different encoding
presets222https://trac.ffmpeg.org/wiki/Encode/H.264
https://trac.ffmpeg.org/wiki/Encode/H.265, i.e., slow TQ, medium TQ, and fast
TQ. From the fast to slow TQs, the transcoding speed gradually decreases but
the transcoded video quality gradually rises. A higher video quality usually
corresponds to a larger video bit rate. The edge server is deployed with two
kinds of TQs, i.e., medium TQ and fast TQ, considering that the edge server
with limited computing resources hardly supports a slow TQ that requires
plenty of computing resources. As the minimal transcoding unit of video
streams, a group of pictures (GoP) is dispatched to different TQs in the cloud
and edge servers to generate multiple video versions of different qualities.
When a GoP is delivered from the cloud server to a user’s device, there exist
two kinds of transcoding paths, i.e., one-step transcoding path and two-step
transcoding path, according to the queues that the GoP goes through. The
former path refers to the original GoP being transcoded in any cloud TQ and
then directly sent to users. The latter path indicates that the original GoP
is first transcoded in the high or medium cloud TQ, and then dispatched to any
edge TQ for further transcoding. By transcoding the GoP into multiple video
versions, users’ differentiated requests can be satisfied.

To properly select the video transcoding path, two kinds of DTs are
constructed, i.e., cloud transcoding DT (CTDT) and edge transcoding DT (ETDT),
and located at the network controller. Both CTDT and ETDT consist of three
components: historical transcoding data on GoPs’ spatial-temporal information
and servers’ transcoding configurations, a unique TWE model and multiple
virtual TQs. The TWE model and virtual TQs are used to characterize the
transcoding workload of each GoP and the dynamics of physical TQs,
respectively. Specifically, the TWE model adopts a BNN to analyze the
historical transcoding data, which can estimate the specific transcoding
workload for each GoP in different TQs. To ensure the estimation accuracy of
TWE models, actual transcoding data are recorded in the cloud and edge servers
if the transcoding workload bias exceeds a prescribed threshold. The recorded
transcoding data are periodically uploaded to the network controller for the
BNN model update. In addition, the virtual TQs in both CTDT and ETDT are
digital representations of physical TQs in the cloud and edge servers, which
are updated based on selected transcoding paths and estimated transcoding
workloads of GoPs. The cloud and edge servers also send the queue
synchronization message if the TQ length bias exceeds a prescribed threshold.

When users send GoP requests consisting of GoP indexes and corresponding video
bit rates to the network controller, CTDT and ETDT collect GoPs’ transcoding
data, and use respective TWE models to estimate corresponding transcoding
workloads in different TQs. Then, the network controller makes the
transcoding-path selection for each GoP by taking GoP requests, estimated
transcoding workloads, and virtual TQ lengths into account. TWE models and
virtual TQs are updated based on the recorded transcoding data and determined
transcoding paths of GoPs, respectively.

###  II-B DT-Assisted GoP’s TWE

CTDT and ETDT can emulate the dynamics of physical TQs for newly arrived GoPs
by analyzing the GoP’s spatial-temporal information, the servers’ transcoding
configurations, and the computing capability of each TQ.

Spatial and temporal information: The spatial information (SI) and temporal
information (TI) are critical factors that can affect the transcoding workload
[11]. The SI indicates the amount of spatial details in a video frame.
Specifically, each video frame in a GoP is first filtered by the Sobel filter,
and then the standard deviation over the pixels in each filtered frame is
calculated. The maximum standard deviation among filtered frames represents
the GoP’s SI. Therefore, the SI of GoP $k$ can be calculated by

${{\varsigma}_{k}}=\underset{l\in\mathcal{L}}{\mathop{\max}}\,\left\\{\sigma\left[\Theta(F_{k}^{l})\right]\right\\},$
(1)

where $F_{k}^{l}$ is video frame $l$ of GoP $k$, and $\Theta(\cdot)$ is the
Sobel filter operation. In addition, $\sigma$ is the standard deviation of a
filtered frame, and $\mathcal{L}$ is the set of GoP frames.

The TI indicates the amount of temporal changes of a video frame sequence.
Specifically, the pixel difference of each adjacent frames is first
calculated. Then, the standard deviation of each pixel difference is
calculated. The maximum standard deviation is chosen to represent the GoP’s
TI. Therefore, the TI of GoP $k$ can be calculated by

${{\xi}_{k}}=\underset{l\in\mathcal{L}\setminus\\{L\\}}{\mathop{\max}}\,\left\\{\sigma\left[F_{k}^{l+1}-F_{k}^{l}\right]\right\\}.$
(2)

TWE: To emulate the transcoding performance of each TQ for newly arrived GoPs,
we select the BNN to construct TWE models in both CTDT and ETDT. The BNN is a
widely used tool to fit the function by finding the distribution of the
weighting vector and uses the regularization technology to avoid overfitting
[12].

The input of the BNN-based model is a multi-dimension vector, denoted by
${{V}}$, which mainly abstracts a GoP’s spatial-temporal information, servers’
transcoding configurations, and TQs’ computing capabilities, as presented in
Table I.

TABLE I: The input vector of the BNN-based DT model Index | Name | Index | Name
---|---|---|---
1 | Encoding Preset | 5 | Computing Processor
2 | Number of Frames | 6 | Computing Density
3 | Resolution | 7 | Computing Capability
4 | SI | 8 | TI

The output of the BNN-based model is the estimated transcoding workload for a
GoP in TQ $i$, which can be expressed as ${{\Omega}_{i}}\left({{V}}\right)$.
The function ${{\Omega}_{i}}(\cdot)$ is fitted based on the BNN regularization
algorithm and updated periodically to decrease the estimation error. The
recorded actual transcoding workloads are selected as ground truths to
calculate the mean squared error (MSE) for model training and update.

###  II-C DT-Assisted Collaborative Transcoding Model

After estimating transcoding workloads of GoPs in different TQs in the cloud
server and edge server, the network controller will select the appropriate
transcoding path for each GoP based on GoP requests, estimated transcoding
workloads, and virtual TQ lengths. A simplified procedure is shown in Fig. 2.
There exist six transcoding paths, including three one-step transcoding paths,
i.e., path 1: slow cloud-transcoding, path 3: medium cloud-transcoding, and
path 6: fast cloud-transcoding, and three two-step transcoding paths, i.e.,
path 2: slow-medium cloud-edge transcoding, path 4: medium-fast cloud-edge
transcoding, and path 5: slow-fast cloud-edge transcoding.

Figure 2: Collaborative transcoding procedure.

The collaborative transcoding decision is made in each scheduling slot,
indexed by $t$. We consider that there are $K$ GoPs arriving at the cloud
server at scheduling slot $t$, and the corresponding set of GoPs is
$\mathcal{K}$. The sets of cloud TQs and edge TQs are denoted by
${{\Lambda}_{1}}=\left\\{1,2,3\right\\}$ and
${{\Lambda}_{2}}=\left\\{4,5\right\\}$, respectively, which correspond to the
sets of CTDT queues and ETDT queues. To determine the transcoding path for
each GoP, we define a binary transcoding position variable $x_{t,k}^{i}$,
where $x_{t,k}^{i}=1$ indicates that GoP $k$ is dispatched to TQ $i$ for
transcoding at scheduling slot $t$; Otherwise, $x_{t,k}^{i}=0$. The set of
index $i$ ranges from $1$ to $5$, which refers to the slow, medium, and fast
CTDT queues, and the medium and fast ETDT queues, respectively.

CTDT queue dynamics: In scheduling slot $t$, the queue length of TQ $i$ in
CTDT, denoted by, $L_{t}^{i}$, is updated via

$L_{t+1}^{i}={{\left[L_{t}^{i}+\sum\limits_{k=1}^{K}{x_{t,k}^{i}\frac{{{\Omega}_{1}}\left(V_{t}^{k}\right)b_{t}^{k}}{{{f}_{i}}{{\kappa}_{i}}}}-d\right]}^{+}},\forall
i\in{{\Lambda}_{1}},$ (3)

where ${\Omega}_{1}\left(V_{t}^{k}\right)$ is the estimated transcoding
workload for GoP $k$ in TQ 1, and $b_{t}^{k}$ is the original bit rate of GoP
$k$ at scheduling slot $t$. Parameters ${{f}_{i}}$ and ${{\kappa}_{i}}$
represent the computing capability of TQ $i$ and the corresponding computing
density, respectively. Parameter $d$ is the scheduling slot length.

ETDT queue dynamics: To reflect the queue dynamics of ETDT queues, we count
how many GoPs are dequeued in each CTDT queue and sent to each ETDT queue
based on two-step transcoding paths in each scheduling slot. For the medium
ETDT queue, i.e., TQ 4, enqueued GoPs origin from path 2, and the
corresponding queue length is updated via

$L_{t+1}^{4}={{\left[L_{t}^{4}+\frac{\alpha_{t}^{1}\min\left\\{L_{t}^{1},T\right\\}\overline{\Omega}_{t}^{1\to
4}{{f}_{1}}{{\kappa}_{1}}}{\overline{\Omega}_{t}^{1}{{f}_{4}}{{\kappa}_{4}}}-d\right]}^{+}},$
(4)

where $\overline{\Omega}_{t}^{1}$ is the average transcoding workload of GoPs
in the slow CTDT queue at scheduling slot $t$. Here, $\alpha_{t}^{1}$
indicates the ratio of GoPs of path 2 in the slow CTDT queue, i.e., TQ 1,
which is updated via

$\alpha_{t+1}^{1}=\frac{\alpha_{t}^{1}L_{t}^{1}+\sum\nolimits_{k=1}^{K}{x_{t,k}^{4}\frac{{{\Omega}_{1}}\left(V_{t}^{k}\right)b_{t}^{k}}{{{f}_{1}}{{\kappa}_{1}}}}}{L_{t}^{1}+\sum\nolimits_{k=1}^{K}{x_{t,k}^{1}\frac{{{\Omega}_{1}}\left(V_{t}^{k}\right)b_{t}^{k}}{{{f}_{1}}{{\kappa}_{1}}}}}.$
(5)

In addition, $\overline{\Omega}_{t}^{1\to 4}$ is the average transcoding
workload of GoPs of path 2 in TQ 1, which is updated via

$\overline{\Omega}_{t+1}^{1\to
4}=\frac{\sum\nolimits_{k=1}^{K}{x_{t,k}^{1}x_{t,k}^{4}{{\Omega}_{4}}(V_{t}^{k})b_{t}^{k}}+\overline{\Omega}_{t}^{1\to
4}B_{t}^{1\to
4}}{\sum\nolimits_{k=1}^{K}{x_{t,k}^{1}x_{t,k}^{4}s_{t}^{k}}+B_{t}^{1\to 4}},$
(6)

where $B_{t}^{1\to 4}$ is the bit rate of all GoPs belonging to path 2 at
scheduling slot $t$.

For the fast ETDT queue, the enqueued GoPs can be sent from path 4 and path 5.
Therefore, the fast ETDT queue, i.e., TQ 5, is updated via

$\begin{split}L_{t+1}^{5}=&{{\left[L_{t}^{5}+\frac{\alpha_{t}^{2}\min\left\\{L_{t}^{1},T\right\\}\overline{\Omega}_{t}^{1\to
5}{{f}_{1}}{{\kappa}_{1}}}{\overline{\Omega}_{t}^{1}{{f}_{5}}{{\kappa}_{5}}}\right.}}\\\
&{{\left.+\frac{\alpha_{t}^{3}\min\left\\{L_{t}^{2},T\right\\}\overline{\Omega}_{t}^{2\to
5}{{f}_{2}}{{\kappa}_{2}}}{\overline{\Omega}_{t}^{2}{{f}_{5}}{{\kappa}_{5}}}-d\right]}^{+}},\end{split}$
(7)

where $\overline{\Omega}_{t}^{2}$ is the average transcoding workload of GoPs
in the medium CTDT queue, i.e., TQ 2, at scheduling slot $t$. Here,
$\alpha_{t}^{2}$ and $\alpha_{t}^{3}$ refer to the ratio of GoPs of path 4 and
5 in TQ 1 and TQ 2, respectively, which are updated via

$\alpha_{t+1}^{i}=\frac{\alpha_{t}^{i}L_{t}^{i}+\sum\nolimits_{k=1}^{K}{x_{t,k}^{5}\frac{{{\Omega}_{i-1}}(V_{t}^{k}){{b}_{t,k}}}{{{f}_{i-1}}{{\kappa}_{i-1}}}}}{L_{t}^{i-1}+\sum\nolimits_{k=1}^{K}{x_{t,k}^{i-1}\frac{{{\Omega}_{i-1}}(V_{t}^{k}){{b}_{t,k}}}{{{f}_{i-1}}{{\kappa}_{i-1}}}}},\forall
i\in\left\\{2,3\right\\}.$ (8)

In addition, $\overline{\Omega}_{t}^{1\to 5}$ and $\overline{\Omega}_{t}^{2\to
5}$ indicate the average transcoding workloads of GoPs of path 4 and path 5 in
TQ 1 and TQ 2, respectively, which are updated by

$\overline{\Omega}_{t+1}^{i\to
5}\\!=\\!\frac{\sum\nolimits_{k=1}^{K}{x_{t,k}^{i}\\!x_{t,k}^{5}{{\Omega}_{5}}\left(V_{t}^{k}\right){{b}_{t,k}}}+\overline{\Omega}_{t}^{i\to
5}B_{t}^{i\to
5}}{\sum\nolimits_{k=1}^{K}{x_{t,k}^{i}x_{t,k}^{5}{{b}_{t,k}}}+B_{t}^{i\to
5}},\forall i\\!\in\\!\left\\{1,2\right\\}.$ (9)

Service delay: Based on the analysis of queuing dynamics in CTDT and ETDT, we
can estimate the service delay, $D_{t}(X_{t})$, in this cloud-edge
collaborative transcoding system, which is given by

${{D}_{t}}({{X}_{t}})=\frac{1}{3}\sum\nolimits_{i\in{{\Lambda}_{1}}}{L_{t}^{i}}+{{I}_{t}}+\frac{1}{2}\sum\nolimits_{i\in{{\Lambda}_{2}}}{L_{t}^{i}},$
(10)

where ${{X}_{t}}$ is the decision variable set, and
${{X}_{t}}={{\left\\{x_{t,k}^{i}\right\\}}_{i\in\mathcal{I},k\in\mathcal{K}}}$.
Here, ${{I}_{t}}$ is the transmission delay between the cloud server and the
edge server at scheduling slot $t$, which can be estimated based on the round-
trip time (RTT).

User satisfaction: In addition to the service delay, we also introduce user
satisfaction to evaluate transcoding performance. The user satisfaction,
$W_{t}$, refers to the ratio of users’ video requests that can be satisfied
through transcoding operations in scheduling slot $t$, which is depicted as

${{W}_{t}}=\frac{\sum\nolimits_{k=1}^{K}{x_{t,k}^{1}w_{t,k}^{1}+\left(x_{t,k}^{2}+x_{t,k}^{1}x_{t,k}^{4}\right)w_{t,k}^{2}+M_{t,k}}}{\sum\nolimits_{k=1}^{K}{w_{t,k}^{1}+w_{t,k}^{2}+w_{t,k}^{3}}},$
(11)

where $w_{t,k}^{1},w_{t,k}^{2},w_{t,k}^{3}$ represent the number of GoP
requests of high, medium, and low quality, respectively. Here,
$M_{t,k}=\left(x_{t,k}^{3}+x_{t,k}^{1}x_{t,k}^{5}+x_{t,k}^{2}x_{t,k}^{5}\right)w_{t,k}^{3}$.

###  II-D Problem Formulation

Our objective is to maximize long-term user satisfaction within an average
service delay threshold over $T$ scheduling slots. Correspondingly, the
optimization problem is formulated as

$\displaystyle{{\textbf{P}}_{\text{0}}:}\quad$
$\displaystyle\underset{{{X}_{t}}}{\mathop{\max}}\,\underset{T\to\infty}{\mathop{\lim}}\,\frac{1}{T}\sum\limits_{t=1}^{T}{{{W}_{t}}({{X}_{t}})},$
(12) s.t.
$\displaystyle\underset{T\to\infty}{\mathop{\lim}}\,\frac{1}{T}\sum\limits_{t=1}^{T}{{{D}_{t}}({{X}_{t}})}\leq\overline{D},$
(12a) $\displaystyle x_{t,k}^{2}+x_{t,k}^{4}\leq 1,$ (12b) $\displaystyle
x_{t,k}^{3}+x_{t,k}^{5}\leq 1,$ (12c) $\displaystyle
x_{t,k}^{1}-x_{t,k}^{4}\geq 0,$ (12d) $\displaystyle
x_{t,k}^{1}+x_{t,k}^{2}-x_{t,k}^{5}\geq 0,$ (12e) $\displaystyle
x_{t,k}^{i}\in\left\\{0,1\right\\},\forall i\in\mathcal{I}\text{,}$ (12f)

where $\overline{D}$ is the average service delay threshold. Constraint (12a)
represents the average service delay requirement. Constraints (12b,12c) avoid
the repeated transcoding operation for a GoP between cloud TQs and edge TQs.
Constraint (12d) guarantees that a GoP is transcoded from the high-quality
version to the medium-quality version. Constraint (12e) guarantees that a GoP
is transcoded from the high-quality and medium-quality versions to the low-
quality version.

The formulated problem is a CMDP since the state transition is Markovian with
the long-term constraint. However, the CMDP cannot be directly solved by
general DRL algorithms due to the intractable long-term constraint [13].
Therefore, we first transform it into a standard MDP and then develop a DRL
algorithm to solve it.

##  III DRL-Based Transcoding Scheduling Algorithm

###  III-A Problem Transformation

The Lyapunov optimization technique is an effective way to handle a long-term
constraint [14] in an optimization problem. The main idea is to establish a
deficit queue of service delay to characterize the satisfaction status of the
long-term constraint, which can transform the time average constraint problem
into a queue stability problem.

First, we establish the deficit queue of service delay, as follows:

${{Z}_{t+1}}={{\left[{{D}_{t}}\left({{X}_{t}}\right)-\overline{D}+{{Z}_{t}}\right]}^{+}},$
(13)

where ${{Z}_{t}}$ represents the deviation of achieved instantaneous values
from the long-term constraint, whose initial state is set to 0. To
characterize the satisfaction status of the long-term constraint, the Lyapunov
function is defined as $L({{Z}_{t}})=\frac{1}{2}{{({{Z}_{t}})}^{2}}$[14]. A
smaller Lyapunov function value indicates a better satisfaction of the long-
term constraint.

Second, to guarantee that the Lyapunov function can be consistently preserved
within a small value, the one-shot Lyapunov drift is introduced to capture the
variation of Lyapunov function values across two subsequent time slots.
Therefore, the one-shot Lyapunov drift is defined as
$\Delta({{Z}_{t}})=L({{Z}_{t+1}})-L({{Z}_{t}})$. The upper bound of
$\Delta({{Z}_{t}})$ can be derived as

$\small\begin{split}&\Delta({{Z}_{t}})=\frac{1}{2}{{({{Z}_{t+1}})}^{2}}-\frac{1}{2}{{({{Z}_{t}})}^{2}}\\\
&\leq\frac{1}{2}{{\left({{Z}_{t}}+{{D}_{t}}\left({{X}_{t}}\right)-\overline{D}\right)}^{2}}-\frac{1}{2}{{\left({{Z}_{t}}\right)}^{2}}\\\
&={{Z}_{t}}\left({{D}_{t}}\left({{X}_{t}}\right)-\overline{D}\right)+\frac{1}{2}{{\left({{D}_{t}}\left({{X}_{t}}\right)-\overline{D}\right)}^{2}}\\\
&={{Z}_{t}}\left({{D}_{t}}\left({{X}_{t}}\right)-\overline{D}\right)+\frac{1}{2}{{\left(\frac{1}{3}\sum\limits_{i=1}^{3}{L_{t}^{i}}+{{I}_{t}}+\frac{1}{2}\sum\limits_{i=4}^{5}{L_{t}^{i}}-\overline{D}\right)}^{2}}\\\
&\leq{{Z}_{t}}\left({{D}_{t}}\left({{X}_{t}}\right)-\overline{D}\right)+\Gamma,\end{split}$
(14)

where
$\Gamma=\frac{1}{2}{{\left(\frac{1}{3}\sum\limits_{i=1}^{3}{L_{i}^{\max}}+{{I}_{\max}}+\frac{1}{2}\sum\limits_{i=4}^{5}{L_{i}^{\max}}-\overline{D}\right)}^{2}}$
and it is a constant. Here, $L_{i}^{\max}$ and ${{{I}}_{\max}}$ are the
maximum length of TQ $i$ and the maximum RTT, respectively.

Third, the original optimization problem of maximizing long-term user
satisfaction within an average service delay threshold can be transformed to
minimize a drift-plus-cost. Problem $\textbf{P}_{0}$ is reformulated as

$\begin{split}\textbf{P}_{\text{1}}:&\quad\underset{{{X}_{t}}}{\mathop{\min}}\,{{Z}_{t}}\left(D\left({{X}_{t}}\right)-\overline{D}\right)-V{{W}_{t}}({{X}_{t}}),\\\
\text{s}\text{.t}\text{. }&(\ref{1b}),(\ref{1c}),(\ref{1d}),(\ref{1e}),\\\
&x_{t,k}^{i}\in\left\\{0,1\right\\},\forall i\in\mathcal{I}.\end{split}$ (15)

###  III-B Proposed Algorithm

The transformed problem is a nonlinear integer programming problem, which is
hard to be directly solved. Since the state transition is Markovian, we adopt
the dueling DQN (DDQN) algorithm [15] to obtain the transcoding-path selection
for each GoP in each scheduling slot. Specifically, the dueling architecture
constructs two streams of fully connected layers to provide separate
estimations of value and advantage functions, which can learn which TQ lengths
are valuable without having to learn the effect of transcoding-path selection
of each GoP.

The state includes TQs’ lengths, the deficit queue length, estimated
transcoding workloads, and bit rates of GoPs, i.e.,
${{s}_{t}}=\left\\{{{\left\\{L_{t}^{i}\right\\}}_{i\in\mathcal{I}}},Z_{t},{{\left\\{{{\Omega}_{i}}\left(V_{t}^{k}\right)\right\\}}_{i\in\mathcal{I},k\in\mathcal{K}}},{{\left\\{{{b}_{t,k}}\right\\}}_{k\in\mathcal{K}}}\right\\}$.
The action includes all transcoding decisions in Eq. (15), i.e.,
${{a}_{t}}={{\left\\{x_{t,k}^{i}\right\\}}_{{}_{i\in\mathcal{I},k\in\mathcal{K}}}}$.
The reward at step $t$ is the opposite value of drift-plus-cost, i.e.,
${{r}_{t}}=V{{W}_{t}}({{X}_{t}})-{{Z}_{t}}\left(D\left({{X}_{t}}\right)-\overline{D}\right)$.

In the DDQN algorithm, the target value at step $t$, denoted by
${{y}_{t}^{T}}$, is calculated as follows

$y_{t}^{T}={{r}_{t+1}}+\gamma{{\max}_{{{a}_{t}}}}Q\left({{s}_{t+1}},{{a}_{t}};{{\theta}^{T}}\right),$
(16)

where $\gamma$ is the discount factor, and ${{\theta}^{T}}$ is the network
parameters of the target network. Function $Q\left(\cdot\right)$ is calculated
based on value function $U\left(\cdot\right)$ and advantage function
$A\left(\cdot\right)$, which can be expressed as

$\begin{split}Q\left({{s}_{t}},{{a}_{t}};{{\theta}_{1}},{{\theta}_{2}},{{\theta}_{3}}\right)&=U\left({{s}_{t}};{{\theta}_{1}},{{\theta}_{2}}\right)+A\left({{s}_{t}},{{a}_{t}};{{\theta}_{1}},{{\theta}_{3}}\right)\\\
&-\frac{1}{\left|\mathcal{A}\right|}\sum\nolimits_{a_{t}^{{}^{\prime}}\in\mathcal{A}}{A\left({{s}_{t}},a_{t}^{{}^{\prime}};{{\theta}_{1}},{{\theta}_{3}}\right)},\end{split}$
(17)

where ${{\theta}_{1}}$ indicates the parameters of convolutional layers. Here,
${{\theta}_{2}}$ and ${{\theta}_{3}}$ refer to the parameters of fully-
connected layers for the advantage function and the value function,
respectively. Here, $\mathcal{A}$ is the set of actions. The proposed DT-DDQN
algorithm is presented in Algorithm 1.

1 Initialize the primary network, the target network and the replay memory.

2for  _each episode_ do

3 Reset GoP requests, and CTDT and ETDT queues;

4 for  _each step $t\in\\{1,...,t_{max}\\}$_ do

5

6 Estimate the transcoding workload for each GoP in CTDT and ETDT queues;

7 Update ETDT and CTDT queue lengths based on action, where
${{a}_{t}}=\arg{{\max}_{{{a}_{t}}\in\mathcal{A}}}Q\left({{s}_{t}},{{a}_{t}};\theta\right)$;

8 Obtain the reward ${{r}_{t}}$ and update the state ${{s}_{t+1}}$;

9 Store $\left\\{{{s}_{t}},{{a}_{t}},{{r}_{t}},{{s}_{t+1}}\right\\}$ into the
replay memory, and sample a random mini-batch;

10 Calculate the target value by Eq. (16);

11 Calculate the loss value by
$\mathcal{L}\left(\theta\right)=\mathbb{E}\\!\left[\\!{{\left({{r}_{t}}\\!+\\!\gamma{{\max}_{a}}Q\left({{s}_{t+1}},{{a}_{t+1}};{{\theta}^{T}}\right)\\!-\\!Q\left({{s}_{t}},{{a}_{t}};\theta\right)\right)}^{2}}\right]$;

12 Update the parameters with Adam optimizer;

13 end for

14

15 end for

Algorithm 1 DT-DDQN

##  IV Performance Evaluation

In this section, simulation results are presented to evaluate the performance
of the proposed DT-DDQN algorithm. We adopt a video dataset including twenty
1080p videos and eight 720p videos.333 https://media.xiph.org/video/derf/ The
Matlab neural fitting toolbox with Bayesian
regularization444https://www.mathworks.com/help/deeplearning/gs/fit-data-with-
a-neural-network.html is utilized to generate the BNN-based DT model for
transcoding workload estimation. All input variables in Table I are normalized
to values with a mean of 0 and a standard deviation of 1. In addition, users’
request statistics are sampled from the public dataset [16]. The main
simulation parameters are presented in Table II.

TABLE II: Simulation Parameters Parameter | Value | Parameter | Value
---|---|---|---
Primary DQN learning rate | $10^{-4}$ | $\gamma$ | 0.99
Target DQN learning rate | $10^{-3}$ | $\mathcal{R}$ | 5000
Number of episodes | 1000 | $\overline{D}$ | 1.8 s
Number of steps | 100 | $\kappa$ | [1, 10] MFlops
$L_{i}^{\max},I_{max}$ | 1.5, 0.5 s | $f$ | [5, 20] GHz

We compare the proposed DT-DDQN algorithm with the following benchmark
schemes.

  * $\bullet$

Round Robin (RR): GoPs are placed to TQs in equal portions and in a circular
order.

  * $\bullet$

Proportional Fair (PF): GoPs are sequentially placed in TQs based on
scheduling priority considering GoP requests and transcoding workloads.

  * $\bullet$

Utility-based Multi-dimensional and Multi-choice KnaPsack (UMMKP) [17]: Each
TQ is seen as a knapsack, while each GoP is viewed as a good. Each GoP is
sequentially placed to TQs based on the user satisfaction-based utility value.

  * $\bullet$

DDQN: Only DDQN is adopted for the transcoding-path selection. The TWE model
is a universal one without DT.

Figure 3: BNN-based DT TWE model performance: (a) error histogram, and (b)
MSE.

In Figure 3, 7976 and 1994 instances (samples) are utilized for training and
testing the BNN-based TWE models in DTs, respectively. It can be observed that
the estimation errors are concentrated in [-0.0163, 0.02057], and more than
half instances have fairly small errors. In Fig. 3, the MSE can achieve a fast
degradation in the initial stage and gradually converges to $10^{-4}$ when up
to 1000 episodes. Therefore, the BNN-based TWE model can well extract
historical transcoding data for accurate TWE.

Figure 4: Transcoding performance comparison: (a) cloud and edge TQ lengths,
and (b) user satisfaction.

Figure 4 shows the average lengths and variance of five TQs in the cloud and
edge servers. The TQ lengths under the proposed DT-DDQN algorithm are
relatively balanced with smaller fluctuations compared with the DDQN
algorithm, and the length of TQ 2 is higher. This is because the proposed DT-
DDQN algorithm can accurately estimate the transcoding workload of each GoP to
help the network controller balance transcoding workloads among TQs, and TQ 2
relieves a part of transcoding workloads of TQ 4. Fig. 4 shows the proposed
DT-DDQN algorithm can achieve the highest user satisfaction in the most of
time. Around $500\,s$, the user satisfaction of all algorithms suffers a
distinct degradation since the GoP requests of high-quality version
dramatically increase, and the slow TQ cannot support soaring trancoding
operations of GoPs and have to transfer a part of GoPs to medium TQ and fast
TQ to satisfy the service delay requirement.

##  V Conclusion

In this paper, we have studied a collaborative transcoding problem for better
user satisfaction in live streaming. We have proposed a DT-assisted
collaborative transcoding model to capture the dynamics of cloud and edge TQs
and developed a DRL-based transcoding scheduling algorithm to enhance long-
term user satisfaction within an average service delay threshold. The proposed
DT-assisted collaborative transcoding model can also be applied to distributed
computation offloading to improve the efficiency of cooperative computation.
For the future work, we will investigate how to efficiently coordinate the DT-
assisted cloud-edge collaborative transcoding process and wireless
transmission to further enhance user satisfaction.

## References

  * [1] Grand View Research, “Global video streaming market share report, 2022-2030,” pp. 1–200, 2021.
  * [2] R. Aparicio-Pardo, K. Pires, A. Blanc, and G. Simon, “Transcoding live adaptive video streams at a massive scale in the cloud,” in _Proc. ACM MMSys_ , 2015, Portland, OR, USA, pp. 49–60.
  * [3] Y. Shi, K. Yang, Z. Yang, and Y. Zhou, “Mobile edge artificial intelligence: Opportunities and challenges,” _Amsterdam, The Netherlands: Elsevier_ , pp. 57–65, 2021.
  * [4] Y. Zheng, D. Wu, Y. Ke, C. Yang, M. Chen, and G. Zhang, “Online cloud transcoding and distribution for crowdsourced live game video streaming,” _IEEE Trans. Circuits Syst. Video Technol._ , vol. 27, no. 8, pp. 1777–1789, 2016.
  * [5] Y. Zhu, Q. He, J. Liu, B. Li, and Y. Hu, “When crowd meets big video data: Cloud-edge collaborative transcoding for personal livecast,” _IEEE Trans. Netw. Sci. Eng._ , vol. 7, no. 1, pp. 42–53, 2020.
  * [6] Y. Shi, J. Dong, and J. Zhang, “Low-overhead communications in IoT networks,” _Singapore: Springer_ , pp. 1–10, 2020.
  * [7] H. Pang, C. Zhang, F. Wang, H. Hu, Z. Wang, J. Liu, and L. Sun, “Optimizing personalized interaction experience in crowd-interactive livecast: a cloud-edge approach,” in _Proc. ACM MM_ , 2018, Seoul, Korea, pp. 1217–1225.
  * [8] A. Erfanian, F. Tashtarian, A. Zabrovskiy, C. Timmerer, and H. Hellwagner, “OSCAR: On optimizing resource utilization in live video streaming,” _IEEE Trans. Netw. Serv. Manage._ , vol. 18, no. 1, pp. 552–569, 2021.
  * [9] X. Shen, J. Gao, W. Wu, M. Li, C. Zhou, and W. Zhang, “Holistic network virtualization and pervasive network intelligence for 6G,” _IEEE Commun. Surveys Tuts._ , vol. 24, no. 1, pp. 1–30, 2022.
  * [10] X. Huang, C. Zhou, W. Wu, M. Li, and H. Wu, “Personalized QoE enhancement for adaptive video streaming: A digital twin-assisted scheme,” _arxiv preprint arXiv: 2205.04014_ , 2022.
  * [11] ITU-T, “Subjective video quality assessment methods for multimedia applications,” pp. 1–50, 2022. [Online]. Available: https://www.itu.int/rec/T-REC-P.910-202207-I/en
  * [12] F. Burden and D. Winkler, “Bayesian regularization of neural networks,” _Artificial Neural Netw._ , pp. 23–42, 2018.
  * [13] W. Wu, C. Zhou, M. Li, H. Wu, H. Zhou, N. Zhang, X. Shen, and W. Zhuang, “AI-native network slicing for 6G networks,” _IEEE Wireless Commun._ , vol. 29, no. 1, pp. 96–103, 2022.
  * [14] M. J. Neely, “Stochastic network optimization with application to communication and queueing systems,” _Synthesis Lectures on Communication Networks_ , vol. 3, no. 1, pp. 1–211, 2010.
  * [15] Z. Wang, T. Schaul, M. Hessel, H. Hasselt, M. Lanctot, and N. Freitas, “Dueling network architectures for deep reinforcement learning,” in _Proc. ICML_ , 2016, New York, NY, USA, pp. 1995–2003.
  * [16] F. Loh, F. Wamser, F. Poignée, S. Geißler, and T. Hoßfeld, “Youtube dataset on mobile streaming for internet traffic modeling and streaming analysis,” _Scientific Data_ , vol. 9, no. 1, pp. 1–12, 2022.
  * [17] J. Liu, W. Zhang, S. Huang, H. Du, and Q. Zheng, “QoE-driven HAS live video channel placement in the media cloud,” _IEEE Trans. Multimedia_ , vol. 23, pp. 1530–1541, 2020.