# Differentiable Trajectory Generation for Car-like Robots with Interpolating
Radial Basis Function Networks

Hongrui Zheng, Rahul Mangharam All authors are with the University of
Pennsylvania, Department of Electrical and Systems Engineering, 19104,
Philadelphia, PA, USA. Emails: {hongruiz<EMAIL_ADDRESS>
###### Abstract

The design of Autonomous Vehicle software has largely followed the Sense-Plan-
Act model. Traditional modular AV stacks develop perception, planning, and
control software separately with little integration when optimizing for
different objectives. On the other hand, end-to-end methods usually lack the
principle provided by model-based white-box planning and control strategies.
We propose a computationally efficient method for approximating closed-form
trajectory generation with interpolating Radial Basis Function Networks to
create a middle ground between the two approaches. The approach creates smooth
approximations of local Lipschitz continuous maps of feasible solutions to
parametric optimization problems. We show that this differentiable
approximation is efficient to compute and allows for tighter integration with
perception and control algorithms when used as the planning strategy.

##  I Introduction

Traditionally, the motion planning task for a car-like robot requires
synthesizing trajectories online. The local planning task ultimately searches
for a sequence of feasible control inputs given a desired local goal. This is
naturally expressed as a parametric optimization problem. One can easily
enforce constraints from vehicle dynamics or limits of operation range.
However, even with efficient optimization solvers, solving potentially
hundreds of optimization with different specifications online up to thirty
times a second is still challenging. Massive look-up tables storing
discretized solutions found by running the optimizations offline have been
used to speed up the process online. However, the look-up table can only
provide discrete approximations of the actual optimal control. Moreover, the
trajectory generation remains single-threaded and requires high memory usage
without specialized software implementation.

More recently, different efforts have tried to bypass solving optimal control
online by creating end-to-end planners attempting to directly produce control
inputs using sensor information [1, 2, 3]. End-to-end approaches provide the
benefit of exposing gradient information for upstream and downstream processes
(e.g., a neural network-based perception pipeline). Although these approaches
show potential in efficiently generating local motion plans, they cannot
enforce dynamic constraints without external help or guarantee the solutions’
validity. We propose Interpolating Radial Basis Function Networks (IRBFN), a
differentiable trajectory generation method that produces dynamically feasible
trajectories and can be efficiently parallelized. Unlike existing
differentiable planners, the gradient information is available throughout the
process from local goal to final states on the trajectory. Our approach also
leverages the GPU for highly parallelizable computations for efficiency. One
key contribution of our work is that it preserves dynamic constraints and
approximates solutions from optimal control problems arbitrarily well if
enough training samples are provided. Another key contribution compared to
existing work is that the planner is differentiable with respect to all
possible parameters, e.g. gradients can flow from control or planning loss to
the local goal selection. The scope of this paper only includes providing a
theoretical contribution that enables the use of differentiable planners and
provides an important error bound. As shown in Figure 1, differentiable
planners blend the explainability of traditional modular planners and the
scalability of end-to-end planners.

In the following sections, we’ll present preliminary information on trajectory
generation and Radial Basis Function Networks in Section III-A, define the
interpolating RBFNs in Section III-C, discuss the error bounds of
interpolation in Section III-D, and show benchmarks in Section IV.

Figure 1: Comparison between standard modular planner stack, differentiable
planner stack, and end-to-end planner stack. Gradient information is not
available between modules of a traditional modular planner, and is available
in differentiable and end-to-end planners.

##  II Related Work

###  II-A State Lattice Planners

State lattice based motion planners based on clothoids [4] have been used on
Autonomous Vehicles since the DARPA Urban Challenge [5] and continue to see
success in highly unstructured and competitive environments as well [6]. It
has also been shown that they can be easily integrated as differentiable
planners [7, 8]. The sampling-based scheme that state lattice planners use is
flexible for planning cost evaluation and quick trajectory optimization.

###  II-B Data-driven Motion Planning

One line of work advocates the use of a fully end-to-end trainable Autonomous
driving stack and learning a network-based policy directly from a large amount
of data [1, 2, 3]. These approaches boast the scalability that comes with
using black-box algorithms. However, fully end-to-end approaches lack the
explainability and interpretability of more principled and traditional motion
planning approaches. In addition, verifiability and safety guarantees become
extremely hard to inject into these algorithms due to the use of black-box
approaches. In contrast, another line of work aims to find a middle ground
between end-to-end pipelines and traditional motion planning [8, 9, 10, 11].
These approaches try to preserve interpretability in planning by using
objectives provided by the prediction and detection modules. And use a joint
backbone network that takes sensor observations directly to decision-making to
utilize an end-to-end model.

###  II-C Differentiable Modeling and Motion Planning

More recently, approaches have focused more on applying differentiable
algorithms to Autonomous driving. The idea is to make existing model-based
algorithms differentiable. Differentiable modeling and simulation [12, 13, 14,
15, 16, 17] has been shown to help achieve a better synthesis of controllers
and improve sample efficiency using Reinforcement Learning and other learning-
based methods. Differentiable planning has shown promise in general planning
task [18, 19, 20, 21] as well as vision-based planning tasks and planning
under uncertainty [22, 23, 24]. Lastly, differentiable control [25, 26, 27,
28] has also shown performance on par with traditional control strategies and
provides additional gradient information that could help improve upstream
modules in the software stack. Our approach aligns with this theme closest by
making model-based algorithms differentiable. IRBFNs preserve the
interpretability of planning with state lattice planners and cubic spirals. It
uses a differentiable function approximator that can get arbitrarily close to
the ground truth to make the planning pipeline differentiable. Compared to
approaches similar to [7], the gradient information is available for all
planner parameters, e.g. gradients can flow from control or planning loss to
the local goal selection policy. This enables true end-to-end training of the
planner without sacrificing the properties of model-based algorithms.

###  II-D Function Approximator in Dynamics

The method used in this paper is also related to methods proposed by [29],
where ANNs are used to create differentiable function approximations, and the
approximation properties are studied. Similarly, [30] describes methods for
characterizing non-linear plant models and controllers for these models in
terms of RBFNs. However, neither of these provide an approximation error
bound. Additionally, the training dataset in these works does not preserve the
exact fitting of training data.

##  III Methodology

###  III-A Preliminaries

####  III-A1 Trajectory Generation

We use clothoids as the parameterization of dynamically feasible trajectory
for Ackermann steering vehicles because clothoids are posture continuous and
can be represented as a polynomial. We represent the curvature of the
trajectory as a cubic polynomial of the arc length $s$:

$\kappa(s)=a+bs+cs^{2}+ds^{3}$ (1)

Following [31], we re-formulate the above cubic polynomial such that the
parameters are the curvatures of four equidistant points along the trajectory:

$\begin{split}a&=\kappa_{0}\\\
b&=-\frac{1}{2}\frac{-2\kappa_{3}+11\kappa_{0}-18\kappa_{1}+9\kappa_{2}}{s_{f}-s_{0}}\\\
c&=\frac{9}{2}\frac{-\kappa_{3}+2\kappa_{0}-5\kappa_{1}+4\kappa_{2}}{(s_{f}-s_{0})^{2}}\\\
d&=-\frac{9}{2}\frac{-\kappa_{3}+\kappa_{0}-3\kappa_{1}+3\kappa_{2}}{(s_{f}-s_{0})^{3}}\end{split}$
(2)

To incorporate dynamic constraints into the generation process, we use the
following kinematic dynamic model of the vehicle:

$\begin{split}\mathbf{x}=\left[x,y,\theta,\kappa\right]&~{}~{}~{}\mathbf{u}=\left[v,\dot{\delta}\right]\\\
\dot{x}(t)&=v(t)\cos\theta(t)\\\ \dot{y}(t)&=v(t)\sin\theta(t)\\\
\dot{\theta}(t)&=\kappa(t)v(t)\\\
\dot{\kappa}(t)&=\frac{\dot{\delta}(t)}{L}\end{split}$ (3)

where $\mathbf{x}$ is the state vector consisting of the pose of the vehicle
and the current curvature, $\mathbf{u}$ is the input vector consisting of the
vehicle’s velocity and steering velocity input, and $L$ is the wheelbase of
the vehicle. At low steering angles, we can approximate the curvature of the
vehicle with the steering angle:
$\kappa=\frac{\tan\delta}{L}\approx\frac{\delta}{L}$. Additionally, we
introduce constraints on the initial and final states of the trajectory:

$\begin{split}\mathbf{x}(s_{0})&=\left[x_{0},y_{0},\theta_{0},\kappa_{0}\right]\\\
\mathbf{x}(s_{f})&=\left[x_{g},y_{g},\theta_{g},\kappa_{g}\right]\end{split}$
(4)

The initial state constraint is trivial since it’s determined by the current
state of the vehicle. The final state is determined by the goal pose and goal
curvature of the vehicle. Next, we rewrite the formula using the substitution
$v=\frac{ds}{dt}$ to find the ODEs with respect to the arc length. By dividing
$v$ on both sides, the first three terms in the dynamics then become the
following equations.

$\begin{split}\frac{dx}{ds}=\cos\theta(s),~{}~{}~{}~{}~{}~{}&x(s)=\int\cos\theta(s)\\\
\frac{dy}{ds}=\sin\theta(s),~{}~{}~{}~{}~{}~{}&y(s)=\int\sin\theta(s)\\\
\frac{d\theta}{ds}=\kappa(s),~{}~{}~{}~{}~{}~{}&\theta(s)=\int\kappa(s)ds\end{split}$
(5)

Then by using Equation 1 as the curvature, we perform the following
integrations to find the states on the trajectory.

$\begin{split}\kappa(s)&=a+bs+cs^{2}+ds^{3}\\\
\theta(s)&=as+bs^{2}+cs^{3}+ds^{4}\\\
x(s)&=\int_{0}^{s_{f}}\cos\left(as+\frac{bs^{2}}{2}+\frac{cs^{3}}{3}+\frac{ds^{4}}{4}\right)ds\\\
y(s)&=\int_{0}^{s_{f}}\sin\left(as+\frac{bs^{2}}{2}+\frac{cs^{3}}{3}+\frac{ds^{4}}{4}\right)ds\end{split}$
(6)

We can then set up an optimization where the objective is to minimize the
Euclidean distance between the goal pose $(x_{g},y_{g},\theta_{g})$ of the
trajectory and the integrated pose at the final arc length:

$\text{minimize}~{}~{}~{}~{}~{}|x(s_{f})-x_{g}|^{2}+|y(s_{f})-y_{g}|^{2}+|\theta(s_{f})-\theta_{g}|^{2}$
(7)

Along with Equation 6 as the constraints for optimization, we can find the
optimization variable $q=[\kappa_{0},\kappa_{1},\kappa_{2},\kappa_{3},s_{f}]$.
Following [4], the position quadrature gradients and Hessians can be
efficiently calculated with Simpson’s rule. Thus Newton’s method can be used
for optimization. Alternatively, Powell’s method [32] could also be used to
find the solution without the use of gradients. To enforce dynamic
constraints, we clip the curvature allowed in Equation 6 to the actual
physical limits.

####  III-A2 Look-up Tables

Using the optimization outlined in the previous section, a grid of local goals
in the car’s frame can create a look-up table (LUT) that stores the optimized
parameters. However, efficient online planning using the LUT requires a high
resolution of the look-up grid. In addition, points between the grid points
can’t be interpolated accurately using the stored values.

####  III-A3 Radial Basis Function Networks:

Radial Basis Function Networks (RBFNs) [33] use a smooth function of the
distance of an input to an origin in place of a sigmoidal function as a neuron
in sigmoidal neural networks. We define RBFNs following the standard
definition as follows.

######  Definition III.1 (Radial Basis Function Networks).

An RBFN consists of two layers, a hidden layer with multiple RBF neurons and a
linear layer. Each of the RBF neurons is centered around a predefined or
trainable center where the distances are calculated.

$\Phi(\mathbf{x})=\sum_{i=1}^{M}k_{i}\rho(||\mathbf{x}-c_{i}||)$ (8)

$M$ is the number of centers, $c_{i}$s are the centers for the hidden RBF
layers, $k_{i}$ are the weights for the linear layer. $\rho$ is the smooth
activation function chosen as the Radial Basis Function. In our use case, we
use an inverse quadratic function as the kernel function:
$\rho(\mathbf{z})=\frac{1}{1+z^{2}}$. Usually, during training, the centers of
the RBFN are chosen as the available data points in the training dataset for
better approximations. We left the centers as trainable parameters during our
training process.

###  III-B Problem Definition

Given a local goal $\left[x_{g},y_{g},\theta_{g},\kappa_{g}\right]$ for a car-
like robot, we define the trajectory generation task as creating a sequence of
feasible control inputs $\left\\{\mathbf{a},\mathbf{\delta}\right\\}$, and
corresponding poses and velocity profiles
$\left\\{\mathbf{x,y,\theta,\kappa}\right\\}$ in the workspace, that takes the
robot from its current pose to the local goal. We further formulate a
parametric optimization problem to describe the generation process. From
Section III-A1, we can fully describe a single trajectory with the polynomial
parameters $\left[\kappa_{0},\kappa_{1},\kappa_{2},\kappa_{3},s_{f}\right]$.
Thus the parametric optimization problem is:

$\begin{split}\text{minimize}&~{}~{}~{}|x(s_{f})-x_{g}|^{2}+|y(s_{f})-y_{g}|^{2}+|\theta(s_{f})-\theta_{g}|^{2}\\\
\text{subject
to}&~{}~{}~{}x(s)=\int_{0}^{s_{f}}\cos\left(as+\frac{bs^{2}}{2}+\frac{cs^{3}}{3}+\frac{ds^{4}}{4}\right)ds\\\
&~{}~{}~{}y(s)=\int_{0}^{s_{f}}\sin\left(as+\frac{bs^{2}}{2}+\frac{cs^{3}}{3}+\frac{ds^{4}}{4}\right)ds\\\
&~{}~{}~{}\theta(s)=as+bs^{2}+cs^{3}+ds^{4}\\\ &~{}~{}~{}a=\kappa_{0}\\\
&~{}~{}~{}b=-\frac{1}{2}\frac{-2\kappa_{3}+11\kappa_{0}-18\kappa_{1}+9\kappa_{2}}{s_{f}-s_{0}}\\\
&~{}~{}~{}c=\frac{9}{2}\frac{-\kappa_{3}+2\kappa_{0}-5\kappa_{1}+4\kappa_{2}}{(s_{f}-s_{0})^{2}}\\\
&~{}~{}~{}d=-\frac{9}{2}\frac{-\kappa_{3}+\kappa_{0}-3\kappa_{1}+3\kappa_{2}}{(s_{f}-s_{0})^{3}}\end{split}$
(9)

Then, ultimately, we frame it as a function approximation problem. Where we
approximate a function $f_{\operatorname{opt}}$ such that:

$f_{\operatorname{opt}}\left(\left[x_{g},y_{g},\theta_{g},\kappa_{g}\right]\right)=\left[\kappa^{*}_{0},\kappa^{*}_{1},\kappa^{*}_{2},\kappa^{*}_{3},s^{*}_{f}\right]$
(10)

Where the function being approximated takes in the local goal and outputs the
optimized parameters that describe the clothoid. In the next section, we’ll
introduce our proposed function estimator.

###  III-C Interpolating Radial Basis Function Networks

We introduce a modified RBFN that performs interpolation over a uniform finite
grid of the domain. Together with a smooth indicator function, interpolating
approximations from multiple RBFNs produces bounded approximation error. The
interpolating RBFN consists of multiple RBFNs where each RBFN performs
approximation over an orthotope partitioned by a uniform grid over the domain
of $f_{\operatorname{opt}}$ (shown in Figure 2). We define a smooth and
differentiable indicator function that returns a scalar between 0 and 1 on
each region in each dimension of the input vector.

Figure 2: Network architecture of IRBFNs in the trajectory generation use
case. The network takes a set of local goals as input and outputs parameters
of polynomials describing the desired clothoids.

######  Definition III.2 (Smooth Indicator Function).

A smooth indicator function is a function
$\gamma_{m}:\mathbb{R}^{N}\rightarrow[0,1]$ defined on the orthotope
$r_{m}=[l_{1},u_{1}]\times\cdots\times[l_{N},u_{N}]\in\mathcal{X}$.

$\begin{split}\gamma_{n}(\mathbf{x})=\prod_{m=1}^{M}\left(\frac{\sigma\left(\zeta\left(\mathbf{u}_{m,n}-x_{m}\right)\right)+1}{2}\right)\\\
\times\left(\frac{\sigma\left(\zeta\left(x_{m}-\mathbf{l}_{m,n}\right)\right)+1}{2}\right)\end{split}$
(11)

Where $M$ is the dimensions of the input vector, $N$ is the number of
orthotopes defined, $\sigma$ is the $\mathtt{tanh}$ activation function,
$\zeta$ is a scalar parameter, and $\mathbf{u}$ and $\mathbf{l}$ are the
bounds of the intervals that define the orthotopes.

The bounds of the intervals are defined to be multiples of $\mathbf{s}$, which
is the spacing of the partitioning grid over $\mathcal{X}$. In addition, we
define a controllable parameter $\delta\in[0,1]$ as a function of $N$ and
$\zeta$. An important property of the indicator function emerges by carefully
choosing the values of $\delta$ and $\zeta$. At any point
$\mathbf{x}\in\mathcal{X}$, the value of $\sum_{n}\gamma_{n}(\mathbf{x})$ can
be made arbitrarily close to 1, thus creating smooth interpolation over all
the orthotopes. Intuitively, the indicator function is designed such that the
output of the interpolating RBFNs at the boundaries of the orthotopes in each
dimension is from exactly one-half of the outputs of each neighboring RBF. The
interpolation is performed before the final linear layer in the network. Thus,
we define the Interpolating Radial Basis Function Network (IRBFN) as:

######  Definition III.3 (Interpolating Radial Basis Function Network
(IRBFN)).

$\Phi_{\mathrm{interp}}(\mathbf{x})=\sum_{i=1}^{N}k_{i}\rho(||\mathbf{x}-c_{i}||)\gamma_{i}(\mathbf{x})$
(12)

###  III-D Bounded Interpolation Error

An important property of the interpolating RBFNs is bounded interpolation
error. In the following section, we provide proof with the following sketch:
we first show that $f_{\operatorname{opt}}$ is a computable function, hence
continuous. Then we show that approximation of continuous function using
interpolating Radial Basis Functions on a uniform grid provides a bounded
error. First, we provide some important definitions following [34].

######  Definition III.4 (Oracle Turing Machine).

An Oracle Turing Machine (TM) is an ordinary Turing Machine $M$ equipped with
an additional query tape and two additional states: the query state and the
answer state. When the machine enters the query state, the oracle, a function
$\phi$, replaces the current string $s$ in the query tape by the string
$\phi(s)$, moves the tape head back to the first cell of the query tape, and
puts the machine $M$ in the answer state.

######  Definition III.5 ($k$-oracle Turing Machine).

A $k$-oracle TM $M$ is an Oracle Turing Machine which uses $k$ oracle
functions $\phi_{1},\ldots,\phi_{k}$. M can make queries to the $i$-th
function $\phi_{i},1\leq i\leq k$ by writing down $\langle i,n\rangle$ on its
query tape and the $i$-th oracle will answer by writing $\phi_{i}(n)$ on the
tape.

######  Definition III.6 (Computability of functions).

A real function $f:[0,1]^{k}\rightarrow R$ is computable if there exists a
$k$-oracle TM $M$ such that for all $x_{1},\ldots,x_{k}\in[0,1]$ and all
$\phi_{i}\in CF_{x_{i}},1\leq i\leq k$, $M^{\phi_{1}},\ldots,M^{\phi_{k}}$
halts and outputs a dyadic rational $d$ such that
$|d-f(x_{1},\ldots,x_{k})|\leq 2^{-n}$. Where $CF_{x_{i}}$ denotes the Cauchy
function or the set of all functions binary converging to $x_{i}$.

######  Definition III.7 (Local Lipshitz Continuity).

A function $f$ is locally Lipschitz continuous over a bounded domain
$\mathcal{X}$ if for $\mathbf{u,v}\in\mathcal{X}$ there exists a finite $K$
such that

$||f(\mathbf{u})-f(\mathbf{v})||\leq K||\mathbf{u}-\mathbf{v}||$ (13)

Since $f_{\operatorname{opt}}$ generates unique solutions with gradient
descent given the same initial condition in polynomial time, and the solution
exists universally, $f_{\operatorname{opt}}$ is computable. Note that
computable real functions were first formally defined by Grzegorczyk [35]. The
Oracle TMs definitions are equivalent to the original definition. From [35,
36, 34], computable functions preserve Lipschitz continuity. Hence we can
derive the following approximation bound for the interpolating RBFNs following
[37, 38].

######  Theorem III.1 (Bounded Interpolation Error on Finite Uniform Grid).

The maximum interpolation error is uniformly bounded by the following
equation:

$\begin{split}||\Phi_{\mathrm{interp}}(\mathbf{x})-&f_{\operatorname{opt}}(\mathbf{x})||_{\infty}<\frac{1}{N^{\alpha}}\left[L2^{\frac{\alpha}{2}+1}s^{\alpha}\right.\\\
&\left.+2^{\alpha/2}s^{\alpha}||\Phi_{\mathrm{interp}}||_{\infty}+||f_{\operatorname{opt}}||_{\infty}\right]\end{split}$
(14)

$N$ is the number of training samples, $\alpha$ is the Hölder order, $L$ is
the Hölder constant, and $s$ is the spacing between training samples,
equivalent to the grid spacing we’ve defined. Note that when the function is
locally Lipschitz continuous in our case, $\alpha=1$. It is clear that as the
limit of $s$ goes to zero, implying a finer grid defined for the look-up
table, the interpolation error goes to zero. By increasing the number of
samples, the error can be made arbitrarily small.

###  III-E Differentiability

Using a differentiable region indicator function (Equation 11), the
interpolating RBFN is fully differentiable. In the integration step, we use
Autograd [39] to make the gradient available from the sampled local goals to
the states on the trajectories. The availability of gradients through the
trajectory generation pipeline in a white-box model provides benefits over
black-box in Model-based Reinforcement Learning [17]. Moreover, incorporating
neural components significantly improves computational efficiency when
synthesizing controllers for dynamic systems [16].

##  IV Experiments

###  IV-A Software Implementation and Training Dataset

The trajectory generation pipeline is written in JAX [40] and FLAX [41] and
can be found online at https://github.com/hzheng40/irbfn. The implementation
makes use of just-in-time (JIT) compilation to speed up mathematical
calculations. The model utilizes automatic vectorization (vmap) to create the
multi-headed structure of the network. Finally, during the integration step,
the implementation makes use of the Haskell-like type signature scan to
eliminate for-loops with carryovers. The training dataset is generated using
Newton’s method using the Jacobians and Hessians found in [4]. The resolution
of the look-up table is specified in Table I. All points are used as the
training set since overfitting the available data is desired. The
interpolating RBFNs use 100 trainable centers for each RBFN and 880 regions
(RBFNs). The length of the intervals that defines the orthotopes (regions) is
1.0 meter in $x$, 1.6 meters in $y$, and $0.39$ radians in $\theta$. $\zeta$
used in the indicator function are 15 for $x$, 15 for $y$, and 100 for
$\theta$. Finally, the network is trained using the Adam [42] optimizer with a
learning rate of 0.001, MSE loss, and batch size of 2000. At 400 epochs, the
average training loss over the entire dataset is 0.03107.

TABLE I: Training Data (Look-up Table) Resolution Dimension | Value | Resolution
---|---|---
Min $x$ (m) | 1.0 | 0.1
Max $x$ (m) | 10.0
Min $y$ (m) | -6.0 | 0.1
Max $y$ (m) | 6.0
Min $\theta$ (rad) | $-\pi/2$ | 0.1
Max $\theta$ (rad) | $\pi/2$
# Points | 3213142

###  IV-B Benchmarks and Trajectory Generation Errors

We run all benchmarks on a system with an NVIDIA RTX 2070 Super GPU and an AMD
Ryzen 9 3900X CPU. We profile the generation of 500 trajectories with
different goal points. The peak VRAM usage using our approach was 273.58 MiB.
In 1000 evaluations with different random noises at each evaluation added to
the goal points, our approach was able to achieve an update frequency of
230.08 Hz. Compared to optimizing for the polynomial solutions online at 3.25
Hz, our approach is a 70x+ speed up only using a small amount of VRAM. Figure
3 shows example outputs with goals set at $x=5$ meters with various $y$ and
$\theta$ values.

We also measure the average error between generated trajectories’ endpoints
and the given goal points. Across 500 trajectories spanning a region of 2 to 6
meters in $x$, -4 to 4 meters in $y$, and -0.3 to 0.3 radians in $\theta$. We
compare the experimental error and the theoretical error bounds at the
endpoints of trajectories in Table II. The theoretical error bounds are
obtained by using Equation 14, then propagated through the dynamics
integration in Equation 6. We note that while the $x$ and $y$ experimental
average errors are within the theoretical bound, the $\theta$ error is above
the theoretical bound. This could be due to intrinsic sensitivity in the
$\theta$ dimension. Since the absolute values for $\theta$ are much smaller
than those of $x$ and $y$ while having the same resolution in the look-up
table used as training data. Additionally, since our training error is not
zero, the prediction error adds to the interpolation error, which gets
compounded through integration to generate the states on the trajectories.
Another noticeable decay in the accuracy of trajectory generation is as the
desired goal moves towards the edge of the available training data.

TABLE II: Trajectory Generation Endpoint Errors Dimension | Experimental Error | Theoretical Error
---|---|---
$x$ (m) | 0.0264 | 0.0410
$y$ (m) | 0.0365 | 0.0410
$\theta$ (rad) | 0.0110 | 0.0008
Figure 3: Example trajectory output from IRBFN at $x=5$m, and at various $y$
and $\theta$ values.

##  V Limitations and Conclusions

Limitations: One of the limitations is the difficulty of filtering out invalid
training data during the offline generation of the look-up table. In our
experiments, we reject invalid trajectories by comparing the arclength and the
corresponding goal’s $x$ and $y$ coordinates. The percentage of valid
trajectories can be improved by increasing the iteration limits in the
optimization. However, this step has a tradeoff between quality and
computation time. Another limitation is that this work doesn’t show the
possible improvement in planning using the gradient information provided by
the differentiable pipeline. Future can utilize a trainable goal selection
policy and the IRBFNs in an end-to-end pipeline, i.e., model-based RL. Lastly,
calculating the gradient information is noticeably more computationally
intensive than inferencing on the trained network. Future work could benchmark
and improve the efficiency of gradient calculations.

Conclusions: In this paper, we proposed a differentiable trajectory generation
pipeline for car-like robots with interpolating Radial Basis Function
Networks. Though using IRBFNs to implement a complete planner stack is out of
the scope of this paper, our approach provides an important theoretical
contribution that shows success in trajectory generation and a uniformly
bounded interpolation error. In terms of computation efficiency, our
implementation achieves a 70x+ speed up at 230+ Hz when generating 500
trajectories simultaneously compared to existing methods while only using a
small amount of VRAM. In addition, the gradient information of every parameter
is available throughout the pipeline.

## Acknowledgments

We thank Matthew O’Kelly for the initial ideation of the project and Joshua P.
Reddy for their contribution to the initial experiments of the project.

## References

  * [1] M. Bojarski, D. Del Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang, _et al._ , “End to end learning for self-driving cars,” _arXiv preprint arXiv:1604.07316_ , 2016.
  * [2] A. E. Sallab, M. Abdou, E. Perot, and S. Yogamani, “Deep reinforcement learning framework for autonomous driving,” _arXiv preprint arXiv:1704.02532_ , 2017.
  * [3] A. Jain, L. Del Pero, H. Grimmett, and P. Ondruska, “Autonomy 2.0: Why is self-driving always 5 years away?” _arXiv preprint arXiv:2107.08142_ , 2021\.
  * [4] A. Kelly and B. Nagy, “Reactive nonholonomic trajectory generation via parametric optimal control,” _The International Journal of Robotics Research_ , vol. 22, no. 7-8, pp. 583–601, 2003.
  * [5] C. Urmson, J. Anhalt, D. Bagnell, C. Baker, R. Bittner, M. Clark, J. Dolan, D. Duggins, T. Galatali, C. Geyer, and others, “Autonomous driving in urban environments: Boss and the urban challenge,” _Journal of Field Robotics_ , vol. 25, no. 8, pp. 425–466, 2008, publisher: Wiley Online Library.
  * [6] T. Stahl, A. Wischnewski, J. Betz, and M. Lienkamp, “Multilayer graph-based trajectory planning for race vehicles in dynamic scenarios,” in _2019 IEEE Intelligent Transportation Systems Conference (ITSC)_ , 2019, pp. 3149–3154.
  * [7] P. Karkus, B. Ivanovic, S. Mannor, and M. Pavone, “Diffstack: A differentiable and modular control stack for autonomous vehicles,” _arXiv preprint arXiv:2212.06437_ , 2022.
  * [8] W. Zeng, W. Luo, S. Suo, A. Sadat, B. Yang, S. Casas, and R. Urtasun, “End-to-end interpretable neural motion planner,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 8660–8669.
  * [9] W. Zeng, S. Wang, R. Liao, Y. Chen, B. Yang, and R. Urtasun, “Dsdnet: Deep structured self-driving network,” in _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXI 16_. Springer, 2020, pp. 156–172.
  * [10] A. Cui, S. Casas, A. Sadat, R. Liao, and R. Urtasun, “Lookout: Diverse multi-future prediction and planning for self-driving,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 16 107–16 116.
  * [11] S. Casas, A. Sadat, and R. Urtasun, “Mp3: A unified model to map, perceive, predict and plan,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 14 403–14 412.
  * [12] E. Heiden, D. Millard, H. Zhang, and G. S. Sukhatme, “Interactive differentiable simulation,” _arXiv preprint arXiv:1905.10706_ , 2019.
  * [13] J. Degrave, M. Hermans, J. Dambre, _et al._ , “A differentiable physics engine for deep learning in robotics,” _Frontiers in neurorobotics_ , p. 6, 2019.
  * [14] C. Song and A. Boularias, “Learning to slide unknown objects with differentiable physics simulations,” _arXiv preprint arXiv:2005.05456_ , 2020\.
  * [15] F. de Avila Belbute-Peres, K. Smith, K. Allen, J. Tenenbaum, and J. Z. Kolter, “End-to-end differentiable physics for learning and control,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [16] E. Heiden, D. Millard, E. Coumans, Y. Sheng, and G. S. Sukhatme, “Neuralsim: Augmenting differentiable simulators with neural networks,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 9474–9481.
  * [17] M. Lutter, J. Silberbauer, J. Watson, and J. Peters, “Differentiable physics models for real-world offline model-based reinforcement learning,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 4163–4170.
  * [18] A. Tamar, Y. Wu, G. Thomas, S. Levine, and P. Abbeel, “Value iteration networks,” _Advances in neural information processing systems_ , vol. 29, 2016.
  * [19] N. Nardelli, G. Synnaeve, Z. Lin, P. Kohli, P. H. Torr, and N. Usunier, “Value propagation networks,” _arXiv preprint arXiv:1805.11199_ , 2018.
  * [20] L. Lee, E. Parisotto, D. S. Chaplot, E. Xing, and R. Salakhutdinov, “Gated path planning networks,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 2947–2955.
  * [21] D. S. Chaplot, D. Pathak, and J. Malik, “Differentiable spatial planning using transformers,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 1484–1495.
  * [22] P. Karkus, X. Ma, D. Hsu, L. P. Kaelbling, W. S. Lee, and T. Lozano-Pérez, “Differentiable algorithm networks for composable robot learning,” _arXiv preprint arXiv:1905.11602_ , 2019.
  * [23] P. Karkus, D. Hsu, and W. S. Lee, “Qmdp-net: Deep learning for planning under partial observability,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [24] S. Gupta, J. Davidson, S. Levine, R. Sukthankar, and J. Malik, “Cognitive mapping and planning for visual navigation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 2616–2625.
  * [25] M. Okada, L. Rigazio, and T. Aoshima, “Path integral networks: End-to-end differentiable optimal control,” _arXiv preprint arXiv:1706.09597_ , 2017\.
  * [26] B. Amos, I. Jimenez, J. Sacks, B. Boots, and J. Z. Kolter, “Differentiable mpc for end-to-end planning and control,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [27] M. Pereira, D. D. Fan, G. N. An, and E. Theodorou, “Mpc-inspired neural network policies for sequential decision making,” _arXiv preprint arXiv:1802.05803_ , 2018.
  * [28] S. East, M. Gallieri, J. Masci, J. Koutnik, and M. Cannon, “Infinite-horizon differentiable model predictive control,” _arXiv preprint arXiv:2001.02244_ , 2020.
  * [29] S. Ferrari and R. F. Stengel, “Smooth function approximation using neural networks,” _IEEE Transactions on Neural Networks_ , vol. 16, no. 1, pp. 24–38, 2005.
  * [30] M. Pottmann and D. E. Seborg, “A nonlinear predictive control strategy based on radial basis function models,” _Computers & chemical engineering_, vol. 21, no. 9, pp. 965–980, 1997, publisher: Elsevier.
  * [31] T. M. Howard _et al._ , “Adaptive model-predictive motion planning for navigation in complex environments,” 2009.
  * [32] M. J. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” _The computer journal_ , vol. 7, no. 2, pp. 155–162, 1964.
  * [33] D. S. Broomhead and D. Lowe, “Radial basis functions, multi-variable functional interpolation and adaptive networks,” Royal Signals and Radar Establishment Malvern (United Kingdom), Tech. Rep., 1988.
  * [34] K.-I. Ko, _Complexity theory of real functions_. Birkhauser Boston Inc., 1991.
  * [35] A. Grzegorczyk, “Computable functionals,” _Fundamenta Mathematicae_ , vol. 42, no. 19553, pp. 168–202, 1955.
  * [36] ——, “On the definitions of computable real continuous functions,” _Fund. Math._ , vol. 44, pp. 61–71, 1957.
  * [37] D. Costarelli, “Sigmoidal functions approximation and applications,” _Dottorato di Ricera in Matematica, Departiment di Matematica e Fisica, Sezione di Matematica, Roma, UNIVERSITA DEGLI STUDI_ , 2014.
  * [38] A. Bejancu Jr, “Local accuracy for radial basis function interpolation on finite uniform grids,” _Journal of approximation theory_ , vol. 99, no. 2, pp. 242–257, 1999.
  * [39] D. Maclaurin, D. Duvenaud, and R. P. Adams, “Autograd: Effortless gradients in numpy,” in _ICML 2015 AutoML workshop_ , vol. 238, no. 5, 2015.
  * [40] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, “JAX: composable transformations of Python+NumPy programs,” 2018. [Online]. Available: http://github.com/google/jax
  * [41] J. Heek, A. Levskaya, A. Oliver, M. Ritter, B. Rondepierre, A. Steiner, and M. van Zee, “Flax: A neural network library and ecosystem for JAX,” 2023\. [Online]. Available: http://github.com/google/flax
  * [42] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.


# Multi-Loss Convolutional Network with Time-Frequency Attention for Speech
Enhancement

1st Liang Wan  School of Communication and Information Engineering
Chongqing University of Posts and Telecommunications
Chongqing, China
<EMAIL_ADDRESS>2nd Hongqing Liu  School of Communication
and Information Engineering
Chongqing University of Posts and Telecommunications
Chongqing, China
<EMAIL_ADDRESS>3rd Yi Zhou  School of Communication and
Information Engineering
Chongqing University of Posts and Telecommunications
Chongqing, China
<EMAIL_ADDRESS>4rd Jie Jia  vivo Mobile Commun co Ltd
China
<EMAIL_ADDRESS>
###### Abstract

The Dual-Path Convolution Recurrent Network (DPCRN) was proposed to
effectively exploit time-frequency domain information. By combining the DPRNN
module with Convolution Recurrent Network (CRN), the DPCRN obtained a
promising performance in speech separation with a limited model size. In this
paper, we explore self-attention in the DPCRN module and design a model called
Multi-Loss Convolutional Network with Time-Frequency Attention(MNTFA) for
speech enhancement. We use self-attention modules to exploit the long-time
information, where the intra-chunk self-attentions are used to model the
spectrum pattern and the inter-chunk self-attention are used to model the
dependence between consecutive frames. Compared to DPRNN, axial self-attention
greatly reduces the need for memory and computation, which is more suitable
for long sequences of speech signals. In addition, we propose a joint training
method of a multi-resolution STFT loss and a WavLM loss using a pre-trained
WavLM network. Experiments show that with only 0.23M parameters, the proposed
model achieves a better performance than DPCRN.

###### Index Terms:

Speech enhancement, axial self-attention, multiple losses, time-frequency
domain

##  I Introduction

Excessive noise and reverberation can significantly impair the performance of
automatic speech recognition (ASR) systems and reduce speech intelligibility
during communication. To improve speech intelligibility and perceptual
quality, speech enhancement techniques are generally employed to separate
clean speech from background interferences. This separation helps eliminate
the effects of noise and reverberation, resulting in clearer and more
understandable speech.

Deep neural networks (DNNs) have become increasingly popular in audio
processing applications, and have been shown to be highly effective in speech
enhancement tasks. In terms of domain used, DNN-based speech enhancement
methods can be broadly divided into two categories: time-frequency (TF)
domain[1, 2, 3, 4, 5] and time domain[6, 7, 8] methods. Methods operating in
the time domain directly predict clean speech waveforms by end-to-end
training, avoiding the challenge of estimating phase information in the TF
domain. Conv-Tasnet[9], as a representative time-domain method, employs a 1-D
convolutional neural network to encode the waveform into effective
representations for accurate speech estimation, and then recovers the waveform
by a decoder consisting of transposed convolutional layers. However, modeling
extremely long sequences is challenging in the time domain, hence deep
convolutional layers such as wave-u-net[7] are often needed for feature
compression. TF-domain methods operate on the spectrogram, under the
assumption that the detailed structures of speech and noise can be more easily
distinguished and separated through time-frequency representations. In TF-
domain methods, training targets can be broadly categorized into two groups:
masking-based targets and mapping-based targets. Masking-based targets such as
ideal binary mask (IBM)[10], ideal ratio mask (IRM)[2], and spectral magnitude
mask (SMM)[11], only consider the magnitude differences between clean speech
and mixture speech while ignoring phase information. Subsequently, complex
ratio mask (CRM)[12] was proposed as a masking-based target, which enhances
both the real and imaginary components of the division between clean speech
and mixture speech spectrograms, allowing for a better speech reconstruction.

For TF domain method, Dual-Path Convolution Recurrent Network (DPCRN)[13]
combining DPRNN[14] and CRN[15], it is possible to obtain a well-behaved
model. On the basis of DPCRN, in this work, a new model called Multi-Loss
Convolutional Network with Time-Frequency Attention(MNTFA) for speech
enhancement is proposed. It replace DPRNN with an Axial Self-Attention
(ASA)[16] block consisting of two ASAs. ASA improves network ability to
capture long-range relations between features in time domain and frequency
domain, respectively. Compared to DPRNN, ASA greatly reduces the need for
memory and computation, which is more suitable for long sequences of speech
signals. The high-level features extracted by encoder are fed into ASA blocks
for further processing. Subsequently, the decoder reconstructs the low-
resolution features to the output.

While DNN-based speech enhancement has made rapid progress, it still faces
challenges in real-world applications due to low signal-to-noise ratio (SNR),
high reverberation, and far-field pickup. These issues remain significant
obstacles to achieving high-quality speech enhancement in practical settings.
Previous studies have found that blindly feeding the output of a speech
enhancement system (SE) to an ASR system can result in degraded performance.
However, incorporating a pre-trained ASR system during the SE model training
stage has been shown to improve ASR performance. For the proposed model, we
employ a multi-loss strategy formed by SE loss and an ASR loss to
simultaneously ensure good SE and ASR performance. The SE loss consists by a
mean square error (MSE) loss for noise suppression and multi-resolution STFT
loss[17] for speech preservation. The ASR loss[18] is formed by a pretrained
WavLM[19] network.

The proposed system is evaluated using PESQ and STOI for SE performance and
word error rate (WER) for ASR performance. On the test datasets, our model
achieves competitive results as DPCRN while having a smaller number of
parameters and computational complexity.

Figure 1:  Architecture of the proposed MNTFA

##  II Network

###  II-A Problem formulation

The observed noisy speech can be represented as a sum of the clean speech and
noise signals in the time domain, which is $y(t)=s(t)+n(t)$, where $y(t)$,
$s(t)$, and $n(t)$ refer to the noisy, clean, and noise signals, respectively.
The equation can be converted to the time-frequency domain by utilizing the
Short-Time Fourier Transform (STFT), given by

$Y(t,f)=S(t,f)+N(t,f),$ (1)

where $Y(t,f)$, $S(t,f)$, and $N(t,f)$ represent the time-frequency bins of
the noisy, clean, and noise speech spectrograms, respectively, at time frame
$t$ and frequency index $f$. To recover clean speech from a noisy mixture, a
mask $M(t,f)$ is typically estimated and multiplied with the noisy speech
spectrogram $Y(t,f)$[20]. To better utilize phase information, we adopt
complex ratio mask (CRM), denoted by $M(t,f)=M_{r}(t,f)+iM_{i}(t,f)$, where
$M_{r}(t,f)$ and $M_{i}(t,f)$ represent the real and imaginary parts of the
mask. CRM better reconstructs speech by enhancing both real and imaginary
spectrogram[21] simultaneously. The denoising process can be represented as
the complex product of mask and noisy speech in the form of:

$\widetilde{S}(t,f)=Y(t,f)\bigodot M(t,f),$ (2)

where $\bigodot$ denotes element-wise multiplication and $\widetilde{S}(t,f)$
is the enhanced speech.

###  II-B Model architecture

The Dual-Path Convolution Recurrent Network combining DPRNN and CRN achieves a
state-of-the-art (SOTA) performance in single-channel speech separation task.
Similar to the original DPRNN, DPCRN consists of an encoder, a dual-path RNN
module, and a decoder. The structure of the encoder and decoder is similar to
CRN. DPCRN replaces the RNN part of the CRN with the DPRNN modules. It is
believed that modeling the dependence of frequency is beneficial to speech
enhancement due to the harmonic structures of speech. Instead of learning the
dependence in the time domain, the intra-chunk RNNs are applied to model the
spectral patterns in a single frame and inter-chunk RNNs model the time
dependence of a certain frequency.

However, the DPCRN suffers from a high computational complexity. In this
paper, we propose a novel model with an Axial Self-Attention block consisting
of two ASAs to model long-time dependencies. Figure 1 shows the overall
network structure of the proposed model. Axial Self-Attention (ASA) was
introduced in [16] as a more efficient approach that reduces the demand for
memory and computation resources. ASA is particularly useful for processing
lengthy speech signal sequences. The ASA structure, illustrated in Figure 2,
consists of input channels represented by $C_{i}$ and attention channels
represented by $C$. ASA calculates attention scores along the frequency and
time axis, which are referred to as F-attention and T-attention, respectively.
These attention score matrices enhance the model’s ability to capture long-
range dependencies in both the time and frequency domains simultaneously.
Similar to DPCRN, we construct an ASA block by combining two ASAs. Figure 3
shows the structure of ASA block, and as shown in figure 3, the input and
output dimensions of the first ASA are T, F. We rearrange the output
dimensions of the first ASA to F, T, and pass it to the second ASA. In MNTFA,
we use two ASA blocks to capture long-range relations between features. To
alleviate the problem of gradient vanishing, a residual connection[22] was
utilized between the intra-chunk ASA and inter-chunk ASA. We found that this
development enables the model to better learn both temporal and spectral
contextual information simultaneously, which is particularly relevant in the
context of speech processing.

###  II-C Training target

When training, MNTFA estimates CRM and is optimized by signal approximation
(SA)[23]. Given the complex-valued STFT spectrogram of clean speech S and
noisy speech $Y$, CRM is defined as

$CRM=\frac{Y_{r}S_{r}+Y_{i}S_{i}}{{Y_{r}}^{2}+{Y_{i}}^{2}}+j\frac{Y_{r}S_{i}-Y_{i}S_{r}}{{Y_{r}}^{2}+{Y_{i}}^{2}},$
(3)

where $Y_{r}$ and $Y_{i}$ denote the real and imaginary parts of the noisy
complex spectrogram, respectively. The real and imaginary parts of the clean
complex spectrogram are represented by $S_{r}$ and $S_{i}$. By multiplying the
spectrogram of noisy speech $Y=Y_{r}+Y_{i}$ with the estimated mask
$M=M_{r}+M_{i}$, we obtain the enhanced spectrogram in the form of:
$\widetilde{S}=Y_{r}M_{r}–Y_{i}M_{i}+i(Y_{r}M_{i}+X_{i}M_{r})$, which is
converted back to the waveform using inverse STFT (iSTFT), given by

$\widetilde{s}=iSTFT(\widetilde{S}).$ (4)

###  II-D Loss function

We use three loss functions to measure the estimation error from the signal,
perceptual, and ASR aspects, respectively. To suppress the noise for the low-
SNR time-frequency bins, we employ a metric called mean square error (MSE) of
the spectrogram, which is

$\begin{split}L_{MSE}=&\log(MSE(S_{r},\widetilde{S_{r}})+MSE(S_{i},\widetilde{S_{i}})\\\
&+MSE(|S|,|\widetilde{S}|)).\end{split}$ (5)

The MSE loss is composed of three components that individually quantify the
differences in the real part, imaginary part, and magnitude between the
predicted spectrogram and the ground truth.

To enhance the clarity and intelligibility of speech following the denoising
process, we used a multi-resolution STFT auxiliary loss[17]. Similar to prior
research[24], we define a single Short-Time Fourier Transform (STFT) loss as

$L_{s}=L_{sc}(s,\hat{s})+L_{mag}(s,\hat{s}),$ (6)

where $s$, $\hat{s}$ are the clean and estimated time-domain waveform.
$L_{sc}$ and $L_{mag}$ denote spectral convergence and log STFT magnitude
loss, respectively, which are [25]

$L_{sc}(s,\hat{s})=\frac{\lVert\left|STFT(s)\right|-\left|STFT(\hat{s})\right|\rVert_{F}}{\lVert\left|STFT(s)\right|\rVert_{F}},$
(7)

$L_{mag}(s,\hat{s})=\frac{1}{N}\lVert\log\left|STFT(s)\right|-\log\left|STFT(\hat{s})\right|\rVert_{1},$
(8)

where $\lVert.\rVert_{F}$ and $\lVert.\rVert_{1}$ denote the Frobenius and
$L_{1}$ norms, respectively; $|STFT(.)|$ and N denote the STFT magnitudes and
number of elements in the magnitude, respectively. The multi-resolution STFT
loss is obtained by summing the STFT losses computed with different analysis
parameters such as FFT size, window size, and frame shift. If there are $M$
STFT losses, the multi-resolution STFT auxiliary loss ($L_{aux}$) is

$L_{aux}=\frac{1}{M}\sum_{m=1}^{M}L_{s}^{(m)}.$ (9)

When using STFT for time-frequency representation of signals, the balance
between time and frequency resolution is crucial. Increasing the window size
improves the frequency resolution but reducing the temporal resolution. By
using multiple STFT losses with different analysis parameters, the model can
better learn the time-frequency characteristics of speech and avoid becoming
overfit to a fixed STFT representation, which could result in suboptimal
performance in the waveform-domain.

Moreover, an ASR-oriented loss $L_{ASR}$ is used to reduce the speech
distortion of enhanced speech and reduce WER for a given ASR system. More
specifically, KL divergence distance between clean speech ASR embedding and
estimated speech ASR embedding is employed as our ASR loss. The ASR embeddings
are extracted by the pre-trained WavLM.

Our final loss function for the model is defined as a linear combination of
the MSE loss, Multi-resolution STFT auxiliary loss, and ASR loss, with a equal
weighting, given by

$L=L_{MSE}+L_{aux}+L_{ASR},$ (10)

Figure 2: Axial self attention module Figure 3: ASA blocks TABLE I: PESQs and
STOI of different methods.

Metrics PESQ STOI (%) SNR(dB) -5 0 5 Avg. -5 0 5 Avg. Noisy 1.58 1.84 2.14
1.85 65.50 73.25 79.53 72.76 DCCRN 2.30 2.68 3.02 2.67 74.52 81.32 86.53 80.79
DPCRN 2.36 2.67 2.91 2.65 75.08 81.60 86.28 80.98 MNTFA 2.42 2.72 3.02 2.72
75.64 81.97 86.78 81.46

##  III Experiments

###  III-A Datasets

The clean speech and noise datasets from the DNS Challenge (ICASSP 2022)[26]
are used for our experiments. By randomly pairing speech and noise signals, we
synthesize a training set with 500 hours samples and a validation set with 50
hours, in both of which the signal-to-noise ratio (SNR) is randomly sampled
between -5 dB and 5 dB. Note that all speech and noise signals are randomly
truncated to 10 seconds before mixing. We additionally use the synthetic test
set released by DNS Challenge for testing. All the audios used is sampled at
16 kHz.

TABLE II: The model size, MACS and WER of different models. Model | Para.(M) | MACs(GFLOPS) | WER
---|---|---|---
Noisy |  |  | 13.12
DCCRN | 3.67 | 14.38 | 15.23
DPCRN | 0.8 | 8.45 | 13.71
MNTFA | 0.23 | 1.89 | 12.55

###  III-B Result

Multiple metrics are employed to measure the speech enhancement performance,
including perceptual evaluation speech quality (PESQ)[27], short-time
objective intelligibility (STOI)[28] and word error rate (WER). A pretrained
model based on the conformer architectures[29] are used to compute the WER.

We compare MNTFA with DCCRN and DPCRN on the test set. From the comparison
results shown in Table I, our approach produces a superior performance to
DCCRN and DPCRN in terms of PESQ and STOI. In table II, ”Para.”, ”MACs”, and
”WER” respectively represent the parameter amount of the model, computation
complexity, and word error rate for ASR performance. With 0.23 M parameters,
1.89 Giga floating-point operations per second (GFLOPs) and without any future
information, our model outperform DPCRN with 0.8M parameters.

TABLE III: Ablation study on test set. System | Loss | PESQ | STOI | WER
---|---|---|---|---
Noisy | - | 1.85 | 72.76 | 13.13
MNTFA | $L_{MSE}$ | 2.67 | 81.15 | 13.79
MNTFA | $L_{MSE}+L_{ASR}$ | 2.66 | 81.27 | 12.74
MNTFA | $L_{MSE}+L_{aux}$ | 2.77 | 81.64 | 14.35
MNTFA | $L_{MSE}+L_{ASR}+L_{aux}$ | 2.72 | 81.46 | 12.55

At the same time, we conducted ablation experiments on the proposed model to
show the importance of the different losses. Table III shows the ablation
results. After removing ASR loss, PESQ, STOI, WER increased by 0.05, 0.18, and
1.8, respectively. After removing multi-resolution STFT auxiliary loss, PESQ
and STOI decreased by 0.06 and 0.37, respectively, but WER increased 0.19. The
results show that MNTFA with MSE loss and multi-resolution STFT auxiliary loss
improved the performance in PESQ and STOI, but not the WER. We found that
there is no absolute direct correlation between PESQ and WER, as a higher PESQ
score does not necessarily guarantee a lower WER, which means the higher PESQ
should not always be the final goal.

## CONCLUSION

In this paper, we propose a multi-loss network with Time-Frequency attention
for speech enhancement. In the proposed method, the self-attention blocks are
utilized to model the long-time relations, and we exploit the T-attention and
F-attention to model the temporal and frequency decencies. With 0.23 M
parameters, our model achieves competitive results and we also found that
there is no absolute direct correlation between PESQ and WER, as a higher PESQ
score does not necessarily guarantee a lower WER. In the future, we will
investigate a better way to ensure that speech recognition accuracy while also
preserving a good noise reduction performance.

## References

  * [1] S. Srinivasan, N. Roman, and D. Wang, “Binary and ratio time-frequency masks for robust speech recognition,” _Speech Communication_ , vol. 48, no. 11, pp. 1486–1501, 2006, robustness Issues for Conversational Interaction. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0167639306001129
  * [2] A. Narayanan and D. Wang, “Ideal ratio mask estimation using deep neural networks for robust speech recognition,” in _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2013, pp. 7092–7096.
  * [3] Y. Zhao, D. Wang, I. Merks, and T. Zhang, “Dnn-based enhancement of noisy and reverberant speech,” in _2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2016, pp. 6525–6529.
  * [4] Y. Xu, J. Du, L.-R. Dai, and C.-H. Lee, “An experimental study on speech enhancement based on deep neural networks,” _IEEE Signal Processing Letters_ , vol. 21, no. 1, pp. 65–68, 2014.
  * [5] D. Yin, C. Luo, Z. Xiong, and W. Zeng, “Phasen: A phase-and-harmonics-aware speech enhancement network,” _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 05, pp. 9458–9465, Apr. 2020. [Online]. Available: https://ojs.aaai.org/index.php/AAAI/article/view/6489
  * [6] S.-W. Fu, T.-W. Wang, Y. Tsao, X. Lu, and H. Kawai, “End-to-end waveform utterance enhancement for direct evaluation metrics optimization by fully convolutional neural networks,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 26, no. 9, pp. 1570–1584, 2018.
  * [7] S. Daniel, E. Sebastian, and D. Simon, “Wave-u-net: A multi-scale neural network for end-to-end audio source separation,” _arXiv preprint arXiv:1806.03185_ , 2018.
  * [8] L. Zhang, Z. Shi, J. Han, A. Shi, and D. Ma, “Furcanext: End-to-end monaural speech separation with dynamic gated dilated temporal convolutional networks,” in _MultiMedia Modeling: 26th International Conference, MMM 2020, Daejeon, South Korea, January 5–8, 2020, Proceedings, Part I 26_. Springer, 2020, pp. 653–665.
  * [9] Y. Luo and N. Mesgarani, “Conv-tasnet: Surpassing ideal time-frequency magnitude masking for speech separation,” _IEEE Transactions on Audio, Speech, and Language Processing_ , pp. 1–1, 2019.
  * [10] D. Wang, “On ideal binary mask as the computational goal of auditory scene analysis,” _Speech separation by humans and machines_ , pp. 181–197, 2005\.
  * [11] Y. Wang, A. Narayanan, and D. Wang, “On training targets for supervised speech separation,” _IEEE/ACM transactions on audio, speech, and language processing_ , vol. 22, no. 12, pp. 1849–1858, 2014.
  * [12] K. Tan and D. Wang, “Complex spectral mapping with a convolutional recurrent network for monaural speech enhancement,” in _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2019, pp. 6865–6869.
  * [13] X. Le, H. Chen, K. Chen, and J. Lu, “Dpcrn: Dual-path convolution recurrent network for single channel speech enhancement,” _arXiv preprint arXiv:2107.05429_ , 2021.
  * [14] Y. Luo, Z. Chen, and T. Yoshioka, “Dual-path rnn: efficient long sequence modeling for time-domain single-channel speech separation,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 46–50.
  * [15] T. Ke and W. DeLiang, “A convolutional recurrent neural network for real-time speech enhancement,” _Interspeech_ , pp. 3229–3233, 2018.
  * [16] G. Zhang, C. Wang, L. Yu, and J. Wei, “Multi-scale temporal frequency convolutional network with axial attention for multi-channel speech enhancement,” in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 9206–9210.
  * [17] R. Yamamoto, E. Song, and J.-M. Kim, “Parallel wavegan: A fast waveform generation model based on generative adversarial networks with multi-resolution spectrogram,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6199–6203.
  * [18] T.-A. Hsieh, C. Yu, S.-W. Fu, X. Lu, and Y. Tsao, “Improving perceptual quality by phone-fortified perceptual loss using wasserstein distance for speech enhancement,” _arXiv preprint arXiv:2010.15174_ , 2020.
  * [19] S. Chen, C. Wang, Z. Chen, Y. Wu, S. Liu, Z. Chen, J. Li, N. Kanda, T. Yoshioka, X. Xiao _et al._ , “Wavlm: Large-scale self-supervised pre-training for full stack speech processing,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 6, pp. 1505–1518, 2022.
  * [20] D. Wang and J. Chen, “Supervised speech separation based on deep learning: An overview,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 26, no. 10, pp. 1702–1726, 2018.
  * [21] Y. Hu, Y. Liu, S. Lv, M. Xing, S. Zhang, Y. Fu, J. Wu, B. Zhang, and L. Xie, “Dccrn: Deep complex convolution recurrent network for phase-aware speech enhancement,” _arXiv preprint arXiv:2008.00264_ , 2020.
  * [22] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [23] ——, _Delving deep into rectifiers: Surpassing human-level performance on imagenet classification_. Proceedings of the IEEE international conference on computer vision 1026–1034, 2015.
  * [24] R. Yamamoto, E. Song, and J.-M. Kim, “Probability density distillation with generative adversarial networks for high-quality parallel waveform generation,” _arXiv preprint arXiv:1904.04472_ , 2019.
  * [25] S. Ö. Arık, H. Jun, and G. Diamos, “Fast spectrogram inversion using multi-head convolutional neural networks,” _IEEE Signal Processing Letters_ , vol. 26, no. 1, pp. 94–98, 2018.
  * [26] C. K. Reddy, V. Gopal, R. Cutler, E. Beyrami, R. Cheng, H. Dubey, S. Matusevych, R. Aichner, A. Aazami, S. Braun _et al._ , “The interspeech 2020 deep noise suppression challenge: Datasets, subjective testing framework, and challenge results,” _arXiv preprint arXiv:2005.13981_ , 2020.
  * [27] A. W. Rix, J. G. Beerends, M. P. Hollier, and A. P. Hekstra, “Perceptual evaluation of speech quality (pesq)-a new method for speech quality assessment of telephone networks and codecs,” in _2001 IEEE international conference on acoustics, speech, and signal processing. Proceedings (Cat. No. 01CH37221)_ , vol. 2. IEEE, 2001, pp. 749–752.
  * [28] C. H. Taal, R. C. Hendriks, R. Heusdens, and J. Jensen, “An algorithm for intelligibility prediction of time–frequency weighted noisy speech,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 19, no. 7, pp. 2125–2136, 2011.
  * [29] D. Wu, B. Zhang, C. Yang, Z. Peng, W. Xia, X. Chen, and X. Lei, “U2++: Unified two-pass bidirectional end-to-end model for speech recognition,” _arXiv preprint arXiv:2106.05642_ , 2021.


# The Evolution of Flow and Mass Transport in 3D Confined Cavities

Reem Khojah  Bioengineering Department, University of California Los Angeles,
Los Angeles, California 90095, USA  Darren Lo  Bioengineering Department,
University of California Los Angeles, Los Angeles, California 90095, USA
Fiona Tang  Bioengineering Department, University of California Los Angeles,
Los Angeles, California 90095, USA  Dino Di Carlo<EMAIL_ADDRESS>Bioengineering Department, University of California Los Angeles, Los Angeles,
California 90095, USA

###### Abstract

Flow in channels and ducts with adjoining cavities are common in natural and
engineered systems. Here we report numerical and experimental results of 3D
confined cavity flow, identifying critical conditions in the recirculating
flow formation and mass transport over a range of channel flow properties
($0.1\leq Re\leq 300$) and cavity aspect ratio ($0.1\leq H/X_{s}\leq 1$). In
contrast to 2D systems, a mass flux boundary is not formed in 3D confined
cavity-channel flow. Streamlines directly enter a recirculating vortex in the
cavity and exit to the main channel leading to an exponential increase in the
cavity mass flux when the recirculating vortex fills the cavity volume. These
findings extend our understanding of flow entry and exit in cavities and
suggest conditions where convective mass transport into and out of cavities
would be amplified and vortex particle capture reduced.

Rare cells Microfiltration Tangential-flow Inertial separation

###### pacs:

Valid PACS appear here

When flowing fluid with finite inertia in a channel that suddenly expands in
cross-sectional dimension, an adverse pressure gradient can occur resulting in
the formation of recirculating flow in the expansion region Batchelor (2000);
Macagno and Hung (1967); Baloch  _et al._ (1995). If the cross-sectional
dimension is returned to the initial value further downstream, this creates a
cavity, which can support a recirculating flow up to the cavity size Sinha
_et al._ (1982). The formation of recirculating flows in 3D-confined cavities
can result in selective particle capture from the mainstream channel flow, an
area attracting significant interest in biotechnology and biomedicine where
microcavities have been utilized as well-controlled reaction chambers for
trapped bio-particles in miniaturized diagnostic devices Marcus  _et al._
(2006); Gierahn  _et al._ (2017); Hur  _et al._ (2011); Khojah  _et al._
(2017); Dhar  _et al._ (2018). Moreover, cavity flow capture is a diverse
phenomenon which can also be found in nature such as in the assembly of
biological self-replicating molecules inside microvortices that arise in
porous ocean rocks Sun  _et al._ (2018), as well as platelet aggregates
accumulation in aneurysms Rayz  _et al._ (2008) (Figs 1(a-c)). Thus, a better
understanding of recirculating cavity flow can elucidate trapping phenomena in
natural and physiological flows and enable the engineering of mass transport
in cavities and microwells.

Figure 1:  Fluid mass transport in three-dimensional cavity flow geometries
commonly found in (a) porous ocean rocks, (b) physiological blood flow in
aneurysms, and (c) selective cell trapping in microchambers. Confocal 3D
imaging of fluorophore containing streams shows flow enters and exits a cavity
without circulating at low Reynolds numbers ($Re$) in (d) and fully-developed
circulating cavity flow at higher $Re$ in (e) reveals a set of streamlines
swirling out of the vortex core and returing to the main channel flow (Fig.
S7, video S2).  Figure 2:  Analysis of recirculating vortex flow development
in three-dimensional confined cavities (a) Three dimensional flow diagram in a
square cavity where the leading wall ($x_{1}$), side wall ($x_{2}$), and
lagging wall ($x_{3}$) are equal in length. The formation and growth of
recirculating flow in the vortex or wake bubble (gray area) is characterized
by the transition of the reattachment point ($x_{r}$), normalized by the
length of all stationary cavity walls ($x^{*}=\frac{x_{r}}{\sum x_{i}}$). (b)
Experimental ($\square$) and numerical simulation ($\diamond$) of the
reattachment point ($x^{*}$) logistic growth and transition on the cavity
walls as a function of the channel’s Reynolds number ($Re$). (c) Cavity flow
development at a range of cavity side lengths collapse in a logistic function
presented in a solid red line (-) in (d) when $Re$ is rescaled by the cavity
geometry ratio ($\frac{H}{X_{s}}$). (e) Phase diagram of cavity flow
development states showing regions where phase I: no circulating cavity flow,
phase II: partial cavity flow, and phase III: full cavity flow, as a function
of the cavity aspect ratio ($\frac{H}{X_{s}}$) versus the channel Re. The
dashed red line is a power-law fit of $Re_{m}\sim(\frac{H}{X_{s}})^{-2}$.

There are similarities between the initial stages of vortical flow formation
in a confined 3D cavity and the flow through a sudden expansion or backward-
facing step at low channel Reynolds numbers Biswas  _et al._ (2004). Flow
separates at the sudden expansion region forming an internal recirculating
vortical flow at the leading wall. However, without the existence of a
trailing wall to confine the recirculating flow, secondary and tertiary
separation bubbles occur downstream in the later stages of such flows Durst
_et al._ (1993); Armaly  _et al._ (1983). In the presence of the trailing
wall, flow separation forms a recirculating vortex that grows and occupies the
entire cavity.

Although 2D cavity flow has been extensively studied, many aspects of
recirculating flow formation and mass transport remains unknown in 3D confined
cavities. For example, a defining feature in incompressible and unconfined 2D
cavity flow is the formation of a material flux boundary, or separatrix,
between the main channel flow and open cavity flow Maull and East (1963);
O’Brien (1972); Ghia  _et al._ (1982); Horner  _et al._ (2002). In contrast,
for confined 3D conditions, recirculating inertial flow is highly complex and
unfolds in a qualitatively different manner in three-dimensional space Koseff
_et al._ (1990), causing unexpected complex flow patterns like the ones
generated in wall-confined microchannels Amini  _et al._ (2013) and
T-junctions Vigolo  _et al._ (2014). Therefore, more investigations should
address the non-trivial dynamics of mass flux in 3D cavities.

In this letter, we further study the unexpected material flux into and out of
the cavity flow inspired by initial 3D experimental measurements of flow
streams exiting from the core of a cavity to the mainstream channel flow
without a mass flux boundary (Fig. 1(d-e)). We found no evidence of a
separatrix forming in 3D cavity flow, matching previous numerical studies
indicating that the recirculating wake could not be thought of as containing a
separatrix Torczynski and OHern (1993); Haddadi and Di Carlo (2017). We
further show that the streamline linkage, i.e. locations of entry and exit, to
the main flow is strongly dependent on the development of the cavity flow.
First, we introduced a dimensionless parameter to generalize flow development
stages for different cavity geometries. Then we investigated the mass flux
through the recirculating cavity flow based on the cavity flow development
parameter. Our investigation uncovered an unexpected evolution in the nature
of the channel mass flux to the cavity flow as it develops which has
implications for the ability to transport particles and cells into cavities.

###  .1 Recirculating cavity flow formation

Experimental observations in a channel-cavity system with fixed channel flow
conditions ($Re=uD\textsubscript{h}/\nu$), where $u$ is the average inlet flow
velocity, $D\textsubscript{h}$, hydraulic diameter of the channel, and $\nu$
is the kinematic viscosity of the fluid, demonstrate the development of
recirculating flow in a square cavity as a function of the cavity height ($H$)
and the cavity side length ($X_{s}$) (Fig. S1 and video S1). Given the
complexity in geometry, a Reynolds number alone is not sufficient to define or
generalize cavity flow conditions nor predict recirculating flow formation and
development in microcavities. One can use a Reynolds number along with a
number of non-dimensional geometric parameters instead. Alternatively, one
common defining characteristic shared by all cavity geometries is a transition
in the location of the reattachment line along the cavity wall at which
streamlines diverge between circulating in the cavity, and flowing into the
cavity but returning to the main channel without circulating (Fig. 2(a))
Cherry  _et al._ (1984). The reattachment point $x_{r}$ as observed from a
top-down perspective in the x-y plane, is used to measure the development of
recirculating flow in cavities (Fig. S8). To generalize the transition of the
reattachment point $x_{r}$ in different cavity geometries, (${\textstyle
x^{*}}{\textstyle=}\frac{x_{r}}{\sum x_{i}}$) is used as a measure of the
reattachment point transition distance along the cavity walls normalized to
the total wall lengths of a square cavity with side length $X_{s}$ (scaling
method in supplementary material).

Figure 3:  Cavity-channel mass flux dynamics as a function of cavity vortex
flow development. (a) An exponential increase of cavity flow mass flux
$Q^{*}=\frac{Cavityflux(Q_{ca})}{Channelflux(Q)}$ is observed when the vortex
occupies cavity volume at ($x^{*}>0.9$) and does not significantly expand
further in size. The red line is the power-fit law of $Q^{*}\sim x^{*30}$. (b)
Entry (red) and exit (black) streamline regions into and out of the
recirculating flow are measured at a fixed location upstream and downstream of
the cavity, respectively. For a full cavity flow, entry and exit regions
switch locations from the side to the middle of the channel wall adjacent to
the cavity, and vice versa. The size of the entry region expands towards
higher velocity regions of the main channel flow. (c) The evolution of a
subset of channel streamlines that enter (and exit) the cavity, color mapped
by the by the initial z-height in the channel (Fig S10).

The evolution of the reattachment point in a microcavity ($X_{s}=250\mu m$,
$H=70\mu m$) varies with the entry channel’s Reynolds number ($Re=0.1-100$)
(Fig. 2(b)). At a low Reynolds number, fluid flow passes with fore-aft
symmetry through the cavity, as expected for Stokes flow, with no
recirculation ($x^{*}$ = 0) (fig. 1(d)). Flow separates at higher Reynolds
number creating a separation bubble with a circulating wake. The vortical flow
first evolves over an order of magnitude in Reynolds number, with the
reattachment point $x^{*}$ remaining on the cavity-leading wall ($x_{1}$) as
the Reynolds number increases. Then the reattachment point rapidly migrates
along the cavity sidewall ($x_{2}$), over a small range of Reynolds numbers.
Lastly, recirculating flow separation in the cavity reaches an asymptotic
phase where the reattachment point stagnates at the end of the cavity trailing
wall ($x_{3}$) ($x^{*}$ $\approx$ 1). The evolution of the reattachment point
as a function of the Reynolds number can be described by a logistic function:

$x^{\ast}=\frac{1}{1+e^{-k(Re-a)}}$ (1)

in which $x^{*}$ is the estimated reattachment point location, $k$ is the
curve steepness, and $a$ is the value of the sigmoid’s midpoint. Both
simulation and experimental results reflect the logistic function behavior
(Fig. 2(b)).

The evolution of the reattachment point follows similar forms for different
cavity side lengths ($X_{s}=400-1000\mu m$) (Fig. 2(c)). The vortex growth
curves for different cavity side lengths collapse into (Eq. 1) logistic
relation (red-line in Fig. 2(d)) when it is rescaled by $Re(\frac{H}{X}_{s})$,
where $k=1.31$ and $a=5.5$ (more data in supplementary material Fig. S2-4).
Using this scaling and the logistic function, one can define a condition at
$Re=a$ when the vortical flow occupies approximately half of the total cavity
volume and the reattachment point reaches the middle of the cavity, i.e.
$Re_{m}=Re_{(x^{*}=0.5)}$. $Re_{m}$ is substituted in terms of the logistic
function (Eq. 1) to derive a universal expression that estimates full cavity
flow at $2Re_{m}$ and no cavity circulating flow at $Re_{m}/2$ for any cavity
side length, height, and channel Reynolds number.

$Re_{m\;}=\;a(\frac{X_{s}}{H})$ (2)

To generalize cavity flow transition behavior at any channel-cavity geometry,
we studied cavity flow development as a function of the channel $Re$ along
with the cavity non-dimensional geometric parameter $H/X_{s}$. Experimental
observations are summarized in a phase diagram (Fig. 2(e)) showing the main
cavity flow developmental phases on a range of cavity aspect ratios $0.1\leq
H/X_{s}\leq 1$. Phase I: No cavity flow ($x^{*}=0$) with no recirculation in
the cavity. Phase II: partial cavity flow ($0.1\leq x^{*}<0.9$) after the
formation of a recirculating wake. Phase III: full cavity flow ($x^{*}\geq
0.9$) when the recirculating vortex fills the cavity volume. The red dotted
line represents the power-law fit of $Re_{m}$ at $x^{*}=0.5$. The phase
diagram suggests an inverse relation between recirculating flow development
and cavity aspect ratio. Cavities with low aspect ratio reach full
recirculating cavity flow (phase III) at low channel inertial conditions and
vice versa. Thus, channel-cavity geometry and flow rate can be designed in a
way to either prevent or accommodate recirculating flow.

###  .2 Mass transport in cavity flow

Given the stages in the development of a recirculating vortex, we asked
whether the mass transfer between the main channel and the cavity is dependent
on the vortex reattachment point. Three-dimensional finite element method
simulations of the incompressible Navier-Stokes equations were used to
identify streamlines from the main channel that contribute to the
recirculating cavity flow. We plot in Fig. 3(a) the recirculating volumetric
flux ($Q^{*}$) in the channel versus the reattachment point ($x^{*}=0-0.98$)
in the cavity. The recirculating flux ($Q^{*}$ = $uA$), where $u$ is the
average velocity in $A$ the cross-sectional area, is normalized by the whole
channel flux ($Q$) (Fig. S9). Therefore, this represents the percentage of the
main channel flow that enters (and exits, -conserving mass) the cavity.
Interestingly, the fraction of the main channel flow entering (and exiting)
gradually increases with vortex size in the cavity. In particular, an
unexpected exponential growth of the channel mass flux contribution to cavity
flow is observed when the vortex asymptotically approaches its maximum size
after reaching the cavity volume ($x^{*}\approx 1$).

Coincident with the increase in mass transport, we observe an unexpected
switch in the location of entering and exiting streamlines into the
recirculating cavity flow. Leveraging the flow symmetry in the channel cross-
section, we tracked the evolution of the entrance (red) and exit (black)
locations across a quarter of the channel (Fig. 3(b)). Before the cavity
completely fills ($0.1\leq x^{*}<0.9$), fluid enters the vortex from the
regions near the upper and lower walls and exit along the middle of the
channel side wall downstream of the cavity. Surprisingly, a new entry region
grows at the middle of the cavity as the re-attachment point increases
eventually leading to the entry and exit regions switching locations as the
recirculating flow fills the cavity space further ($x^{*}=0.8-0.9$) (Fig. 3(c)
and Fig. S10). At full cavity flow ($x^{*}>0.9$) the topology of the entering
streamlines switches completely such that streamlines enter at the middle of
the channel sidewall. Coincident with the mass flux increasing, the entry
region area grows towards the center of the channel as the re-attachment point
and Reynolds number increases further. This increase in the entry area
spanning high-velocity regions of the channel cross-section explains the
exponential increase of the recirculating flux in the cavity.

Confocal imaging experiments of selectively dyed streams also reflect the mass
flux into and out of the recirculating flow in the cavity. In our experimental
set-up, we increment the location in the channel cross-section of
fluorescently labeled streams by controlling flow rates: Q1 (rhodamine: red
dye), Q2 center (fluorescein: green dye), and Q3 (water: colorless), and with
Q1=Q2=Q3 set at the same flow rate (Fig. S5). A significant decrease in the
contribution of Q2 streamlines in the recirculating flow is observed in the
transition between ($x^{*}=0.8-0.9$) over which flow entry streamlines are
modeled to change location from the center to the sidewall. For a full cavity
flow ($x^{*}>0.9$), we visualize in 3D the flow exit from the cavity at the
top and bottom of the channel leaving from the core of the vortex. Here
streamlines are selectively stained as follows: Q2 (fluorescein: green dye),
Q1, and Q3(water: colorless) (Fig.1(e) and Fig. S6 and video S2). In a new
cavity design with a notch or step partially blocking the cavity outlet
region, fluorescently labeled streamlines are observed leaving from the core
while diverting from the notch and returning to the mainstream channel flow.
Thus, the flow exit behavior remained consistent with geometric changes at the
cavity outlet region (Fig. S7).

Our numerical and experimental results in sum indicate that there is direct
mass transport from the channel flow into the recirculating flow in the cavity
with no mass flux boundary or separatrix. Besides showing that streamlines
from the main flow enter the recirculating flow and later leave, we found no
evidence of closed streamlines in the recirculating flow that are not
connected to the main flow. Previous theoretical investigation suggests the
same manner of continuous mass exchange between the channel and the cavity
when a full cavity flow is reached Cherry  _et al._ (1984). Our results shed
new light on the increase in the relative mass flux (as a function of main
channel flow rate) as the cavity flow reaches a fully-filled configuration, as
well as the shifting locations of influx and efflux from the cavity. The
enhanced transport between the cavity and the channel as well as shifting
streamline locations that leave the cavity can contribute to the depletion of
particles captured in microvortices at later stages of cavity flow as observed
in previous studies Khojah  _et al._ (2017). This change in fluid flux
direction and magnitude can modulate the limit cycle trajectory of trapped
particles, also experiencing inertial lift forces, to exit the vortex as the
recirculating flow envelops the entire cavity Haddadi and Di Carlo (2017).

These findings expand our understanding of mass transfer in 3D cavity flow
with finite inertia -a widely observed geometry across nature and engineered
systems. Our findings can inform the design of cavities to avoid dead zones of
low mass transport, or enhance transport from specific entry locations in the
main flow to e.g. isolate particles. The understanding of how mass transport
is not limited across a separatrix also has implications for the evolution of
biochemical processes in microvortices inside ocean rocks, physiological flow
and micro-clot accumulation in aneurysms, and engineering microwells for cell
trapping in miniaturized diagnostic devices.

## Acknowledgements

We Thank Dr. Dan Stoecklein and Dr. Kaytlien Hood for helpful discussion and
input. We gratefully acknowledge the support of the Advanced Light
Microscopy/Spectroscopy Center (ALMS) at the California NanoSystems Institute
(CNSI) in the University of California, Los Angeles.

## References

  * Batchelor (2000) G. K. Batchelor,  _An introduction to fluid dynamics_ (Cambridge university press, 2000).
  * Macagno and Hung (1967) E. O. Macagno and T.-K. Hung, Journal of fluid Mechanics 28, 43 (1967).
  * Baloch  _et al._ (1995) A. Baloch, P. Townsend, and M. Webster, Computers & Fluids 24, 863 (1995).
  * Sinha  _et al._ (1982) S. Sinha, A. Gupta, and M. Oberai, AIAA journal 20, 370 (1982).
  * Marcus  _et al._ (2006) J. S. Marcus, W. F. Anderson, and S. R. Quake, Analytical chemistry 78, 956 (2006).
  * Gierahn  _et al._ (2017) T. M. Gierahn, M. H. Wadsworth II, T. K. Hughes, B. D. Bryson, A. Butler, R. Satija, S. Fortune, J. C. Love, and A. K. Shalek, Nature methods 14, 395 (2017).
  * Hur  _et al._ (2011) S. C. Hur, A. J. Mach, and D. Di Carlo, Biomicrofluidics 5, 022206 (2011).
  * Khojah  _et al._ (2017) R. Khojah, R. Stoutamore, and D. Di Carlo, Lab on a Chip 17, 2542 (2017).
  * Dhar  _et al._ (2018) M. Dhar, J. N. Lam, T. Walser, S. M. Dubinett, M. B. Rettig, and D. Di Carlo, Proceedings of the National Academy of Sciences 115, 9986 (2018).
  * Sun  _et al._ (2018) J. Sun, Y. Li, F. Yan, C. Liu, Y. Sang, F. Tian, Q. Feng, P. Duan, L. Zhang, X. Shi, _et al._ , Nature communications 9, 2599 (2018).
  * Rayz  _et al._ (2008) V. Rayz, L. Boussel, M. Lawton, G. Acevedo-Bolton, L. Ge, W. Young, R. Higashida, and D. Saloner, Annals of biomedical engineering 36, 1793 (2008).
  * Biswas  _et al._ (2004) G. Biswas, M. Breuer, and F. Durst, Journal of Fluids Engineering 126, 362 (2004).
  * Durst  _et al._ (1993) F. Durst, J. Pereira, and C. Tropea, Journal of Fluid Mechanics 248, 567 (1993).
  * Armaly  _et al._ (1983) B. F. Armaly, F. Durst, J. Pereira, and B. Schönung, Journal of fluid Mechanics 127, 473 (1983).
  * Maull and East (1963) D. Maull and L. East, journal of Fluid Mechanics 16, 620 (1963).
  * O’Brien (1972) V. O’Brien, The Physics of Fluids 15, 2089 (1972).
  * Ghia  _et al._ (1982) U. Ghia, K. N. Ghia, and C. Shin, Journal of computational physics 48, 387 (1982).
  * Horner  _et al._ (2002) M. Horner, G. Metcalfe, S. Wiggins, and J. Ottino, Journal of Fluid Mechanics 452, 199 (2002).
  * Koseff  _et al._ (1990) J. R. Koseff, A. K. Prasad, C. Perng, and R. L. Street, Physics of Fluids A: Fluid Dynamics 2, 619 (1990).
  * Amini  _et al._ (2013) H. Amini, E. Sollier, M. Masaeli, Y. Xie, B. Ganapathysubramanian, H. A. Stone, and D. Di Carlo, Nature communications 4, 1826 (2013).
  * Vigolo  _et al._ (2014) D. Vigolo, S. Radl, and H. A. Stone, Proceedings of the National Academy of Sciences 111, 4770 (2014).
  * Torczynski and OHern (1993) J. Torczynski and T. OHern,  _Numerical simulations of flow in a three-dimensional cavity-channel geometry_ , Tech. Rep. (Sandia National Labs., Albuquerque, NM (United States), 1993).
  * Haddadi and Di Carlo (2017) H. Haddadi and D. Di Carlo, Journal of Fluid Mechanics 811, 436 (2017).
  * Cherry  _et al._ (1984) N. Cherry, R. Hillier, and M. Latour, Journal of Fluid Mechanics 144, 13 (1984).
  * Haddadi  _et al._ (2018) H. Haddadi, H. Naghsh-Nilchi, and D. Di Carlo, Biomicrofluidics 12, 014112 (2018).