# Solitons in media with mixed, high-order dispersion and cubic nonlinearity

Y. Long Qiang1,Tristram J. Alexander1,C. Martijn de Sterke1,2 1 Institute of
Photonics and Optical Science (IPOS), School of Physics, The University of
Sydney, NSW 2006, Australia.
2 The University of Sydney Nano Institute (Sydney Nano), The University of
Sydney, NSW 2006, Australia<EMAIL_ADDRESS>
###### Abstract

Although most soliton research has traditionally considered dominant quadratic
dispersion, the recent discovery of pure-quartic solitons has inspired
analysis of soliton solutions with large higher orders of dispersion. Here we
present analytic expressions for families of bright soliton solutions at
arbitrary dispersion orders and practical methods to obtain the associated
dispersion relations. These results provide a framework for considering higher
order dispersion solitons and show the potential for further investigation of
solitons in higher order dispersion systems.

  *

##  1 Introduction

Solitons are a striking phenomenon that are unique to physical systems with
nonlinear effects [1, 2, 3]. In optics, solitons can occur from the balance of
group velocity dispersion and the nonlinear Kerr effect, resulting in pulses
which maintain their shape upon propagation. They have been used in the study
and development of various optical applications including telecommunications
[4, 5, 6], supercontinuum generation [7, 8] and ultrafast lasers [9, 10].

Traditionally, the majority of bright soliton studies and applications has
involved negative quadratic dispersion, corresponding to the lowest order
dispersion permitting soliton formation, whereas higher order terms were
treated as a nuisance which needed to be suppressed or managed [11, 12, 13,
14]. More recently, there has been experimental confirmation of the existence
of solitons that rely on the balance of the Kerr effect and pure, negative
higher-order even dispersion [15, 16, 17, 18]. While initial experiments were
carried out in a waveguide with pure negative $4^{\rm th}$ order dispersion,
the development of a fibre laser with programmable dispersion has allowed
demonstration of solitons with pure negative high-order dispersion of up to
$10^{\rm th}$ order [16, 17, 18].

Theoretical studies carried out in parallel have aimed to find analytic
solutions in the presence of high-order dispersion. While no analytic
solutions are known for pure, negative, even dispersion of order $4$ or higher
and a Kerr nonlinearity, some solutions in the presence of mixed orders of
dispersion have been reported. Karlsson and Höök reported a stable stationary
exact solution in the presence of both $2^{\rm nd}$ and $4^{\rm th}$ order
dispersion and a pure Kerr nonlinearity in the form of a hyperbolic secant
squared [19], subsequently confirmed by Piché et al. [20]. More recently, we
reported a solution in the form of a stationary hyperbolic secant cubed in the
presence of $2^{\rm nd}$, $4^{\rm th}$ and $6^{\rm th}$ orders of dispersion
[21]. In 2020 Kudryashov reported a more general set of such solutions, both
stationary and moving, in the form of hyperbolic secant solutions up to power
$6$, as well as solutions consisting of the sum of two terms, each of which is
an integer power of a hyperbolic secant [22]. Kudryashov’s method requires the
solution of a set of algebraic equations, but these become increasingly
complicated and unwieldy as the powers and the number of terms, and
simultaneously the number of dispersion orders, increase. As a consequence
some of the underlying properties of the solutions were not apparent. While
Kudryashov considered a pure Kerr nonlinearity, in subsequent work this was
generalized to systems with more general nonlinearities [23, 24, 25, 26].

Here we develop a systematic method for finding stationary analytic solutions
in the form of superpositions of an arbitrary number of hyperbolic secant
functions of arbitrary integer power, i.e. superpositions of $\rm
sech^{p}(\alpha\tau)$, where $\tau$ is the time in a co-moving frame, $\alpha$
is a width parameter and $p$ is a positive integer. Since our approach is a
systematic one it is straightforward to generalise the solutions reported by
Kudryashov [22]. We fix the form of the solution, and find the associated
dispersion relation. We show that for solutions consisting of a single term,
the associated dispersion relation is associated with a polynomial with even
power terms and $2p$ roots at integer values of $\alpha$, immediately giving
the associated dispersion relation. Since this method cannot be applied to
solutions that are superpositions of integer powers of hyperbolic secants, we
develop another systematic method in which we cast the problem in the language
of linear algebra, and exploit the linear dependence of the equations we
obtain. In this way we straightforwardly generate the relevant algebraic
equations, which can then be solved by standard methods and find solitons not
reported by Kudryashov [22]. We also show that the poles of the spectra of the
solutions lie on the imaginary axis at integer or half-integer multiples of
$i\alpha$.

The outline of this paper is as follows. In Section 2, we introduce the form
of the nonlinear Schrödinger equation to which we find analytic soliton
solutions. In Section 3, we consider solutions consisting of a single term,
reporting the form of these solutions for arbitrary powers and outline how to
find the associated dispersion relation. In Section 4, we consider solutions
consisting of a superposition of the powers of two hyperbolic secant terms. We
also present the systematic method involving linear algebra which we use to
find the associated dispersion relation. In Sec. 5, we generalise the analytic
solution for a superposition to an arbitrary number of terms, and the linear
algebra method to find the associated dispersion relation. In each section, we
also show the unique pole structure in the spectral domain for each of our
solutions.

##  2 Generalised Nonlinear Schrödinger Equation

We consider a variety of analytic soliton solutions in the presence of high
even orders of dispersion. The evolution of these solutions satisfies the
generalised nonlinear Schrödinger equation

$\displaystyle i\frac{\partial\psi}{\partial
z}+\sum_{r=1}^{n}(-1)^{r}\frac{\beta_{2r}}{(2r)!}\frac{d^{2r}\psi}{d\tau^{2r}}+\gamma|\psi|^{2}\psi=0.$
(1)

where $\psi(z,\tau)$ is the pulse envelope, $z$ the propagation coordinate,
$\tau$ the local time, $\gamma$ the nonlinear parameter, which we take to be
positive, and $\beta_{n}$ is the highest order dispersion coefficient. To find
stationary solutions, we set $\psi(z,\tau)=u(\tau;\mu)e^{i\mu z}$, so that $u$
satisfies

$\displaystyle-\mu
u+\sum_{r=1}^{n}(-1)^{r}\frac{\beta_{2r}}{(2r)!}\frac{d^{2r}u}{d\tau^{2r}}+\gamma
u^{3}=0.$ (2)

We note in particular, that since we aim to observe bright, pulse-like
solutions, we require $\beta_{n}<0$ and $\mu>0$ [21].

We briefly consider the linear wave solutions to Eq. (2) in the limit in which
the wave amplitude is small, allowing us to drop the $\gamma u^{3}$ term. The
resulting equation has solutions given as the linear combination of $n$ terms
of the form $\exp(\lambda\tau)$ [17]. The $\lambda$ satisfy the algebraic
equation

$-\mu+\sum_{r=1}^{m}(-1)^{r}\frac{\beta_{2r}}{(2r)!}\lambda^{2r}=0,$ (3)

to which we will refer to as the characteristic polynomial for each soliton
solution.

##  3 Analytic Solution: Single term sechp

Analytic stationary solutions to problems with nonlinearity and dispersion are
fairly rare. However, solutions of the type

$u=A~{}{\rm sech}^{p}(\alpha\tau)$ (4)

are known. The value $p=1$ corresponds to conventional nonlinear Schrödinger
solitons with negative quadratic dispersion [27, 28]. Karlsson and Höök showed
that $p=2$ corresponds to a solution for mixed $2^{\rm nd}$ and $4^{\rm th}$
order dispersion [13]. Similarly, a solution for $p=3$, in the presence of
$2^{\rm nd}$, $4^{\rm th}$ and $6^{\rm th}$ order dispersion has been found
[22, 21]. An obvious question that arises is whether such analysis can be
extended to higher values of $p$, and whether this analysis can be carried out
in a systematic way. We address both of these questions in Section 3.2. Before
doing so, we note that the way we have cast the problem is to assert a
solution (of the form (4)), and then to find the associated dispersion
relation that solves Eq. (2). This contrasts with the usual approach, where,
given a dispersion relation, the associated solution needs to be found.

###  3.1 Sech3 solution

Here we search for the dispersion coefficients $\beta_{n}$ and $\mu$ for which
Eq. (2) has solutions of the form (4) for arbitrary integer $p$. This problem
can be solved by using the relation

$\displaystyle\frac{d^{2}{\rm sech}^{r}(\tau)}{d\tau^{2}}=r^{2}{\rm
sech}^{r}(\tau)-r(r+1){\rm sech}^{r+2}(\tau).$ (5)

Applied recursively, this expression allows even derivatives of hyperbolic
secants to be written as a sum of hyperbolic secant terms. For the highest
order dispersion term with an $n$ power derivative in Eq. (2), we can obtain
hyperbolic secant terms of power $r$ up to $r+n$, and the power of the terms
which can be generated for all the different dispersion terms are shown in
Table 1. In this table, each row corresponds to a term in Eq (2), whereas each
column groups terms of the same power. The existence of a stationary analytic
solution requires the entries in each column to cancel.

| Power of Hyperbolic Secant terms
---|---
| generated by terms in the general NLS
$\mu$ |  $p$ |  |  |  |
$\beta_{2}$ |  $p$ |  $p+2$ |  |  |
$\beta_{4}$ |  $p$ |  $p+2$ |  $p+4$ |  |
$\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\ddots$ |
$\beta_{n}$ |  $p$ |  $p+2$ |  $p+4$ |  $\cdots$ |  $p+n$
$\gamma$ |  |  |  |  |  $3p$
Table 1: Rows show the powers of hyperbolic secant terms generated by each
term through Eq. (4). Columns group the hyperbolic secant terms of equal power

We first note that the nonlinear term needs to cancel the highest power term
generated by the dispersion, i.e., $p+n=3p$, or $n=2p$ (see last column of
Table 1). This enforces the relationship between $n$ and $p$ in the general
case. Choosing the coefficients such that the entries in this column cancel
gives a relation between $\beta_{2p}$, $\gamma A^{2}$, and $\alpha$. The
remaining dispersion coefficients and $\mu$ can be found by solving $p$
simultaneous linear equations. These equations can be written as a triangular
matrix (see Table 1) so they are straightforward to solve. For $p=3$, the
parameters are given in Table 2.

$\gamma A^{2}$ |  $28|\beta_{6}|\alpha^{6}$
---|---
$\alpha^{2}$ |  $-\frac{30}{83}\frac{\beta_{4}}{|\beta_{6}|}$
$\mu$ |  $\frac{245|\beta_{6}|}{16}\alpha^{6}$
$\beta_{2}$ |  $\frac{1891|\beta_{6}|}{360}\alpha^{4}$
Table 2: Exact expressions for the unknowns of Eq. (4) for $p=3$

This approach and the results obtained here are consistent with results found
earlier [22, 21]. Clearly, the number of simultaneous equations that need to
be solved increases with $p$. Whereas this becomes increasingly tedious as $p$
increases, below we present an alternative approach that makes this
unnecessary.

###  3.2 Sechp solution

To extend our method to arbitrary dispersion order, we consider the last
column of Table 1. As discussed, it represents the cancellation of the
nonlinear term in Eq. (2) and the highest power hyperbolic secant from the
dispersion. This term can be generated by recursively applying Eq. (5). Using
this it is found that

$\displaystyle\gamma
A^{2}={\beta_{2p}\over(2p)!}p(p+1)\ldots(3p-1)\alpha^{2p}={1\over
3}{(3p)!\over p!(2p)!}\beta_{2p}\alpha^{2p}.$ (6)

The amplitude in Eq. (4) is thus expressed in the dispersion order and
$\alpha$. For $p=2$ and $p=3$, this generalised expression matches both the
Karlsson and Höök solution and known $6^{\rm th}$ order analytic solution,
given in row 1 of Table 2, respectively [21].

To analyse the linear components for arbitrary order, we note that
characteristic polynomials of the Karlsson and Höök solution and known $6^{\rm
th}$ order solution, have a distinct root structure. We find the roots to be
positioned at $\pm 2\alpha,\pm 4\alpha$, and $\pm 3\alpha,\pm 5\alpha,\pm
7\alpha$, respectively, on the real axis. The roots are thus positioned at
integer values of $\alpha$. This can be understood as follows: the asymptotic
expansion of Eq. (4) is a superposition of terms of the type $e^{\pm
p\alpha\tau}$, $e^{\pm(p+2)\alpha\tau}$, $e^{\pm(p+4)\alpha\tau}$, etc. Since
$e^{\pm 3p\alpha\tau}$ is the third harmonic of the first term, it can be
generated nonlinearly, and similarly for higher powers. Terms between powers
$p$ and $3p-2$ cannot be generated nonlinearly, and must be associated with
the solutions of the characteristic polynomial.

We can thus construct the characteristic polynomial by recognising that it
must have $2p$ roots at $\pm p\alpha,\pm(p+2)\alpha,\ldots,\pm(3p-2)\alpha$.
It must thus be of the form

$\displaystyle\left(\\!{\lambda^{2}}\\!-(\alpha
p)^{2}\right)\\!\left(\\!{\lambda^{2}}\\!-(\alpha(p+2))^{2}\right)\\!\ldots\\!\left(\\!{\lambda^{2}}\\!-(\alpha(3p-2))^{2}\right)\\!=0.$
(7)

apart from a multiplicative constant.

Working backwards, we replace $\lambda$ by $\partial/\partial t$ and find

$\displaystyle\left(\\!{\partial^{2}\over\partial t^{2}}\\!-(\alpha
p^{2})\right)\\!\left(\\!{\partial^{2}\over\partial
t^{2}}\\!-(\alpha(p+2))^{2}\right)\\!\ldots\\!\left(\\!{\partial^{2}\over\partial
t^{2}}\\!-(\alpha(3p-2))^{2}\right)\\!u=0.$ (8)

This expression is equivalent to the linear part of the stationary nonlinear
Schrödinger equation Eq. (2). Multiplying this by $|\beta_{2p}|/(2p)!$ and
comparing with the linear terms of Eq. (2) allows us to obtain the parameters,
and linear dispersion coefficients for the general solution Eq. (4). This
produces expressions which are much easier to find than using the method
outlined in Section 3.1.

We now use this approach to obtain a general expression for $\alpha$. We take
the two highest order equations generated using the method described above,
which provides equations for $\beta_{2p}$ and $\beta_{2p-2}$ in terms of
$\alpha$. Solving these simultaneously, gives

$\displaystyle\alpha^{2}=\frac{(2p-1)}{\frac{13}{24}(2p)^{2}-2p+\frac{1}{3}}\frac{\beta_{2p-2}}{\beta_{2p}}.$
(9)

This expression, in conjunction with the amplitude from Eq. (6) provides
expressions for the unknowns in Eq. (4) at arbitrary order. Note that Eq. (9)
is consistent with row 2 of Table 2 for $p=3$, as required.

With $A$ and $\alpha$ known, we require expressions for the linear dispersion
relation and the nonlinear phase shift. We can match the prefactors between
the expanded polynomial Eq. (8) and the linear part of the stationary NLS (Eq.
(2)), since they are equivalent expressions, which give us the linear
dispersion terms and nonlinear phase term, as functions of $\alpha$,
$\beta_{2p}$ and $\beta_{2p-2}$. Thus, by choosing values, the full linear
dispersion relation and nonlinear phase shift become fully defined. For $p=3$,
the polynomial is

$\displaystyle\frac{\partial^{6}u}{\partial\tau^{6}}-83\alpha^{2}\frac{\partial^{4}u}{\partial\tau^{4}}+1891\alpha^{4}\frac{\partial^{2}u}{\partial\tau^{2}}-11025\alpha^{6}u=0.$
(10)

Drawing an equivalence with the linear part of the associated stationary NLS,
results in the same parameters obtained as in Table 2.

An increase in dispersion order corresponds to an increase in the number of
constraints (number of columns in Table 1), whereas the degrees of freedom
($A$, $\alpha$, $\mu$) remains equal to $3$. Thus, for every increase in
constraints, an additional dispersion term is required. In fact, apart from
$\beta_{2p}$ and $\beta_{2p-2}$, all dispersion terms are prescribed.
Geometrically, the dimension of the corresponding phase space increases with
dispersion order, whereas the dimension of the solution space for these exact
analytic solutions remains the same. For example, the solutions for $p=3$
correspond to a single curve in a $3$-dimensional space when $\beta_{6}$ and
$\gamma$ are fixed [21].

###  3.3 Properties of the solutions

It is straightforward to show that a pulse with an intensity of ${\rm
sech}^{2p}(\alpha\tau)$ has a full-width at half-maximum (FWHM)

$\displaystyle
w=\frac{2}{\alpha}\ln\left(2^{\frac{1}{2p}}+\sqrt{2^{\frac{1}{p}}-1}\right)\approx\frac{2}{\alpha}\sqrt{\frac{\ln
2}{p}},$ (11)

where the last result holds when $p\gg 1$. In this limit $\alpha^{2}$
decreases proportionally with $p$, which cancels with the narrowing of ${\rm
sech}^{2p}(x)$ as $p$ increases; as a consequence, the FWHM takes a limiting
value. Using Eq. (9) we find, as $p\rightarrow\infty$

$\displaystyle w\approx\sqrt{\frac{13\ln
2}{3}\frac{|\beta_{2p}|}{\beta_{2p-2}}}.$ (12)

Figure 1:  Power of solutions of the form (4), with peak power
$P=1.3222~{}{\rm W}$, $\alpha=0.6012~{}{\rm ps^{-1}}$ and $\gamma=1~{}{\rm
W}^{-1}{\rm mm}^{-1}$ kept constant, for $p=3$ (blue: $\beta_{4}=-1\ {\rm
ps^{4}~{}mm}^{-1}$, $\beta_{6}=-1\ {\rm ps^{6}~{}mm}^{-1}$, ${\rm
FWHM}=1.63~{}{\rm ps}$), $p=4$ (red: $\beta_{6}=-1.1985\ {\rm
ps^{6}~{}mm}^{-1}$, $\beta_{8}=-0.8597\ {\rm ps^{8}~{}mm}^{-1}$, ${\rm
FWHM}=1.5152~{}{\rm ps}$), $p=15$ (yellow: $\beta_{28}=-3.0689\ {\rm
ps^{28}~{}mm}^{-1}$, $\beta_{30}=-0.5378\ {\rm ps^{30}~{}mm}^{-1}$, ${\rm
FWHM}=0.9789~{}{\rm ps}$), $p=50$ (purple: $\beta_{98}=-8.7665\ {\rm
ps^{98}~{}mm}^{-1}$, $\beta_{100}=-0.4516\ {\rm ps^{100}~{}mm}^{-1}$, ${\rm
FWHM}=0.5835~{}{\rm ps}$) on (a) a linear scale, and (b) a log scale. Same
pulses versus frequency on (c) linear scale, and (d) log scale. The solutions
on the linear scale have been shifted for clarity.

Figures 1(a) and (b) show the analytic solutions (Eq. (4)) with the same peak
power $P=1.322~{}{\rm W}$ and same width parameter $\alpha=0.6012~{}{\rm
ps^{-1}}$, for $p=3,4,15,50$, on linear and logarithmic scales respectively.
The figures show that the tails of the pulses steepen as $p$ increases,
corresponding to the increasing strength of the dominant term ($e^{\pm
p\alpha\tau}$) of the asymptotic expansion of Eq. (4). The corresponding plots
in the frequency domain in Figs. 1(c) and (d) show an opposite effect, with
steeper solutions in the temporal domain corresponding to broader solutions in
the spectral domain, as expected.

We find that considering the solutions in the spectral domain reveals
intriguing properties, and start by considering in this section the Fourier
transform for a single term analytic solution. The Fourier transform of a
hyperbolic secant is [29]

$\displaystyle{\rm FT}\left({\rm
sech}(\alpha\tau)\right)=\frac{\pi}{\alpha}{\rm
sech}\left(\frac{\pi\omega}{2\alpha}\right),$ (13)

and is thus also a hyperbolic secant. It has simple poles that are evenly
spaced infinitely along the imaginary axis; for our particular function, poles
are found at $\omega=\pm i\alpha,\pm 3i\alpha,\pm 5i\alpha,...$, as indicated
in Fig. 2(a).

Figure 2:  The positions of the poles for the Fourier Transform for (a) ${\rm
sech}(\alpha\tau)$, (b) ${\rm sech^{3}}(\alpha\tau)$ and (c) ${\rm
sech^{2q+1}}(\alpha\tau)$.

The general analytic form of the Fourier transform for odd powers of the
hyperbolic secant is known [29], and can also be shown by induction to be

$\displaystyle{\rm FT}\left({\rm
sech^{2q+1}}(\alpha\tau)\right)=\frac{4^{q}\pi}{(2q)!\alpha}{\rm
sech}\left(\frac{\pi\omega}{2\alpha}\right)\prod^{q}_{r=1}\frac{1}{4\alpha^{2}}\left(\omega^{2}+\alpha^{2}\left(2r-1\right)^{2}\right).$
(14)

The resultant function is the product of Eq. (13) and a polynomial of order
$q^{2}$. This polynomial has roots at $\omega=\pm
i\alpha,...,\pm(2q-1)i\alpha,\pm(2q+1)i\alpha$. These roots supplant the
simple poles of the hyperbolic secant function present at these points. The
position of the poles for the Fourier transform of a cubic hyperbolic secant
($q=1$), and the general hyperbolic secant with odd powers are plotted in Fig.
2(b) and (c) respectively, showing the decreasing number of poles as they are
supplanted by the increasing roots. Thus, for every additional dispersion
order, the number of poles decreases by two, at frequencies of increasing
magnitude.

An analogous argument can be made for the Fourier transform of even power
hyperbolic secants [29]. Their Fourier transforms have the same properties
except that the hyperbolic secant is replaced by a hyperbolic cosecant
multiplied by an odd power polynomial. This solution retains similar
properties to the odd order case, where the positions of the poles are evenly
spaced along the imaginary axis but at multiples of $2i\alpha$ instead
($\omega=\pm 2i\alpha,\pm 4i\alpha,\pm 6i\alpha,...$). Thus the same
principles apply for even power, with simple poles supplanted by roots at the
same frequency. This unique property of these analytic solutions in the
spectral domain becomes increasingly useful as we explore additional analytic
solutions in Sec. 4 and 5.

##  4 Analytic Solution: Two terms sechp+sechp-2

Here we consider a solution composed of the sum of two hyperbolic secants, an
analytic solution in the form

$\displaystyle u=A_{1}\ {\rm sech}^{p-2}(\alpha\tau)+A_{2}\ {\rm
sech}^{p}(\alpha\tau)$ (15)

for $p\geq 4$. The condition to satisfy the NLS remains $n=2p$, for the
highest power dispersion term and nonlinear term to cancel.The key to finding
such solutions in a systematic way is to use the linear dependence of the
equations corresponding to the terms at powers $p-2$ and $p$. Table 3 is the
analogue of Table 1 for the two term solutions (15). The two columns
corresponding to powers $p-2$ and $p$ both have hyperbolic secant terms
generated by each of the linear terms in Eq. (2), enforcing a linear
dependence between the equations at these two powers. This leads to a
relationship between the two amplitudes in Eq. (15), and ultimately allows us
to separate the problem of finding the solution to Eq. (15) from the problem
of finding the associated dispersion relation. As we show in Sec. 4.2, we can
then find the parameters in Eq. (15) for arbitrary order. We then separately
develop a method, which differs from that in Sec. 3, to find the dispersion
relation.

| Power of Hyperbolic Secant terms
---|---
| generated by terms in the general NLS
$\mu$ |  $p-2$ |  $p$ |  |  |  |  |  |
$\beta_{2}$ |  $p-2$ |  $p$ |  $p+2$ |  |  |  |  |
$\beta_{4}$ |  $p-2$ |  $p$ |  $p+2$ |  $p+4$ |  |  |  |
$\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\ddots$ |  $\ddots$ |  $\ddots$ |
$\beta_{n}$ |  $p-2$ |  $p$ |  $p+2$ |  $\cdots$ |  $p+n-6$ |  $p+n-4$ |  $p+n-2$ |  $p+n$
$\gamma$ |  |  |  |  |  $3p-6$ |  $3p-4$ |  $3p-2$ |  $3p$
Table 3: As Table 1 but for substituting Eq. (15) into Eq. (2). The colours of
the columns correspond to the terms of the same colour in Fig. 3.

The reason the method from Sec. 3.2 for finding the dispersion relation does
not apply here is that the substitution of Eq. (15) into the Kerr nonlinear
term leads to four different hyperbolic secant powers (see the last four
columns of Table 3). This increased number of terms due to the nonlinearity
prevents the tail roots being easily predictable, since they no longer all
follow from linear equations as in Section 3.2. As shown by Kudryashov [22],
it is possible to simultaneously solve the equations corresponding to the
columns of Table 3 [22], however this is increasingly tedious as $p$
increases. We therefore look for an alternative method.

Figure 3:  Equations generated at each hyperbolic secant power for $p=5$ with
colours corresponding Table 3. Terms are grouped as follows: (a) from $A\ {\rm
sech}^{3}(\alpha\tau)$; (b) from $B\ {\rm sech}^{5}(\alpha\tau)$; (c)
nonlinear phase shift $\mu u$ and nonlinear terms $\gamma u^{3}$.

###  4.1 Sech3+Sech5 solution

We now take the analytic solution with $p=5$ as an example, giving

$\displaystyle u=A_{1}\ {\rm sech}^{3}(\alpha\tau)+A_{2}\ {\rm
sech}^{5}(\alpha\tau),$ (16)

and identify a process to obtain the unknowns. Figures 3(a)-(c) show all of
the terms generated by a substitution of Eq. (16) into Eq. (2), with the
colours corresponding to the columns of Table 3 for $p=5$.

These correspond to the lowest power columns in Table 3 and all of the
corresponding terms are shown in rows of the same power in Fig. 3. These
produce proportional equations which we simultaneously solve, giving the ratio
$A_{2}=3A_{1}/4$, consistent with known results [22]. We then substitute this
ratio into the unused equations at powers $7,9,11,13,15$, removing one
unknown, and solve the equations for powers $13$ and $15$ simultaneously,
giving exact expressions for $\alpha$ and the amplitude $A_{1}$(see first two
rows in Table 4). Freely choosing values for $\beta_{10}$, $\beta_{8}$ and
$\gamma$ and substituting these into the remaining unused equations, and
simultaneously solving them allows us to find the remaining dispersion terms
and the nonlinear phase shift. Exact expressions for these unknowns can be
found in the last four rows of Table 4. These results mirror existing results
and methods by Kudryashov [22], but encounter the same practical issues as
expressed in Sec. 3 for higher order dispersion. Similar to the single term
case, the number of simultaneous equations to solve increases with dispersion
order, and becomes increasingly tedious. We therefore exploit the equations at
the two lowest hyperbolic secant powers for the general solution Eq. (15), to
develop a different approach.

$\alpha^{2}$ |  $\frac{90}{1157}\frac{\beta_{8}}{\beta_{10}}$
---|---
$\gamma A_{1}^{2}$ |  $-\frac{16016}{9}\alpha^{10}\beta_{10}$
$\mu$ |  $-\frac{97433}{33}\alpha^{10}\beta_{10}$
$\beta_{2}$ |  $\frac{17938}{21}\alpha^{8}\beta_{10}$
$\beta_{4}$ |  $\frac{9448}{33}\alpha^{6}\beta_{10}$
$\beta_{6}$ |  $\frac{27379}{360}\alpha^{4}\beta_{10}$
Table 4: Exact expressions for the unknowns of Eq. (16)

###  4.2 Sechp-2+Sechp solution

To generalise solutions in the form of Eq. (15) to arbitrary consecutive
powers, we start with the first two columns of Table 3. For Eq. (15) to be a
solution, it is necessary that the corresponding equations, can be described
as linearly dependent vectors with the same basis. Thus the entries in each
row must be proportional, such that the prefactor of each ${\rm sech}^{p-2}$
term is linearly proportional to each corresponding ${\rm sech}^{p}$ term. As
an example, the ratio of the prefactors for $\mu$ ($A_{1}$/$A_{2}$) must be
equal to the ratio of the prefactors for $\beta_{2}$
($A_{1}(p-2)^{2}$/($A_{2}p^{2}-A_{1}(p-2)(p-1)$). Solving this equivalence
gives

$\displaystyle A_{1}=\frac{4}{p-2}A_{2},$ (17)

reducing the number of unknowns. We note that the same ratio is found from any
of the other dispersion terms.

We then consider the last two columns of Table 3, corresponding to terms at
powers $3p-2$ and $3p$. We substitute Eq. (17) into these equations and
simplify using some series analysis and factorisation. Solving these
simultaneously gives

$\displaystyle\alpha^{2}=$
$\displaystyle\frac{n(n-1)}{\frac{n}{24}(13n^{2}-24n+8)-2(n-2)+\frac{6(3n-4)(3n-2)}{n-4}}~{}\frac{\beta_{n-2}}{\beta_{n}},$
(18) $\displaystyle\gamma A_{1}^{2}$
$\displaystyle=-\frac{\beta_{n}\alpha^{n}}{3}\left(\frac{8}{n-4}\right)^{2}\frac{(3n/2)!}{n!(n/2)!},$
(19)

where $n=2p$. These results are consistent with the results for $p=5$ in Table
4. Choosing values for $\beta_{n}$, $\beta_{n-2}$ and $\gamma$, the remaining
dispersion terms and nonlinear phase shift can be found fom the remaining
equations. Just like for $p=5$, solving the remaining equations simultaneously
gives exact expressions for the dispersion terms and nonlinear phase shift,
however we have the same practical issues as existing methods as $p$
increases.

The general expressions for the solution are given in Eq. (18) and Eq. (19),
requires only the amplitude ratio Eq. (17) and two highest power equations to
find, thus the majority of the effort then goes into finding the dispersion
coefficients using the remaining equations. Since our main analytic problem is
to solve a set of coupled linear equations, we aim to describe the problem in
matrix form and use methods from linear algebra to simplify the process.

This begins by redefining Eq. (2) using a matrix formalism that incorporates
Eq. (5). This relationship means that every second derivative term enforces
two outcomes, either a hyperbolic secant term at the same order with prefactor
$r^{2}$ or a hyperbolic secant term two orders higher with prefactor
$-r(r+1)$. For our solution Eq. (15), the lowest power hyperbolic secant is
$p-2$. In matrix form, these derivative terms can be expressed as

$\displaystyle D_{2}=\begin{bmatrix}(p-2)^{2}&0&0&0&\ldots\\\
-(p-2)(p-1)&p^{2}&0&0&\ldots\\\ 0&-p(p+1)&(p+2)^{2}&0&\ldots\\\
0&0&-(p+2)(p+3)&(p+4)^{4}&\ldots\\\ \vdots&\vdots&\vdots&\ddots&\ddots\\\
\end{bmatrix}$ (20)

for a $(p+2)\times(p+2)$ matrix. The elements $(j,k)$ give the possible
prefactors from Eq. (5). Column $j$ corresponds to the power of the hyperbolic
secant term which Eq. (5) acts upon, spanning from $p-2$ to $p+n$ in multiples
of two. Row $k$ corresponds to the power of the hyperbolic secant terms
resulting from Eq. (5), which span the same range of hyperbolic secant powers
as the columns. This matrix thus describes the action of a $2^{\rm nd}$ order
derivative term on hyperbolic secant terms from order $p-2$ to $p+n$, and
since the process is linear, we can denote higher order derivative effects as
$D_{r}=D_{2}^{r/2}$.

Next, we describe our analytic solution Eq. (15) as a vector

$\displaystyle U^{T}=\begin{bmatrix}A_{1}&A_{2}&0&0&\ldots&\end{bmatrix}$ (21)

where $A_{1}$ and $A_{2}$ are the amplitudes of our solution from Eq. (15).
This length of the vector is $p+2$ to match the the matrix $D_{r}$. We then
account for the prefactors of each of the linear terms from Eq. (2), including
the dispersion terms, the nonlinear phase shift, and the frequency term from
Eq.(15), expressing it in vector form as

$\displaystyle
B^{T}=\begin{bmatrix}-\mu&(-1)\frac{\beta_{2}}{2!}\alpha^{2}&(-1)^{2}\frac{\beta_{4}}{4!}\alpha^{4}&\ldots&(-1)^{p}\frac{\beta_{2p}}{(2p)!}\alpha^{2p}\end{bmatrix}$
(22)

a vector of length $p+1$. Finally, the nonlinear term $\gamma u^{3}$ can be
expressed as

$\displaystyle
N^{T}=\gamma\begin{bmatrix}\ldots&0&A_{1}^{3}&3A_{1}^{2}A_{2}&3A_{1}A_{2}^{2}&A_{2}^{3}\end{bmatrix},$
(23)

a vector of length $p+2$, where the final four rows correspond to the rows
where nonlinear terms exist, at power $3p-6,3p-4,3p-2,3p$, and where $A_{2}$
can be eliminated by Eq. (17).

Now we can express all terms in Eq. (2) for the solution using matrix
notation, as

$\displaystyle\begin{bmatrix}U&D_{2}U&D_{4}U&\ldots&D_{n}U\end{bmatrix}B=-N.$
(24)

Here the expression in square brackets represents a $(p+2)\times(p+1)$ matrix,
the columns of which consists of the $p+1$ vectors $U,D_{2}U,\dots,D_{n}U$,
and represent the $p+2$ linear equations we are aiming to solve. The resultant
matrix is not square, but this issue is resolved by recalling that the two
linear equations at power $p-2$ and $p$ are proportional, allowing us to
remove the first row of Eq. (24) without losing information. This is true as
long as the first two entries of the nonlinear vector $N$ are also zero. This
again, is only true for $p\geq 4$. This leaves us with a $(p+1)\times(p+1)$
upper triangular square matrix for $p\geq 4$, which can be easily solved as a
set of simultaneous linear equations.

Once this set of equations is solved, we can choose values for the two highest
order dispersion terms ($\beta_{n}$ and $\beta_{n-2}$) and the nonlinear
parameter ($\gamma$), and we are left with a simple set of equations which
describe all the remaining unknowns. This includes $\alpha$ and $A_{1}$, using
the two highest power equations obtained using this matrix process, which is
consistent with Eq. (18) and Eq. (19). This approach using matrices outlines a
systematic method to obtain the corresponding unknowns of the solution. Using
this method for $p=5$ leads to the same values as in Table 4. In A, we see an
example for $p<4$ where a systematic approach cannot be applied since the
linear dependence of the two lowest power equations does not hold.

Figure 4: Power of analytic solutions (15), with peak power $P=0.2484W$,
$\alpha=0.3944~{}{\rm ps^{-1}}$ and $\gamma=1~{}{\rm W}^{-1}{\rm mm}^{-1}$
kept constant for $p=5$ (blue line: $\beta_{8}=-1\ {\rm ps^{8}~{}mm}^{-1}$,
$\beta_{10}=-0.5\ {\rm ps^{10}~{}mm}^{-1}$, ${\rm FWHM}=2.1948~{}{\rm ps}$),
$p=6$ (red line: $\beta_{10}=-1.3280\ {\rm ps^{10}~{}mm}^{-1}$,
$\beta_{12}=-0.7078\ {\rm ps^{12}~{}mm}^{-1}$, ${\rm FWHM}=1.9167~{}{\rm
ps}$), $p=7$ (yellow line: $\beta_{12}=-1.6730\ {\rm ps^{12}~{}mm}^{-1}$,
$\beta_{14}=-0.8965\ {\rm ps^{14}~{}mm}^{-1}$, ${\rm FWHM}=1.7280~{}{\rm
ps}$), $p=10$ (purple line: $\beta_{18}=-2.7830\ {\rm ps^{18}~{}mm}^{-1}$,
$\beta_{20}=-1.3271\ {\rm ps^{20}~{}mm}^{-1}$, ${\rm FWHM}=1.3917~{}{\rm ps}$)
on (a) a linear scale, and (b) a logarithmic scale. Associated spectra on (c)
a linear scale, and (d) a logarithmic scale. The solutions on the linear scale
have been shifted for clarity. Figure 5:  The positions of the poles for the
Fourier Transform for (a) $A_{1}{\rm sech^{3}}(\alpha\tau)+$ $A_{2}{\rm
sech^{5}}(\alpha\tau)$ and (b) $A_{1}{\rm sech^{p-2}}(\alpha\tau)+A_{2}{\rm
sech^{p}}(\alpha\tau)$.

###  4.3 Properties of the solutions

We now compare the properties of the solutions (15), to previous solutions Eq.
(4). The principles of the tail argument outlined in Sec. 3 remains the same:
the smallest real part of the tail solution dominates. For Eq. (15), the
smallest real root corresponds to the first term and so we expect the dominant
term to be $e^{\pm(p-2)\alpha\tau}$ as $|\tau|\rightarrow\infty$. In contrast,
solutions at the same dispersion order for Eq. (4) have a dominant term at
power $p$. The exponential tails for Eq. (15) again become steeper as the
dispersion order increases, by the same argument as in Sec. 3.2. In Fig. 4, we
show solutions (15)) for fixed $\alpha$ and $p=5,6,7,10$. Figure. 4(b) shows
the increasing steepness as the order increases, with the opposite effect
shown in the spectra (Fig. 4(d)).

We now analyze the Fourier transform of Eq. (15), using Eq. (17). After some
rearrangement we find that

$\displaystyle{\rm FT}\left(A_{1}\ {\rm sech^{p-2}}(\alpha\tau)+A_{2}\ {\rm
sech^{p}}(\alpha\tau)\right)=\frac{A_{1}}{4\alpha^{2}}\frac{(\omega^{2}+\alpha^{2}p^{2})}{(p-1)}{\rm
FT}({\rm sech^{p-2}}(\alpha\tau)),$ (25)

from which we can immediately infer the pole structure. Recall from Sec. 3.2
that the Fourier transform of a hyperbolic secant of order $p-2$ has simple
poles at $\omega=\pm(p-2)i\alpha,\pm pi\alpha,\pm(p+2)i\alpha,\ldots$.
However, the simple poles at $\omega=\pm pi\alpha$ are supplanted by the root
at the same point which we observe in the RHS of Eq. (25). Thus the ratio Eq.
(17) effectively removes one of the simple poles (in fact we find this ratio
can be found by enforcing the condition that the simple pole is supplanted).
This changes the pole structure, as shown in Fig. 5.

##  5 Analytic solution: Multiple terms

With the family of analytic solutions established for a sum of two hyperbolic
secant terms, we consider possible families, for an arbitrary number of terms,
of the form

$\displaystyle u=A_{1}\ {\rm sech}^{p-v}(\alpha\tau)+A_{2}\ {\rm
sech}^{p-v+2}(\alpha\tau)+...+A_{y}\ {\rm sech}^{p}(\alpha\tau)$ (26)

for $v$ even, with $p-v$ defining the lowest order term, and the total number
of terms $y=v/2+1$. Table 5 is similar to Tables 1 and 3, showing the
hyperbolic secant terms generated from Eq. (26). We can use the same matrix
method developed in Sec. 4.2 for the general case after some straightforward
modifications. This starts with an adjustment of the dimensions of the matrix
$D_{2}$ to $(p+y)\times(p+y)$, with element (1,1) changed to $p-v$ rather than
$p-2$, and all of the remaining terms adjusted accordingly. The dimensions of
the vectors $U,B$ and $N$ are also changed to $p+y$, with $U$ and $N$
potentially containing additional nonzero elements. We also redefine the
analytic solution $U$ as

$\displaystyle
U^{T}=\begin{bmatrix}A_{1}&A_{2}&\ldots&A_{y}&0&\ldots&\end{bmatrix}$ (27)

for a vector of length $p+y$ containing $y$ nonzero elements.

| Power of Hyperbolic Secant terms
---|---
| generated by terms in the general NLS
$\mu$ |  $p-v$ |  $\cdots$ |  $p-4$ |  $p-2$ |  $p$ |  |  |  |  |  |  |
$\beta_{2}$ |  $p-v$ |  $\cdots$ |  $p-4$ |  $p-2$ |  $p$ |  $p+2$ |  |  |  |  |  |
$\beta_{4}$ |  $p-v$ |  $\cdots$ |  $p-4$ |  $p-2$ |  $p$ |  $p+2$ |  $p+4$ |  |  |  |  |
$\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\vdots$ |  $\ddots$ |  $\ddots$ |  $\ddots$ |  $\ddots$ |
$\beta_{n}$ |  $p-v$ |  $\cdots$ |  $p-4$ |  $p-2$ |  $p$ |  $p+2$ |  $\cdots$ |  $3(p-v)$ |  $3(p-v)+2$ |  $\cdots$ |  $p+n-2$ |  $p+n$
$\gamma$ |  |  |  |  |  |  |  |  $3(p-v)$ |  $3(p-v)+2$ |  $\cdots$ |  $3p-2$ |  $3p$
Table 5: The rows and columns of the table are grouped in the same way as
Table 1 and 3 for a substitution of Eq. (26) into Eq. (2).

In Sec. 4.2 the two lowest power columns were proportional and thus defined
the ratio between amplitudes $A_{1,2}$. In Table 5, the first $y$ columns are
similarly proportional, and allow us to derive a relationship between all the
amplitude terms. We follow the same approach as in Sec. 4, where the ratio of
prefactors for $\mu$ $(A_{1}/A_{r})$ must be equal to the ratio of prefactors
for $\beta_{2}$
($A_{1}(p-v)^{2}/(A_{r}(p-2(y-r))^{2}-A_{r-1}(p-2(y-r)-2)(p-2(y-r)-1))$). This
gives the general ratio between any two successive amplitudes

$\displaystyle\frac{A_{r}}{A_{r-1}}=\frac{(p-2(y-r)-2)(p-2(y-r)-1)}{(p-2(y-r))^{2}-(p-v)^{2}},$
(28)

for arbitrary amplitude $A_{r}$. This allows us to find and inductively define
each of the amplitudes in the general case, in terms of $A_{1}$.

Then the nonlinear term $N$ requires us to cube an arbitrary number of terms,
giving

$\displaystyle N^{T}=\gamma\begin{bmatrix}\ldots\ 0&A_{1}^{3}&3A_{1}^{2}A_{2}\
\ldots\ 6A_{1}A_{2}A_{3}\ \ldots\ 6A_{y-2}A_{y-1}A_{y}\ \ldots\
3A_{y-1}A_{y}^{2}&A_{y}^{3}\end{bmatrix}$ (29)

where the prefactors of each term is determined by the multinomial
coefficient. This results in $3y-2$ nonzero terms in the vector $N$,
corresponding to hyperbolic secant terms between powers of $3(p-v)$ and $3p$.

With all the matrices and vectors from Sec. 4.2 now known for the general
case, we use Eq. (24) to follow the same process, generating a triangular
matrix and solving a set of linear simultaneous equations. After choosing
values for the two highest order dispersion terms ($\beta_{n}$ and
$\beta_{n-2}$) and the nonlinear parameter ($\gamma$), we can obtain all
remaining unknowns. Using the first $y$ columns, we have obtained a
relationship between all the unknown amplitudes allowing us to write all of
the amplitudes in terms of $A_{1}$, and similar to Sec. 4.2, we then solve the
two highest power linear equations simultaneously, to obtain values for
$\alpha$ and $A_{1}$. Solving the remaining unused equations gives us the full
dispersion relation and the nonlinear phase shift for Eq. (26).

This method can be applied provided the nonlinear terms do not overlap with
the set of linearly dependent columns at powers $p-v,\ldots,p$. This condition
can be enforced by ensuring that the size of the elements $(p+y)$, is larger
than or equal to the sum of equivalent columns $y$ and nonlinear terms $3y-2$,
giving the general condition $p\geq 3y-2$. Using this, we see when the power
of the lowest hyperbolic secant term is $1$ (i.e., $p-v=1$), this inequality
does not hold, consistent with our analysis in Sec. 4.

Figure 6: Solutions in the form Eq. (26) for $n=20$ are plotted with
$\alpha=0.3106~{}{\rm ps^{-1}}$, peak power $P=0.6975~{}{\rm mW}$ and
$\gamma=1~{}{\rm W}^{-1}{\rm mm}^{-1}$ kept constant, for $v=0$ (blue curve:
$\beta_{18}=-1\ {\rm ps^{18}~{}mm}^{-1}$, $\beta_{20}=-1\ {\rm
ps^{20}~{}mm}^{-1}$, ${\rm FWHM}=1.7056~{}{\rm ps}$), $v=2$ (red:
$\beta_{18}=-0.5778\ {\rm ps^{18}~{}mm}^{-1}$, $\beta_{20}=-0.4444\ {\rm
ps^{20}~{}mm}^{-1}$, ${\rm FWHM}=1.7688~{}{\rm ps}$), $v=4$ (yellow:
$\beta_{18}=-0.2490\ {\rm ps^{18}~{}mm}^{-1}$, $\beta_{20}=-0.1624\ {\rm
ps^{20}~{}mm}^{-1}$, ${\rm FWHM}=1.8782~{}{\rm ps}$), $v=6$, (purple:
$\beta_{18}=-0.07277\ {\rm ps^{18}~{}mm}^{-1}$, $\beta_{20}=-0.04281\ {\rm
ps^{20}~{}mm}^{-1}$, ${\rm FWHM}=2.1015~{}{\rm ps}$) on (a) a linear scale,
and (b) a logarithmic scale. Associated spectra on (c) a linear scale, and (d)
a logarithmic. The solutions on the linear scale have been shifted for
clarity. Figure 7: The positions of the poles for the Fourier Transform for
(a) $A_{1}{\rm sech^{p-4}}(\alpha\tau)+$ $A_{2}{\rm sech^{p-2}}(\alpha\tau)$
$+A_{3}{\rm sech^{p}}(\alpha\tau)$ and (b) $A_{1}{\rm
sech^{p-v}}(\alpha\tau)+A_{2}{\rm sech^{p-v+2}}(\alpha\tau)+...$ $+A_{y}{\rm
sech^{p}}(\alpha\tau)$.

For the general case, we can also modify the analytic form of $A$ and
$\alpha$, which we found for the two term case in Eq. (19) and Eq. (18). This
requires a simple modification of these equations by the elements dependent on
the ratio between the amplitudes, and correcting these for the general case
gives us the expression

$\displaystyle\alpha^{2}=$
$\displaystyle\frac{n(n-1)}{\frac{n}{24}(13n^{2}-24n+8)-\frac{n-4}{8}\frac{A_{y-1}}{A_{y}}\left(2(n-2)-\frac{6(3n-4)(3n-2)}{n-4}\right)}~{}\frac{\beta_{n-2}}{\beta_{n}},$
(30) $\displaystyle\gamma A_{1}^{2}$
$\displaystyle=-\frac{\beta_{n}\alpha^{n}}{3}\left(\frac{A_{1}}{A_{y}}\right)^{2}\frac{(3n/2)!}{n!(n/2)!},$
(31)

Here we normalised the ratio-dependent terms from Eq. (19) and Eq. (18) and
replaced these with the appropriate ratios we know for the general case. We
find these analytic results to be consistent with the results we obtain
numerically using the matrix method.

###  5.1 Properties of the solutions

In Fig. 6, we now plot four solutions with the same peak dispersion order
$n=20$, using our generalised method to find the unknowns for the solutions.
The principles of the tail argument outlined in Sec. 3 once again apply, where
the smallest real part of the tail solution dominates. Since we expect the
solution of Eq. (26) to have a dominant term at $e^{\pm(p-v)\alpha\tau}$, we
expect the steepest tail of these solutions to be for $v=0$, with an
increasing value of $v$ corresponding to decreasing steepness, which we
observe clearly in Fig. 6(b). The solution in the spectral domain in Fig. 6(d)
shows an opposite effect, where the steeper solutions in the temporal domain
correspond to shallower solutions in the spectral domain, which is what we
would expect.

We also now analyse the Fourier transform for the general case, where a
substitution of the appropriate amplitude using Eq. (28), allows us to
rearrange the general Fourier transform into the form

$\displaystyle{\rm FT}\left(A_{1}\ {\rm sech^{p-v}}(\alpha\tau)+A_{2}\ {\rm
sech^{p-v+2}}(\alpha\tau)+...+A_{y}\ {\rm sech^{p}}(\alpha\tau)\right)=$
$\displaystyle\frac{A_{1}}{(4\alpha^{2})^{v/2}}\frac{(\omega^{2}+\alpha^{2}(p-v+2)^{2})(\omega^{2}+\alpha^{2}(p-v)^{2})...(\omega^{2}+\alpha^{2}p^{2})}{(p-v+1)(p-v+3)...(p-1)}~{}{\rm
FT}({\rm sech^{p-v}}(\alpha\tau)).$ (32)

We recall from Sec. 3.3, that for the Fourier transform of the basic
hyperbolic secant of power one, there are simple poles at $\omega=\pm
i\alpha,\pm 3i\alpha,\pm 5i\alpha,...$, however for an arbitrary hyperbolic
secant at power $r$, all simple poles below this power are supplanted. We then
observed in Sec. 4.3 that for an additional hyperbolic secant term in the
temporal solution, the Fourier transform had an additional simple pole
supplanted at $\omega=\pm pi\alpha$, producing a gap in the simple pole
structure. For the Fourier transform of the general case, we now expect simple
poles at $\omega=\pm(p-v)i\alpha,\pm(p-v+2)i\alpha,\pm(p-v+4)i\alpha,...$, but
now with additional poles supplanted at $\omega=\pm(p-v+2)i\alpha,...,\pm
pi\alpha$, due to the amplitude ratios (Eq. (28)) from the additional roots
generated by the temporal terms. We see this pole structure shown in Fig. 7,
and note that the change in the analytic solution in the temporal domain again
corresponds to a clear change in the structure of the Fourier transform, with
an increasing number of simple poles supplanted corresponding to an increasing
number of terms in Eq. (26). In fact, knowing the structure of the Fourier
transform provides a method of identifying the number and power of the terms
for the solution in the temporal domain. We discuss this method further in
Sec. 6.

##  6 Discussion and Conclusion

The results we have presented provide an analytical framework for exact
solutions for higher order dispersion systems with anomalous dispersion, and
have the potential to be expanded upon both analytically and experimentally.
Recent studies have shown the capability of fibre lasers with programmable
dispersion demonstrating soliton solutions up to $10^{\rm th}$ order
dispersion [16, 17], so experimental observation of the analytic solutions
found in this paper is realistic. Further study of their properties may help
to identify applications.

These analytic solutions can be used to improve the numerical methods used to
find additional solitons as solutions of high-order differential equations. An
effective method to find solutions of Eq. (2) is the Newton Conjugate-Gradient
method [30, 31]. It ideally requires an exact soliton solution at the same
order of dispersion, corresponding to the same order of the differential
equation, as the desired one. The solutions found in this work can therefore
serve as starting points for the numerical method. This then allows us to move
adiabatically through solution space without changing the order of the
differential equation until the solution at the desired parameters is found.
Observing the rate of convergence also allows us to examine the effectiveness
of this numerical method for higher order dispersion.

We have presented results of Eq. (2), which is a stationary equation, but this
does not mean that these solutions are stable. Determining stability requires
the Eq. (1). We carried out simulations of the full evolution of the solutions
we have presented using the split step Fourier method [32]. Although not
carried systematically, all of the solutions we have tested, which include up
to three terms and dispersion up to $14^{\rm th}$ order, appear to be stable.
However, a full stability analysis is outside the scope of this paper.

A consequence of the pole structure of the multi-term solutions, and the
removal of some of the poles, is that the asymptotic expansion of the
solutions in Eq. (26) can be written as

$B_{1}~{}e^{\pm(p-v)\alpha\tau}+B_{y+1}~{}e^{\pm(p+2)\alpha\tau}+B_{y+2}~{}e^{\pm(p+4)\alpha\tau}\ldots,$
(33)

where the $B_{i}$ are linear combinations of the $A_{i}$ from Section 5, and
the terms $B_{2},\ldots,B_{y}$ vanish. In other words, the number of missing
terms in the expansion equals the number of poles that have been removed. One
may conclude that as $v$ becomes very large, then the solution increasingly
resembles the function $e^{-(p-v)\alpha|\tau|}$; however, since this function
has a discontinuous derivative at the origin it is unlikely that it can ever
be reached. Nonetheless, result (28) for the ratio of amplitudes of
consecutive terms in the expression for the soliton can also be found by
requiring that the $2^{\rm nd}$ derivative of Eq. (26) is proportional to the
equation itself. Thus, even though the positions of the simple poles are easy
to predict, the implications of some of the simple poles being supplanted and
the significance of a reduced number of simple poles as dispersion order
increases, remain open questions.

The analytic solutions we construct consist of hyperbolic secant terms with
powers in successive multiples of two. For solutions including successive
terms with a larger difference than two, such as

$\displaystyle u=...+A_{r-d}\ {\rm sech}^{c}(\alpha\tau)+A_{r}\ {\rm
sech}^{c+2d}(\alpha\tau)+...,$ (34)

for cases where $d>1$, we find that the method applied in this paper is not
applicable. The reason for this is simple. In Sec. 5 we found the ratio of
successive amplitudes; for a solution with a gap between successive powers
larger than two, $A_{r-1}=0$ is enforced, for which the only solution is
$r=1$, which prohibits solutions with more than one hyperbolic secant term.
Thus we cannot easily find the ratio between amplitudes, preventing us from
applying the method we developed.

In the temporal domain, we exploited the root structure for the single term
solutions. However for the solutions with multiple hyperbolic secant terms,
some of the characteristic tail roots remain analytically unpredictable due to
the presence of nonlinear terms. Numerical calculations shows that all roots
remain on the real axis for the two-term solitons, with complex roots
introduced for solutions with more terms, and that the roots linearly scale
with $\alpha$. Being able to analytically predict the tail root structure for
an arbitrary solution would allow us to apply the same method as in Sec. 3 to
find the associated dispersion relation, which would be even more
straightforward than our current method.

In summary, we present a study of exact analytic stationary solutions
consisting of a superposition of hyperbolic secant terms of different integer
powers, and a corresponding method to find the associated dispersion
relations. In Sec. 3, we considered solutions with a single term, and derived
their form for arbitrary dispersion order using a constructed polynomial to
define the dispersion relation. Then in Sec. 4, we find the form of solutions
with two terms for arbitrary dispersion order, and present a convenient method
using to find the associated dispersion relation, and in Sec. 5, generalise
these results to an arbitrary superposition of terms. These studies all
include a consideration of the distinct pole structure for any unique solution
in the spectral domain.

We acknowledge funding from the Australian Research Project (ARC) Discovery
Project (DP180102234) and the Asian Office of Aerospace R&D (AOARD) grant
(FA2386-19-1-4067).

##  Appendix A Analytic Solution: $Sech+Sech^{3}$

Here we consider a solution, also found by Kudryashov [22], where the
systematic approach we present in Sec. 4.2 cannot be applied, and so the
solution and the associated dispersion relation cannot be separated. We
consider the analytic solution of the form

$\displaystyle u=A_{1}\ {\rm sech(\alpha\tau)}+A_{2}\ {\rm
sech^{3}}(\alpha\tau).$ (35)

In Sec. 4 we argue that for $p<4$, the linear algebra methods described do not
apply, because the columns at power 1 and 3 are not proportional, due to a
nonlinear term at power 3 (see Fig. 8(c)). This prevents us from using the
ratio Eq. (17) to find a relationship between $A_{1}$ and $A_{2}$, leaving
both as unknowns. Finding this analytic solution requires us to simultaneously
solve all five equations without the simplifications from Sec. 4. The terms at
each power which the equations consist of are grouped by hyperbolic secant
power in Fig. 8, for Eq. (35) substituted into Eq. (2).

Figure 8:  Terms at each hyperbolic secant power when Eq. (35) is substituted
into Eq. (2). (a) terms generated by dispersion applied to $A_{1}\ {\rm
sech}(\alpha\tau)$, and (b) $A_{2}\ {\rm sech}^{3}(\alpha\tau)$, and (c) due
to $\mu u$ and $\gamma u^{3}$. The black box shows the nonlinear term at power
$3$. $D(r)\equiv\beta_{r}\alpha^{r}$.

We aim to make a series of substitutions between these equations, in order to
isolate the relationship between the factors $A_{1}$ and $A_{2}$. We achieve
this by eliminating all $\beta_{2,4,6}$ and $\alpha$ terms. We begin by
obtaining an expression for $\mu$ from the equation at power 1, and
substituting this into the equation at power 3. Following this, we make a
series of substitutions of the equation at power 9, into the equations at
power 7 and 5. This allows us to write the dispersion parameters, and
$\alpha$, in terms of $\gamma$ and $A_{1}$ and $A_{2}$, giving

$\beta_{4}\alpha^{4}=-\frac{83\gamma}{840}(2A_{1}A_{2}+A_{2}^{2}),\quad\beta_{2}\alpha^{2}=-\frac{\gamma}{2520}\\!\left(\\!1177A_{1}^{2}+1387A_{1}A_{2}+\frac{1891A_{2}^{2}}{4}\right)\\!.$
(36)

Substituting these two equations and our equation at power 9 into our updated
equation at power 3, gives

$\gamma(1343A_{1}^{3}+3736A_{1}^{2}A_{2}+3600A_{1}A_{2}^{2}+1152A_{2}^{3})=0$
(37)

dividing through by the nonzero factor $\gamma$ we are left with a cubic
equation composed of only the parameters $A_{1,2}$. Solving these equations
gives $A_{1}=-\frac{12}{17}A_{2}$, consistent with Kudryashov [22].
Substitutions of this ratio into prior equations allows us to obtain
relationships between all the parameters, giving the parameters in Table 6.

$\alpha^{2}$ |  $\frac{510}{581}\frac{\beta_{4}}{|\beta_{6}|}$
---|---
$\gamma A_{2}^{2}$ |  $28\alpha^{6}|\beta_{6}|$
$\mu$ |  $\frac{804}{1627}\alpha^{6}|\beta_{6}|$
$\beta_{2}$ |  $-\frac{92659}{104040}\alpha^{4}|\beta_{6}|$
Table 6: Exact expressions for the unknowns of Eq. (35)

The method used to obtain the unknowns in Eq. (35) are significantly more
demanding than solutions for $p\geq 4$, and we predict that solutions in the
general case from Sec. 5 for the condition $p<3y-2$ will similarly be
practically demanding to solve.

By the same argument as for the exact solution in Sec. 3.2, we consider the
asymptotic expansion of Eq. (35), and similarly expect a superposition of
terms of the type $e^{\pm\alpha\tau}$, $e^{\pm 3\alpha\tau}$, $e^{\pm
5\alpha\tau}$, etc. We note that in this case however, the $e^{\pm\alpha\tau}$
term is the only term generated linearly, as all other terms can be generated
from the third harmonic of Eq. (35). This results in a highly unique root
structure, where applying Eq. (3) produces two real roots at $\pm p\alpha$ and
a quartet of complex roots, which have a smaller real part than the real
roots. Nonetheless, the tails are dominated by the real roots. The pole
structure found for solutions in Sec. 3.3 is also absent.

## References

## References

  * [1] Scott A, Chu F and McLaughlin D 1973 Proceedings of the IEEE 61 1443–1483
  * [2] Polturak E, deVegvar P G N, Zeise E K and Lee D M 1981 Phys. Rev. Lett. 46 1588–1591
  * [3] Akhmediev N N and Ankiewicz A 2011 Dissipative solitons: from optics to biology and medicine (Springer) ISBN 9783642062391
  * [4] Nakazawa M 1994 IEEE Communications Magazine 32 34–41
  * [5] Mollenauer L F, Neubelt M J, Haner M, Lichtman E, Evangelides S G and Nyman B M 1991 Electronic Lett. 27 2055–2056
  * [6] Haus H A and Wong W S 1996 Rev. Mod. Phys. 68 423–444
  * [7] Dudley J M, Genty G and Coen S 2006 Reviews of Modern Physics 78 1135–1184
  * [8] Husakou A V and Herrmann J 2001 Physical Review Letters 87
  * [9] Mitschke F M and Mollenauer L F 1987 Optics Letters 12 407
  * [10] Turitsyn S K, Bale B G and Fedoruk M P 2012 Physics Reports 521 135–203
  * [11] Elgin J N 1992 Opt. Lett. 17 1409–1410
  * [12] Kodama Y, Romagnoli M, Wabnitz S and Midrio M 1994 Opt. Lett. 19 165–167
  * [13] Höök A and Karlsson M 1993 Opt. Lett. 18 1388–1390
  * [14] Aceves A B, Angelis C D, Nalesso G and Santagiustina M 1994 Opt. Lett. 19 2104–2106
  * [15] Blanco-Redondo A, de Sterke C M, Sipe J E, Krauss T F, Eggleton B J and Husko C 2016 Nature Comm. 7 10427
  * [16] Runge A F J, Hudson D D, Tam K K K, de Sterke C M and Blanco-Redondo A 2020 Nat. Photonics 14 492–497
  * [17] Runge A F J, Qiang Y L, Alexander T J, Rafat M Z, Hudson D D, Blanco-Redondo A and de Sterke C M 2021 Physical Review Research 3
  * [18] de Sterke C M, Runge A F J, Hudson D D and Blanco-Redondo A 2021 APL Photonics 6 091101 (Preprint https://doi.org/10.1063/5.0059525) URL https://doi.org/10.1063/5.0059525
  * [19] Karlsson M and Höök A 1994 Opt. Commun. 104 303–307
  * [20] Piché M, Cormier J F and Zhu X 1996 Opt. Lett. 21 845–847
  * [21] Qiang Y L, Alexander T J and de Sterke C M 2022 Phys. Rev. A 105(2) 023501 URL https://link.aps.org/doi/10.1103/PhysRevA.105.023501
  * [22] Kudryashov N A 2020 Applied Mathematics and Computation 371 124972
  * [23] Hosseini K, Sadri K, Mirzazadeh M, Chu Y, Ahmadian A, Pansera B and Salahshour S 2021 Results in Physics 23 104035
  * [24] Kudryashov N A 2021 Optik 235 166626
  * [25] Kudryashov N A 2021 Mathematics 9 3024
  * [26] Arnous A H, Zhou Q, Biswas A, Guggilla P, Khan S, Yıldırım Y, Alshomrani A S and Alshehri H M 2022 Physics Letters A 422 127797
  * [27] Zakharov V and Shabat A B 1972 Soviet Journal of Experimental and Theoretical Physics 34 62
  * [28] Ablowitz M J and Segur H 1981 Solitons and the inverse scattering transform (Society For Industrial And Applied Mechanics) ISBN 9780898714777
  * [29] Erdélyi A and Bateman H 1954 Tables of integral transforms. vol. 1. (New York)
  * [30] Yang J 2009 J. Comput. Phys. 228 7007–7024
  * [31] Yang J 2010 Nonlinear waves in integrable and nonintegrable systems (Society For Industrial And Applied Mathematics) ISBN 9780898717051
  * [32] Agrawal G P 2013 Nonlinear fiber optics (Academic Press) ISBN 9780123970237


Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

Proceedings of the First Workshop on Automatic Speech Recognition for
Spontaneous and Prepared Speech and Speech emotion recognition in Portuguese
(SER 2022), co-located with PROPOR 2022. March 21st, 2022 (Online).

<EMAIL_ADDRESS>]

<EMAIL_ADDRESS>]

# Pretrained audio neural networks for Speech emotion recognition in
Portuguese

Marcelo Matheus Gauy  Universidade de São Paulo, Rua do Matão 1010, São Paulo,
Brazil  Marcelo Finger

(2022)

###### Abstract

The goal of speech emotion recognition (SER) is to identify the emotional
aspects of speech. The SER challenge for Brazilian Portuguese speech was
proposed with short snippets of Portuguese which are classified as neutral,
non-neutral female and non-neutral male according to paralinguistic elements
(laughing, crying, etc). This dataset contains about $50$ minutes of Brazilian
Portuguese speech. As the dataset leans on the small side, we investigate
whether a combination of transfer learning and data augmentation techniques
can produce positive results. Thus, by combining a data augmentation technique
called SpecAugment, with the use of Pretrained Audio Neural Networks (PANNs)
for transfer learning we are able to obtain interesting results. The PANNs
(CNN6, CNN10 and CNN14) are pretrained on a large dataset called AudioSet
containing more than $5000$ hours of audio. They were finetuned on the SER
dataset and the best performing model (CNN10) on the validation set was
submitted to the challenge, achieving an $F1$ score of $0.73$ up from $0.54$
from the baselines provided by the challenge. Moreover, we also tested the use
of Transformer neural architecture, pretrained on about $600$ hours of
Brazilian Portuguese audio data. Transformers, as well as more complex models
of PANNs (CNN14), fail to generalize to the test set in the SER dataset and do
not beat the baseline. Considering the limitation of the dataset sizes,
currently the best approach for SER is using PANNs (specifically, CNN6 and
CNN10).

###### keywords:

Speech emotion recognition Pretrained audio neural networks Transfer learning
Transformers

##  1 Introduction

Speech emotion recognition (SER) aims at identifying the emotional aspects of
speech independently from the actual semantic content. SER can be used to
identify the emotions of humans, e.g., when using mobile phones, an ability
that may become crucial in improving human-machine interactions in the future
[1]. Several efforts to acquire speech data classified with different
emotional labels have been undertaken [2, 3, 4]. These datasets are typically
small in size, even for languages such as English. In order to tackle these
datasets, the use of transfer learning and data augmentation techniques may be
instrumental.

Transfer learning is the method of training a network on a particular problem
where there is an abundance of data, with the goal of using the acquired
knowledge to obtain better performance on a related problem with limited data
available. Transfer learning has been effectively used in many fields of deep
learning such as computer vision [5] and language modelling [6]. Data
augmentation is the method of increasing the amount of data available by
slightly modifying copies of the data. This can be done, for example, by
masking parts of the input or by adding Gaussian noise to it.

In this paper, we use transfer learning and data augmentation techniques to
study SER in Brazilian Portuguese speech. We participate in the shared task
SER challenge, a challenge for Brazilian Portuguese speech emotion
recognition. This challenge made available a labeled dataset of $625$ audio
files as training set for SER. Moreover, a dataset of $308$ files was made
available as the test set. The training and test datasets consisted of short
snippets of Brazilian Portuguese speech, usually less than $15\,s$ long,
labeled neutral, non-neutral female and non-neutral male (non-neutral for
audios containing laughs, cries, etc).

For transfer learning, we employ Pretrained audio neural network (PANN) [7],
which are convolutional neural networks trained on a large dataset of audios
(AudioSet [8]), consisting of $1.9$ million audio clips distributed across
$527$ sound classes. By using the pretrained models made available by the
developers, and finetuning on the SER dataset for Brazilian Portuguese speech,
we are able to beat the proposed baselines of prosodic features and wav2vec
features. We achieve (via CNN10) F1-score of $0.73$, up from $0.54$ from the
baselines. During finetuning, we employ a data augmentation technique called
SpecAugment [9].

We also tested the use of Transformer neural networks, pretrained on a large
amount of Brazilian Portuguese audio data [10]. However, we find that, with
the current amount of available data for SER, Transformers do not generalize
their training performance to the validation and test sets. This holds even
while using most common techniques to prevent overfitting. The same behaviour
was also observed for more complex PANNs, such as CNN14.

##  2 Related Work

There is a large literature on SER in English [11, 12, 13, 14, 15, 16, 17,
18]. Moreover, there are a lot of small datasets for SER in English, such as,
RAVDESS [2], SAVEE [3] and IEMOCAP [4]. To the best of our knowledge, the SER
dataset for Brazilian Portuguese speech is the only available dataset on the
language. In addition, English datasets are usually classified in a different
set of labels. RAVDESS [2], for example, has the classes of calm, happy,
angry, sad, fearful, surprise and disgust. This contrasts with the classes of
neutral, non-neutral female and non-neutral male present in the SER dataset
for Brazilian Portuguese speech. As such, comparisons of our work with the
state of the art in English language are not really possible. Nevertheless,
the authors of [18], the most recent work, obtain an average recall on RAVDESS
of $84.3$ percent using wav2vec 2.0 [19]. On IEMOCAP, they obtain an average
recall of $67.2$ percent, also using wav2vec 2.0.

Transfer learning is a very common technique in situations where the dataset
available is small in size. It has been effectively employed in computer
vision [5, 20], language modelling [6, 21] and audio tasks [7, 22, 18]. In the
original PANN paper [7], authors propose several convolutional neural networks
pretrained on AudioSet which can be finetuned on other smaller datasets. In
[18] the authors use wav2vec 2.0 pretrained on Librispeech and finetuned on
either RAVDESS or IEMOCAP for speech emotion recognition. Finally, in [22] the
authors provide a comprehensive review on transfer learning methods used for
speech and language processing tasks.

##  3 Methodology

The code for this paper can be found at GitHub. Below we describe the dataset
and architectures used.

###  3.1 SER Dataset

To perform SER on Brazilian Portuguese speech, we use the training dataset
(CORAA SER version $1.0$) provided for the challenge. This dataset was built
from the C-ORAL-BRASIL I corpus [23], with $625$ audio files, typically less
than $15s$-long, containing informal spontaneous Brazilian Portuguese speech.
These audio files are labeled neutral, non-neutral female, non-neutral male.
An audio is labeled non-neutral male if it is a male speaker and it contains
paralinguistic elements in the speech (such as laughing, crying, etc).
Similarly, an audio is labeled non-neutral female if it is a female speaker
and the speech contains such paralinguistic elements.

We split the official training dataset into training ($80\%$), validation
($10\%$) and test sets ($10\%$). The split was done in an arbitrary way to
ensure that the three datasets were balanced (i.e. contained relatively the
same proportion of neutral, non-neutral female and non-neutral male files).
The training dataset consisted of $500$ files, the validation dataset
consisted of $63$ files and the test set of $62$ files. The results we report
are for the validation and test set performance.

As the official test dataset made available did not have labels, we have
labeled it ourselves, out of curiosity and to enable more consistent tests of
the performance of the networks. While the labels may not be perfect, they
provide a close enough picture, so the performance of the models can be
measured as an average over multiple experiments (as we were observing high
variance). As such, we also provide results for the official test set with our
unofficial labels. We stress that we did not use the test set labels for any
form of model or parameter selection.

Lastly, the PANNs we use have been trained on the AudioSet [8] dataset
containing more than $5000$ hours of audio distributed across $527$ classes.

###  3.2 PANN Architectures

Table 1 describes the three architectures we use. They are named CNN6, CNN10
and CNN14 after the $6$-layer, $10$-layer and $14$-layer CNNs they represent.
These are the same CNN network architectures used in [7]. We take their
pretrained models on AudioSet [8] to allow us to obtain better generalization
performances on the SER dataset.

The audios are preprocessed in the following way: the audios are first
resampled to $32kHz$. After that, we apply short-time Fourier transform [24]
(with a window size of $1024$ frames and hop size of $320$ frames) to the
standard time-domain waveforms to obtain spectrograms. Then, Mel filter banks
are applied to spectrograms, followed by a logarithm operation to obtain log
Mel spectrograms. These preprocessing steps are commonly done when using CNNs
for audio [25, 26].

As described in Table 1, the CNN architectures used are composed of
convolutional layers with kernel $5\times 5$ for CNN6, and $3\times 3$ for
CNN10 and CNN14. Each convolutional layer is followed by batch normalization
[27] and ReLU non-linearity [28] is used to allow for better training
convergence. Each such convolutional block is present $4$ times in CNN6 and,
in between, an average pooling $2\times 2$ layer is applied (average pooling
is observed to be better than max pooling [29]). In CNN10 and CNN14, the
convolutional blocks are always used in pairs before an average pooling layer
is applied. CNN10 contains $8$ such convolutional blocks ($4$ pairs) and CNN14
contains $12$ such convolutional blocks ($6$ pairs). All networks have a
penultimate fully connected layer to add extra representation ability, as well
as a final $527$ units fully connected layer where a sigmoid is applied to
obtain the probabilities for each class. In Table 1, the first line describes
the input of the networks, that is, $n$ frames of Log Mel Spectrogram with
$64$ mel bins for each frame. Each subsequent line represents a layer of the
networks. The numbers following the @ sign represent the quantity of $5\times
5$ or $3\times 3$ feature maps used.

Table 1: PANN architectures. We describe the layers of CNN6, CNN10 and CNN14. CNN6 |  CNN10 |  CNN14
---|---|---
Log Mel Spectrogram $n$ frames $\times$ 64 mel bins
$\binom{5\times 5@64}{BN,ReLU}$ |  $\binom{3\times 3@64}{BN,ReLU}\times 2$ |  $\binom{3\times 3@64}{BN,ReLU}\times 2$
Avg Pooling $2\times 2$ |  Avg Pooling $2\times 2$ |  Avg Pooling $2\times 2$
$\binom{5\times 5@128}{BN,ReLU}$ |  $\binom{3\times 3@128}{BN,ReLU}\times 2$ |  $\binom{3\times 3@128}{BN,ReLU}\times 2$
Avg Pooling $2\times 2$ |  Avg Pooling $2\times 2$ |  Avg Pooling $2\times 2$
$\binom{5\times 5@256}{BN,ReLU}$ |  $\binom{3\times 3@256}{BN,ReLU}\times 2$ |  $\binom{3\times 3@256}{BN,ReLU}\times 2$
Avg Pooling $2\times 2$ |  Avg Pooling $2\times 2$ |  Avg Pooling $2\times 2$
$\binom{5\times 5@512}{BN,ReLU}$ |  $\binom{3\times 3@512}{BN,ReLU}\times 2$ |  $\binom{3\times 3@512}{BN,ReLU}\times 2$
Global Avg Pooling |  Global Avg Pooling |  Avg Pooling $2\times 2$
FC $512$, ReLU |  FC $512$, ReLU |  $\binom{3\times 3@1024}{BN,ReLU}\times 2$
FC $527$, Sigmoid |  FC $527$, Sigmoid |  Avg Pooling $2\times 2$
|  $\binom{3\times 3@2048}{BN,ReLU}\times 2$
|  Global Avg Pooling
|  FC $2048$, ReLU
|  FC $527$, Sigmoid

###  3.3 Transformer Encoder Architecture

In addition to experimenting with the PANNs, we also attempt to extract good
performances from Transformers. The Transformer architecture we use is
equivalent to the Transformer Encoder architecture from [30]. That is, we use
a three-layer Transformer with multi-head self attention. Each encoder layer
is composed of two sub-layers. The first is a multi-head self-attention
network and the second is a fully connected feed-forward layer. Each sub-layer
has a residual connection followed by layer normalization [31]. The encoder
layers and sub-layers produce outputs of dimension $d$ (in experiments $d$ is
either $128$ or $512$). The fully connected feed forward network within each
encoder layer has an inner dimension of $4d$. We feed the Transformer Encoders
the MFCC-gram of the audios, with each token fed to the Transformer
corresponding to a frame of the MFCC-gram [32]. We name these Transformers,
the MFCC-gram Transformers [32]. We use sinusoidal positional encoding so the
Transformer has access to the order of the sequence fed [30, 33]. The input
frames are projected linearly to a hidden layer of dimension $d$, as direct
addition of acoustic features to positional encoding may lead to training
failure [33].

Typically, Transformers undergo two training phases: pretraining and
finetuning. In the pretraining phase, we make use of a technique called Time
alteration [33] to pretrain the Transformer in about $600$ hours of Brazilian
Portuguese audio data (in other words, we use pretrained models from [10]).
Time alteration is a technique that masks random spans of frames of the MFCC-
gram similarly to how time masking functions in SpecAugment (described in
subsection 3.4). During pretraining, the model is trained to reconstruct the
masked frames. For Brazilian Portuguese audio data, we use the corpora of
NURC-São Paulo[34], NURC-Recife [35], ALIP [36], SP2010 [37] and Programa
Certas Palavras [38]. In the experiments, we also show the performance of
Transformers which do not undergo pretraining, that is, which we initialize at
random and do finetuning directly. We name those Transformers the Baseline
MFCC-gram Transformers. After pretraining, the Transformers are finetuned on
the SER dataset.

###  3.4 Data augmentation: SpecAugment

The SER training dataset used for the challenge leans on the small side and
contains about $50$ minutes of audio. To mitigate the potential overfitting
effects of a small training dataset, we perform a common audio data
augmentation technique called SpecAugment [9] on the Mel spectrogram (or MFCC-
gram) of the audio files before feeding it to the network’s layers.
SpecAugment consists in masking random spans of consecutive segments of the
spectrogram of the audios. Masking can be done along the time dimension (that
is, on spans of consecutive frames), or along the frequency dimension (that
is, on spans of consecutive frequency channels).

Following [7], time masking is done by selecting a uniform length $\ell$
(chosen between $0$ and $64$) and a uniform frame start $t$ (chosen between
$0$ and $T-\ell$, where $T$ is the total number of frames of the audio) and
proceeding to mask the frames from $t$ to $t+\ell-1$. We mask two such blocks
of consecutive frames. Frequency masking is similar to time masking but done
along the frequency dimension. So, a random uniform length $\ell$ is chosen
(between $0$ and $8$) and a uniform frequency band $f$ is chosen (between $0$
and $F-\ell$ where $F$ is the total number of Mel frequency bins). The
frequency bands from $f$ to $f+\ell-1$ are masked to zero. As with time
masking, we mask two such blocks of consecutive frequency bands.

##  4 Results and Discussion

We will check the performance of the three proposed PANNs (CNN6, CNN10 and
CNN14) on the SER training and test datasets. In order to take advantage of
the large pretraining done on the AudioSet [8] dataset, we will use the
pretrained models of CNN6, CNN10 and CNN14 made available by the authors of
[7]. These can be found in Zenodo. These pretrained models will be finetuned
on the SER training dataset in order to achieve better performance than the
baseline.

Moreover, to showcase the massive level of transfer learning that is happening
via the pretrained models, we will show the performance of the three networks
(CNN6, CNN10 and CNN14) without the use of a pretrained model, that is,
initializing their weights at random and not making use of the AudioSet [8]
pretraining. We call these three models the Baseline CNN6, the Baseline CNN10
and the Baseline CNN14.

Lastly, we show the performance of three Transformers models. We analyze MFCC-
gram Transformers pretrained on about $600$ hours of Brazilian Portuguese
audio data, as well as, Baseline MFCC-gram Transformers (without pretraining)
containing $512$ and $128$ units per Encoder layer.

As mentioned before, the SER training dataset is split into a training
($80\%$), validation ($10\%$) and test sets ($10\%$). In Table 2, we report
the $F1$ score performance of the nine models in the validation and test
datasets as well as in the official dataset (which was labeled by us). The
results in the table are averaged across $25$ experiments, to better control
the generally high $F1$ score variance between different experiments. Each
experiment consisted of training the model for $100$ epochs for CNNs and $20$
epochs for Transformers111As Transformer does not generalize, no advantage
exists in training it for longer than $20$ epochs. in the training set and the
best validation performance model (checked after each epoch) was saved and
later analyzed on the test set and official test set. The batch size used was
$16$ and the learning rate was $10^{-4}$ for the CNNs and we use a warmup
learning rate schedule according to the formula $d^{-0.5}\times
min(stepnumber^{-0.5},stepnumber\times warmupsteps^{-1.5})$ for the
Transformers as is standard [6]. We use $warmupsteps=4000$.

As can be seen on Table 2, the best results in the test set were attained by
the CNN6 ($0.62$ $F1$ score). Moreover, it seems that the test set built by us
was inherently harder than the official test set. In the official test set,
the best result was obtained by CNN10 ($0.74$ $F1$ score), in line with it
achieving also the best results on the validation set.

We observe that CNN14’s performance was significantly worse both on validation
and test. However, representation ability wise it is the most powerful of the
PANN models. It is likely that the SER dataset being so small meant CNN14
suffered from overfitting.

We also experienced overfitting issues when attempting MFCC-gram Transformers
based models. There, using pretraining techniques did not yield better
performance. This is likely because the pretraining data contained primarily
voice, without laughs or cries, so the important markers were not present in
pretrained data. Moreover, no common technique (such as dropout [39], L1 or L2
regularization [40], data augmentation techniques as SpecAugment [9] and Mixup
[41]) to prevent overfitting yielded good results. It seems that the reduced
size of the SER dataset is currently hindering performance in more complex
networks, so a likely way of dramatically improving results would be to
increase the size of the available dataset.

Lastly, note that the three baseline PANN models are far away from beating the
baselines provided by the challenge. There is noticeable transfer learning
benefit in using the pretrained models on AudioSet [8]. This large difference
illustrates again the fact that the SER dataset is so small ($50$ minutes of
audio) and that these networks suffer to generalize on it.

We have sent for evaluation in the challenge, the model which attained best
test performance (a CNN6 which officially reported $0.66$ $F1$-score) and the
model which attained best validation performance (a CNN10 which officially
reported $0.73$ $F1$-score). Moreover, out of curiosity, we show the confusion
matrix of the CNN10 model sent for evaluation in Table 3. Observe that the
model classifies the vast majority of neutral and non-neutral females files
correctly. Most of the errors are done classifying non-neutral male files
(often wrongly classified as neutral).

Table 2: The mean and standard deviation of the $F1$ score is shown in the table below for the nine models (CNN6, CNN10 and CNN14 and their respective baseline version, i.e., their versions without pretraining on AudioSet [8], as well as MFCC-gram Transformers with and without pretraining and a smaller version of MFCC-gram Transformers). The results shown are for the validation set, the test set and the official test set. Labels for the official test set were created by us. Model |  $F1$ score Validation performance |  $F1$ score Test performance |  $F1$ score Official test performance
---|---|---|---
Baseline CNN6 |  $0.45\pm 0.06$ |  $0.36\pm 0.05$ |  $0.33\pm 0.03$
Baseline CNN10 |  $0.58\pm 0.06$ |  $0.41\pm 0.09$ |  $0.42\pm 0.05$
Baseline CNN14 |  $0.38\pm 0.06$ |  $0.33\pm 0.04$ |  $0.32\pm 0.03$
CNN6 |  $0.78\pm 0.05$ |  $\mathbf{0.62\pm 0.06}$ |  $0.69\pm 0.04$
CNN10 |  $\mathbf{0.80\pm 0.06}$ |  $0.57\pm 0.06$ |  $\mathbf{0.74\pm 0.04}$
CNN14 |  $0.61\pm 0.11$ |  $0.54\pm 0.06$ |  $0.52\pm 0.10$
MFCC-gram Transformers $512$ units |  $0.50\pm 0.04$ |  $0.36\pm 0.06$ |  $0.38\pm 0.03$
Baseline MFCC-gram Transformers $512$ units |  $0.57\pm 0.04$ |  $0.43\pm 0.08$ |  $0.43\pm 0.06$
Baseline MFCC-gram Transformers $128$ units |  $0.60\pm 0.05$ |  $0.45\pm 0.07$ |  $0.44\pm 0.04$
Table 3: We plot the confusion matrix for the CNN10 model which was submitted to the challenge and attained an $F1$ score of $0.73$ (on official labels). Note that the model has the most difficulty classifying non-neutral male files correctly. Confusion Matrix |  predicted neutral |  predicted non-neutral male |  predicted non-neutral female
---|---|---|---
Neutral |  $244$ |  $2$ |  $5$
Non-neutral male |  $14$ |  $8$ |  $2$
Non-neutral female |  $6$ |  $1$ |  $26$

##  5 Conclusion

In this paper, we have effectively used transfer learning to beat the proposed
baselines in the shared task SER challenge in Brazilian Portuguese speech. By
using, the PANNs CNN6 and CNN10, we have attained $F1$ score of $0.73$ up from
$0.54$ from the baselines. We have also observed that more complex networks,
such as CNN14 and Transformers, while being in theory more capable of
attaining better performances, suffer from overfitting. As such, we determine
that probably the best way of improving results is by increasing the size of
the training set.

Future work could involve increasing the size of the training set so that
Transformers and CNN14 generalize their training performances to the test set.
In addition, pretraining Transformers with audio data containing specifically
laughs, cries and so on may prove useful. Moreover, other data augmentation
techniques could be used which might provide additional benefit in terms of
preventing overfitting.

###### Acknowledgements.

This work was supported by FAPESP grant number 2020/16543-7 (POSDOC) and
project 06443-5 (SPIRA). This work was carried out at the Center for
Artificial Intelligence (C4AI-USP), with support by the São Paulo Research
Foundation (FAPESP) (grant #2019/07665-4) and by the IBM Corporation. Marcelo
Finger was partly supported by the São Paulo Research Foundation (FAPESP)
(grants #2015/21880-4, #2014/12236-1); and the National Council for Scientific
and Technological Development (CNPq) (grant PQ 303609/2018-4). This work was
financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível
Superior – Brasil (CAPES) – Finance Code 001.

## References

  * Andre et al. [2004] E. Andre, M. Rehm, W. Minker, D. Bühler,  Endowing spoken language dialogue systems with emotional intelligence,  in: Tutorial and Research Workshop on Affective Dialogue Systems, Springer, 2004, pp. 178–187.
  * Livingstone and Russo [2018] S. R. Livingstone, F. A. Russo,  The ryerson audio-visual database of emotional speech and song (ravdess): A dynamic, multimodal set of facial and vocal expressions in north american english,  PloS one 13 (2018) e0196391.
  * Wang [2010] W. Wang, Machine Audition: Principles, Algorithms and Systems: Principles, Algorithms and Systems, IGI Global, 2010.
  * Busso et al. [2008] C. Busso, M. Bulut, C.-C. Lee, A. Kazemzadeh, E. Mower, S. Kim, J. N. Chang, S. Lee, S. S. Narayanan,  Iemocap: Interactive emotional dyadic motion capture database,  Language resources and evaluation 42 (2008) 335–359.
  * Voulodimos et al. [2018] A. Voulodimos, N. Doulamis, A. Doulamis, E. Protopapadakis,  Deep learning for computer vision: A brief review,  Computational intelligence and neuroscience 2018 (2018).
  * Devlin et al. [2018] J. Devlin, M.-W. Chang, K. Lee, K. Toutanova,  Bert: Pre-training of deep bidirectional transformers for language understanding,  arXiv preprint arXiv:1810.04805 (2018).
  * Kong et al. [2020] Q. Kong, Y. Cao, T. Iqbal, Y. Wang, W. Wang, M. D. Plumbley,  Panns: Large-scale pretrained audio neural networks for audio pattern recognition,  IEEE/ACM Transactions on Audio, Speech, and Language Processing 28 (2020) 2880–2894.
  * Gemmeke et al. [2017] J. F. Gemmeke, D. P. Ellis, D. Freedman, A. Jansen, W. Lawrence, R. C. Moore, M. Plakal, M. Ritter,  Audio set: An ontology and human-labeled dataset for audio events,  in: 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2017, pp. 776–780.
  * Park et al. [2019] D. S. Park, W. Chan, Y. Zhang, C.-C. Chiu, B. Zoph, E. D. Cubuk, Q. V. Le,  Specaugment: A simple data augmentation method for automatic speech recognition,  arXiv preprint arXiv:1904.08779 (2019).
  * Gauy and Finger [2022] M. Gauy, M. Finger,  Acoustic models for brazilian portuguese speech based on neural transformers,  IN PREPARATION (2022).
  * Lech et al. [2020] M. Lech, M. Stolar, C. Best, R. Bolia,  Real-time speech emotion recognition using a pre-trained image classification network: Effects of bandwidth reduction and companding,  Frontiers in Computer Science 2 (2020) 14.
  * Yoon et al. [2019] S. Yoon, S. Byun, S. Dey, K. Jung,  Speech emotion recognition using multi-hop attention mechanism,  in: ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2019, pp. 2822–2826.
  * Xu et al. [2019] H. Xu, H. Zhang, K. Han, Y. Wang, Y. Peng, X. Li,  Learning alignment for multimodal emotion recognition from speech,  arXiv preprint arXiv:1909.05645 (2019).
  * Yoon et al. [2018] S. Yoon, S. Byun, K. Jung,  Multimodal speech emotion recognition using audio and text,  in: 2018 IEEE Spoken Language Technology Workshop (SLT), IEEE, 2018, pp. 112–118.
  * Satt et al. [2017] A. Satt, S. Rozenberg, R. Hoory,  Efficient emotion recognition from speech using deep learning on spectrograms.,  in: Interspeech, 2017, pp. 1089–1093.
  * Issa et al. [2020] D. Issa, M. F. Demirci, A. Yazici,  Speech emotion recognition with deep convolutional neural networks,  Biomedical Signal Processing and Control 59 (2020) 101894.
  * Peng et al. [2021] Z. Peng, Y. Lu, S. Pan, Y. Liu,  Efficient speech emotion recognition using multi-scale cnn and attention,  in: ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2021, pp. 3020–3024.
  * Pepino et al. [2021] L. Pepino, P. Riera, L. Ferrer,  Emotion recognition from speech using wav2vec 2.0 embeddings,  arXiv preprint arXiv:2104.03502 (2021).
  * Baevski et al. [2020] A. Baevski, Y. Zhou, A. Mohamed, M. Auli,  wav2vec 2.0: A framework for self-supervised learning of speech representations,  Advances in Neural Information Processing Systems 33 (2020) 12449–12460.
  * Krizhevsky et al. [2012] A. Krizhevsky, I. Sutskever, G. E. Hinton,  Imagenet classification with deep convolutional neural networks,  Advances in neural information processing systems 25 (2012) 1097–1105.
  * Brown et al. [2020] T. B. Brown, B. Mann, N. Ryder, M. Subbiah, J. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, et al.,  Language models are few-shot learners,  arXiv preprint arXiv:2005.14165 (2020).
  * Wang and Zheng [2015] D. Wang, T. F. Zheng,  Transfer learning for speech and language processing,  in: 2015 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA), IEEE, 2015, pp. 1225–1237.
  * Raso and Mello [2012] T. Raso, H. Mello,  The c-oral-brasil i: reference corpus for informal spoken brazilian portuguese,  in: International Conference on Computational Processing of the Portuguese Language, Springer, 2012, pp. 362–367.
  * Brigham and Morrow [1967] E. O. Brigham, R. Morrow,  The fast fourier transform,  IEEE spectrum 4 (1967) 63–70.
  * Choi et al. [2016] K. Choi, G. Fazekas, M. Sandler,  Automatic tagging using deep convolutional neural networks,  arXiv preprint arXiv:1606.00298 (2016).
  * Kong et al. [2019] Q. Kong, C. Yu, Y. Xu, T. Iqbal, W. Wang, M. D. Plumbley,  Weakly labelled audioset tagging with attention neural networks,  IEEE/ACM Transactions on Audio, Speech, and Language Processing 27 (2019) 1791–1802.
  * Ioffe and Szegedy [2015] S. Ioffe, C. Szegedy,  Batch normalization: Accelerating deep network training by reducing internal covariate shift,  in: International conference on machine learning, PMLR, 2015, pp. 448–456.
  * Nair and Hinton [2010] V. Nair, G. E. Hinton,  Rectified linear units improve restricted boltzmann machines,  in: Icml, 2010.
  * Kong et al. [2019] Q. Kong, Y. Cao, T. Iqbal, Y. Xu, W. Wang, M. D. Plumbley,  Cross-task learning for audio tagging, sound event detection and spatial localization: Dcase 2019 baseline systems,  arXiv preprint arXiv:1904.03476 (2019).
  * Vaswani et al. [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, I. Polosukhin,  Attention is all you need,  Advances in neural information processing systems 30 (2017) 5998–6008.
  * Ba et al. [2016] J. L. Ba, J. R. Kiros, G. E. Hinton,  Layer normalization,  arXiv preprint arXiv:1607.06450 (2016).
  * Gauy and Finger [2021] M. M. Gauy, M. Finger,  Audio mfcc-gram transformers for respiratory insufficiency detection in covid-19,  in: STIL 2021 (), 2021. URL: http://XXXXX/219270.pdf.
  * Liu et al. [2020] A. T. Liu, S.-w. Yang, P.-H. Chi, P.-c. Hsu, H.-y. Lee,  Mockingjay: Unsupervised speech representation learning with deep bidirectional transformer encoders,  in: ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2020, pp. 6419–6423.
  * Castilho and Pretti [1986] A. T. d. Castilho, D. Pretti, A linguagem falada culta na cidade de são paulo: materiais para seu estudo., 1986.
  * Oliviera Jr et al. [2016] M. Oliviera Jr, et al.,  Nurc digital um protocolo para a digitalização, anotação, arquivamento e disseminação do material do projeto da norma urbana linguística culta (nurc),  CHIMERA: Revista de Corpus de Lenguas Romances y Estudios Lingüísticos 3 (2016) 149–174.
  * Gonçalves [2019] S. C. L. Gonçalves,  Projeto alip (amostra linguística do interior paulista) e banco de dados iboruna: 10 anos de contribuição com a descrição do português brasileiro,  Estudos Linguísticos (São Paulo. 1978) 48 (2019) 276–297.
  * Mendes [2013] R. B. Mendes,  Projeto sp2010: Amostra da fala paulistana,  http://projetosp2010. fflch. usp. br>. Acesso em 1 (2013) 2013.
  * Teixeira [1997] C. S. P. Teixeira, Acervo Certas Palavras- Catálogo 1981-1996., Unicamp Cedae, 1997\.
  * Srivastava et al. [2014] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, R. Salakhutdinov,  Dropout: a simple way to prevent neural networks from overfitting,  The journal of machine learning research 15 (2014) 1929–1958.
  * Goodfellow et al. [2016] I. Goodfellow, Y. Bengio, A. Courville, Deep learning, MIT press, 2016.
  * Zhang et al. [2017] H. Zhang, M. Cisse, Y. N. Dauphin, D. Lopez-Paz,  mixup: Beyond empirical risk minimization,  arXiv preprint arXiv:1710.09412 (2017).