# Further Compactifying Linear Optical Unitaries

B.A. Bell<EMAIL_ADDRESS>I.A. Walmsley  Ultrafast Quantum Optics
group, Department of Physics, Imperial College London, London SW7 2AZ, UK

###### Abstract

Quantum integrated photonics requires large-scale linear optical circuitry,
and for many applications it is desirable to have a universally programmable
circuit, able to implement an arbitrary unitary transformation on a number of
modes. This has been achieved using the Reck scheme, consisting of a network
of Mach Zehnder interferometers containing a variable phase shifter in one
path, as well as an external phase shifter after each Mach Zehnder. It
subsequently became apparent that with symmetric Mach Zehnders containing a
phase shift in both paths, the external phase shifts are redundant, resulting
in a more compact circuit. The rectangular Clements scheme improves on the
Reck scheme in terms of circuit depth, but it has been thought that an
external phase-shifter was necessary after each Mach Zehnder. Here, we show
that the Clements scheme can be realised using symmetric Mach Zehnders,
requiring only a small number of external phase-shifters that do not
contribute to the depth of the circuit. This will result in a significant
saving in the length of these devices, allowing more complex circuits to fit
onto a photonic chip, and reducing the propagation losses associated with
these circuits. We also discuss how similar savings can be made to alternative
schemes which have robustness to imbalanced beam-splitters.

††preprint: APS/123-QED

##  I Introduction

Optical quantum computing requires interferometric circuits to process quantum
states of light O’Brien (2007); Kok  _et al._ (2007). In particular,
integrated photonic circuits comprising networks of Mach Zehnder
interferometers (MZIs) have emerged as a compact and versatile solution for
realising reconfigurable linear optics, with applications to linear optical
quantum computing Carolan  _et al._ (2015); Metcalf  _et al._ (2014); Qiang
_et al._ (2018), boson sampling Spring  _et al._ (2013); Carolan  _et al._
(2014); Crespi  _et al._ (2013); Bell  _et al._ (2019), high-dimensional
encodings Wang  _et al._ (2018); Taballione  _et al._ (2019), quantum
simulation Sparrow  _et al._ (2018); Harris  _et al._ (2017), photonic neural
networks Shen  _et al._ (2017), and optical FPGAs Zhuang  _et al._ (2015);
Pérez  _et al._ (2017). Often, universal reconfigurability is essential or at
least desirable, in the sense that a device can be programmed to realise any
unitary transformation between the input and output modes. This can be
achieved using the architecture of Reck et al. (the Reck scheme) Reck  _et
al._ (1994), shown in Fig. 1(a), a triangular mesh with each unit cell
comprising a variable beam-splitter and a variable phase-shift. The
rectangular architecture of Clements et al. Clements  _et al._ (2016), Fig.
1(b), is also universal, and benefits from improved compactness as well as the
balanced loss per channel, since each path has the same number of unit cells,
that improves the overall fidelity of the circuit to the desired operation. In
integrated photonics, the unit cell is usually implemented as an MZI with an
internal phase-shifter in one of the arms, and an external phase-shifter on
one output, as in Fig. 1(c), which we will term an asymmetric MZI (aMZI). The
internal phase-shifter controls the splitting ratio between the two outputs,
while the external controls the relative phase between the outputs. The reason
for using this primitive operation is that it is a universal 2x2 circuit.

Alternatively, Fig. 1(d) shows a symmetric MZI (sMZI) with two internal and no
external phase-shifters. This lacks control of the relative phase between
outputs, yet it has been used in optical FPGAs Zhuang  _et al._ (2015); Pérez
_et al._ (2017) and it has been found that the sMZI can replace the aMZIs in a
Reck scheme without compromising universality Miller (2015); Ribeiro  _et al._
(2016). The sMZI is attractive because it is more compact, without the need
for an external phase-shifter, which can account for a significant fraction of
the length of the circuit. A shorter structure not only occupies less area on
a chip, but also suffers from less propagation loss. Further, there is a
potential advantage in circuit loading from the controller drive signals,
which benefits from the symmetric configuration in terms of e.g. heat
distribution across the circuit. A rectangular architecture such as the
Clements scheme but made up of sMZIs would be particularly beneficial for
quantum integrated photonics, where large-scale circuits are required and
transmission must be kept as high as possible. However sMZIs have not as yet
been utilised in this way, because it is not clear that such a non-universal
2x2 element can generate a universal circuit, or that there is an efficient
algorithm for finding the phases which will implement a target unitary.

Figure 1:  (a) A Reck scheme, where each blue square represents a tunable
beamsplitter and a phase-shifter. (b) The Clements scheme. (c) The tunable
beam-splitters can be realsied by MZIs with one internal phase-shifter. One
external phase-shifter is needed to complete the unit cell. We term this an
asymmetric MZI. (d) The symmetric MZI, with two internal phase shifts, makes
for a more compact unit cell.

In this work, we resolve this issue. We show that sMZIs can replace aMZIs in
the Clements schemes, providing methods of matrix decomposition which generate
the specific set of phases required to implement a given unitary matrix. We
initially find that additional external phase-shifters are required in a
diagonal layer through the middle of the circuit, as well as at a subset of
the inputs and outputs. We subsequently show that the additional mid-circuit
phase-shifts can be moved to otherwise vacant positions at the edge of the
circuit, where they do not contribute to the overall length. This
approximately halves the contribution of the phase-shifters to the length of a
circuit, down to at most external phase-shifters at the inputs and outputs of
the circuit, which are not required for some applications. We show how a
similar advantage can be obtained for more general circuits comprising
alternating layers of beam-splitters and phase-shifters, giving as an example
the design of Fldzhyan et al. Fldzhyan  _et al._ (2020), which does not have a
deterministic method of setting the phases but does appear to heuristically
implement arbitrary unitary matrices with an improved level of robustness to
imbalanced beam-splitters.

The remainder of this article is set out as follows. In Section II, we provide
a matrix decomposition method for a Reck scheme made up of sMZIs. In Section
III, we follow a near-identical method for the Clements scheme, showing how a
layer of external phase-shifters arises mid-circuit. In Section IV, we propose
a variant on the Clements scheme with additional phase-shifters at the edge of
the circuit, which do not add to the overall length, and show that these can
fulfil the same function as the external mid-circuit phase-shifts. In Section
V we discuss how a similar technique can be applied to more general
architectures, such as that of Fldzhyan et al. Fldzhyan  _et al._ (2020). In
Section VI we summarise and draw conclusions.

##  II The Reck Scheme

We consider a Reck scheme Reck  _et al._ (1994) where the unit cell is a sMZI,
with additional external phase-shifts at the input and output of the entire
circuit, on each mode except the first, as depicted in Fig. 2(a). The matrix
transformation implemented by one sMZI, applied to two adjacent modes, can be
written as

$M=\begin{pmatrix}1~{}~{}&&&&&\\\ &\ddots&&&&\\\
&&e^{i\Sigma}~{}\mathrm{sin}~{}\delta&e^{i\Sigma}~{}\mathrm{cos}~{}\delta&&\\\
&&e^{i\Sigma}~{}\mathrm{cos}~{}\delta&-e^{i\Sigma}~{}\mathrm{sin}~{}\delta&&\\\
&&&&\ddots&\\\ &&&&&~{}~{}1\end{pmatrix},$ (1)

where $\Sigma=(\theta_{1}+\theta_{2})/2$, $\delta=(\theta_{1}-\theta_{2})/2$,
and $\theta_{1,2}$ are the values of the two internal phase-shifts. Meanwhile,
an external phase-shift is given by

$P=\begin{pmatrix}1~{}~{}&&&&\\\ &\ddots&&&\\\ &&e^{i\phi}&&\\\ &&&\ddots&\\\
&&&&~{}~{}1\end{pmatrix}.$ (2)

The circuit is divided into diagonals, numbered $j=1$ to $m-1$, where $m$ is
the number of modes in the circuit. The MZIs within each diagonal are numbered
from bottom left to top right as $k=1$ to $j$. The MZI transformations within
the circuit can be identified as $M^{(j,k)}$ with corresponding phase settings
$\Sigma_{j,k}$ and $\delta_{j,k}$, while the input phase-shift operations
associated with a diagonal are denoted $P^{(j)}$ with phase $\phi_{j}$. The
output phase-shift operation applied to mode $j$ is labelled $Q^{(j)}$ with
phase $\zeta_{j}$.

Figure 2:  (a) Layout of a Rech scheme with symmetric MZIs (blue squares) and
external phase-shifters (yellow triangles). The elements are divided into
diagonal rows labelled by $j$, within a row MZIs are labelled by $k$. (b)
Steps (1) to (3) of the decomposition process. Circuit elements circled in red
affect corresponding matrix elements circled in red, and their parameters are
chosen to zero a succession of matrix elements.

The matrix decomposition proceeds by applying the circuit operations to an
auxillary matrix $V$ (initially set to $U^{*}$) so as to successively zero
matrix elements. Fig. 2(b) shows the first few operations used to zero
elements. Multiplying $V$ from the right by an $M^{(j,k)}$ matrix mixes two
columns of $V$ together - specifically the $y=j-k+1$ column is mixed with the
$y+1$ column. A particular element $(x,y)$ can be set to zero by choosing
$\delta_{j,k}$ such that:

$\mathrm{tan}~{}\delta_{j,k}=-V_{x,y}/V_{x,y+1}.$ (3)

This has a real solution for $\delta_{j,k}$ if $V_{x,y}$ and $V_{x,y+1}$ have
the same complex phase, i.e. $\mathrm{arg}(V_{x,y})=\mathrm{arg}(V_{x,y+1})$.
Since $\Sigma_{j,k}$ affects the phase of both columns equally, it cannot be
chosen to achieve this condition - rather the phases need to be equalized by
previous operations. The external phase-shifter $P^{(j)}$ can be used to match
the phases for the first MZI in each diagonal, then for each MZI
$\Sigma_{j,k}$ can be chosen to match the phases for the subsequent MZI. The
detailed order of operations is given as follows:

  1. 1.

Set an auxillary matrix $V\leftarrow U^{*}$.

  2. 2.

For $j=1$ to $m-1$:

     * •

Set $x\leftarrow m$ and $y\leftarrow j$.

     * •

Set $V\leftarrow VP^{(j)}$, choosing
$\phi_{j}=\mathrm{arg}(V_{x,y})-\mathrm{arg}(V_{x,y+1})$.

     * •

For $k=1$ to $j$:

       * –

Set $V\leftarrow VM^{(j,k)}$, choosing $\delta_{j,k}$ such that $V_{x,y}$ is
set to zero.

       * –

Choose $\Sigma_{j,k}$ such that
$\mathrm{arg}(V_{x-1,y-1})=\mathrm{arg}(V_{x-1,y})$. For $k=j$, this choice is
redundant.

       * –

Set $x\leftarrow x-1$ and $y\leftarrow y-1$.

  3. 3.

For $j=2$ to $m$:

     * •

Set $V\leftarrow VQ^{(j)}$, choosing
$\zeta_{j}=\mathrm{arg}(V_{1,1})-\mathrm{arg}(V_{j,j})$

After Step 2, every element of $V$ below the diagonal has been set to zero.
Since $V$ remains unitary, this implies that it is a diagonal matrix, where
the remaining diagonal elements are complex with unit norm. Step 3 then
applies the final phase-shifts such that $V$ is the identity matrix
$\mathds{1}$, up to a global phase. Expanding $V$ as

$U^{*}P^{(1)}M^{(1,1)}...~{}M^{(m-1,m-1)}Q^{(2)}...~{}Q^{(m)}=\mathds{1},$ (4)

it can be seen that

$U=Q^{(m)}...~{}Q^{(1)}M^{(m-1,m-1)}...~{}M^{(1,1)}P^{(1)},$ (5)

where we have used the fact that all of the circuit operations are symmetric
unitaries, and so their complex conjugate is their inverse. Hence $U$ has been
decomposed into the individual circuit operations.

The circuit consists of $\frac{1}{2}m(m-1)$ MZI, and so $m(m-1)$ internal
phase-shifters. In each row of MZI, the final choice of $\Sigma_{j,j}$ is
redundant, and so in these MZI one of the internal phase-shifters could be
omitted. This would leave $(m-1)^{2}$ internal phase-shifters and $2(m-1)$
external. This is equal to the number of free parameters in an $m\times m$
unitary matrix, and hence is the minimum required. We note that in situations
where the output modes are to be connected directly to phase-insensitive
detectors, the external phase-shifts at the output are redundant and could be
omitted. In some situations, the external phases at the inputs could also be
redundant, for instance if single photon input states are used.

##  III The Clements Scheme

We now consider the Clements scheme Clements  _et al._ (2016) circuit shown in
Fig. 3(a). The MZI are organised into diagonals labelled by $j$ as in the
previous section, but for even $j$ the direction of $k$ has been reversed - as
a result, the associated $P^{(j)}$ phase-shifters have been moved from input
to output. In the decomposition, the order in which matrix elements are zeroed
within an even $j$ diagonal is also reversed, and the corresponding circuit
operations are applied to $V$ by left multiplication rather than right
multiplication. As a result these operations mix adjacent rows of the $V$
matrix rather than adjacent columns - the first few steps of this process are
shown in Fig. 3(b). The $Q^{(j)}$ operations appear in a diagonal line through
the middle of the circuit, which is not a usual feature of the Clements scheme
- in the next section we will show how these phase-shifts can be relocated to
positions where they do not add to the overall length. In other respects the
decomposition is identical to that of the Reck scheme.

Figure 3:  (a) Labelling of diagonal rows of MZI and ordering within a row for
the Clements scheme. For even rows, the ordering is reversed, beginning from
the outputs. (b) Choosing the phase-shifts in the first two rows. (i) applies
the MZI in the 1st row by multiplying the auxillary matrix from the right,
mixing the first two columns together to zero the bottom left element. (ii)
and (iii) implement the MZIs in the 2nd row to zero further elements, by
multiplying the auxillary matrix from the left. This has the effect of mixing
two rows instead of two columns. (c) The residual phases are left mid-circuit
rather than at the outputs as with the Reck scheme. However they can be
absorbed into the action of surrounding phase-shifters with minimal changes to
the circuit.

  1. 1.

Set an auxillary matrix $V\leftarrow U^{*}$.

  2. 2.

For $j=1$ to $m-1$:

     * •

if $j$ is odd:

       * –

Set $x\leftarrow m$ and $y\leftarrow j$.

       * –

Set $V\leftarrow VP^{(j)}$, choosing
$\phi_{j}=\mathrm{arg}(V_{x,y})-\mathrm{arg}(V_{x,y+1})$.

       * –

For $k=1$ to $j$:

         * *

Set $V\leftarrow VM^{(j,k)}$, choosing $\delta_{j,k}$ to zero $V_{x,y}$.

         * *

Choose $\Sigma_{j,k}$ so that
$\mathrm{arg}(V_{x-1,y-1})=\mathrm{arg}(V_{x-1,y})$.

         * *

Set $x\leftarrow x-1$ and $y\leftarrow y-1$.

     * •

if $j$ is even:

       * –

Set $x\leftarrow m-j+1$ and $y\leftarrow 1$.

       * –

Set $V\leftarrow P^{(j)}V$, choosing
$\phi^{(j)}=\mathrm{arg}(V_{x,y})-\mathrm{arg}(V_{x-1,y})$.

       * –

For $k=1$ to $j$:

         * *

Set $V\leftarrow M^{(j,k)}V$, choosing $\delta_{j,k}$ to zero $V_{x,y}$.

         * *

Choose $\Sigma_{j,k}$ so that
$\mathrm{arg}(V_{x+1,y+1})=\mathrm{arg}(V_{x,y+1})$.

         * *

Set $x\leftarrow x+1$ and $y\leftarrow y+1$.

  3. 3.

For $j=2$ to $m$:

     * •

Set $V\leftarrow VQ^{(j)}$, choosing
$\zeta_{j}=\mathrm{arg}(V_{1,1})-\mathrm{arg}(V_{j,j})$.

As before, at the end of Step 3, $V$ is the identity matrix up to a global
phase. However, now when expanding $V$:

$\begin{split}...M^{(4,4)}...M^{(2,1)}P^{(2)}~{}U^{*}~{}P^{(1)}M^{(1,1)}...M^{(3,3)}...&~{}Q^{(2)}...Q^{(m)}\\\
&=\mathds{1},\end{split}$ (6)

i.e. the operations with even $j$ are to the left of $U^{*}$ rather than the
right. It can be seen that this naturally places the $Q^{(j)}$ operations
which were applied last in the middle of the circuit:

$U=P^{(2)}M^{(2,1)}...M^{(4,4)}...Q^{(2)}...Q^{(m)}...M^{(3,3)}...M^{(1,1)}P^{(1)}.$
(7)

##  IV Relocating residual phase-shifts

We now propose a variant on the Clements scheme as shown in Fig. 4(a),
consisting of a rectangular network of sMZI. For each vertical layer of MZIs,
there is one odd mode at the edge which is not involved in a MZI (for an even
total number of modes, alternating layers have zero or two modes not involved
in a MZI). The only change compared to the circuit considered in the previous
section is that tunable phase-shifters are added to these sections of
waveguide, instead of placing the $Q^{(k)}$ phase-shifters in the middle of
the circuit. The new phase-shifters do not add to the overall length of the
circuit, since path-length matching dictates that these sections of waveguide
are anyway as long as if they had been involved in a MZI. The external phase-
shifters at the inputs and outputs are not shown here.

Figure 4:  (a) A Clements scheme circuit consisting of symmetric MZIs, with
additional phase-shifters on the waveguides at the edges of the circuit.
(b)(i) A desired phase-shift $\phi$ between layers of MZIs is implemented by
(ii) adding $+\phi$ to both phases in an adjacents MZI; this moves the
required phase to one of opposite sign on the mode immediately below the
desired target phase. (iii) Repeating this operation until the residual phase
is at the edge of the circuit, it can then be implemented by one of the
additional phase-shifters.

Now, if a phase-shift $\phi$ is required on a single waveguide between layers
of MZI, we follow a new procedure, shown in Fig. 4(b). Adding $+\phi$ to both
phase-shifters in the MZI to the left implements the phase, but a $-\phi$
shift is now required to correct the effect on the waveguide below the
original phase-shift. Now an MZI to the right can be used to implement the
$-\phi$, while requiring a $+\phi$ correction to a lower waveguide. This can
be repeated until the residual phase-shift is moved to the edge of the
circuit, where it can be implemented directly with one of the new phase-
shifters. This relies on the fact that the operation of applying the same
phase-shift to two modes will commute with a beam-splitter operation on those
two modes.

This demonstrates that the new circuit is universal - one can follow the
decomposition given in the previous section and then implement the $Q^{(k)}$
phase-shifts using this method. Two of these operations are already at the
edges of the circuit, and so can trivially be implemented by the edge phase-
shifters. One could also follow the original Clements decomposition intended
for aMZIs, and then implement all of the external phase-shifts required by
this method of absorbing them into surrounding sMZIs and shifting them to the
edge of the circuit.

For each layer of MZIs, there are now $m$ phase-shifters. Since adding the
same phase to all the phase-shifters in a layer only adds a global phase,
there is redundancy here, and any one phase-shifter could be removed from each
layer. The circuit would then have the minimum number of degrees of freedom
required to implement an arbitrary unitary. On the other hand, this is quite a
healthy redundancy to have, since any one phase-shifter in each layer can fail
(e.g. due to fabrication error), and the overall operation of the circuit will
be preserved.

##  V Error Tolerant Designs

The Reck and Clements schemes are provably universal on the assumption that
every beam-splitter in the circuit has an exactly balanced splitting ratio. If
some beam-splitters deviate from this, for instance due to uncertainties in
fabrication, then the MZIs no longer have full tunability in their splitting
ratio, and some unitary transformations become inaccessible to the circuit.
Several alternative designs have been proposed which show improved robustness
to randomized beam-splitters, including adding redundant layers of MZI Burgwal
_et al._ (2017), adding permutations of waveguides between MZI layers Pai  _et
al._ (2019), and the design of Fldzhyan et al. Fldzhyan  _et al._ (2020),
which makes use of alternating layers of beam-splitters and phase-shifters in
an arrangement which does not map onto a network of MZIs. For these designs,
there is generally no known deterministic method of decomposing them into
elementary 2x2 unitaries. Rather the phase settings are optimised to minimise
the distance to some target unitary matrix.

For rectangular meshes of MZI similar to the Clements scheme it is fairly
clear where aMZI can be replaced by sMZI, hence we focus on the Fldzhyan
design here. This design is appealing because the depth of the circuit and the
number of elements is identical to the Clements scheme, while the robustness
to beam-splitter imbalance is improved for Haar-random target unitaries. The
circuit layout is as shown in Fig. 5(a) - here, only four layers of beam-
splitters and phase-shifters are shown, with $2m$ such layers required for
universality. Fig. 5(b) shows a compactified design, where every other layer
of phase-shifters has been moved into the preceding layer. It can be seen that
this design is equivalent to the original one, since where-ever a phase is
required in a removed layer, it can be applied by adjusting phases in the
neighbouring layers. Fig. 5(c) shows how a phase $\phi$ in a removed layer is
implemented by applying $+\phi$ to a subset of phase-shifters in the
neighbouring layer to the left, and $-\phi$ to a subset in the neighbouring
layer to the right. As previously, this relies on the fact that equal phase-
shifts applied to both two modes involved in a beam-splitter can be commuted
to the input or output of the beam-splitter. This remains valid regardless of
the splitting ratio of the beam-splitter.

Figure 5:  (a) The Fldzhyan et al. design, consisting of alternating layers of
beam-splitters and phase-shifters. (b) An equivalent design where every other
layer of phase-shifters have been moved to the preceding layer. (c) Any
required phase $\phi$ in a missing layer can be implemented by applying
$\pm\phi$ shifts in neighbouring layers.

##  VI Conclusions

In summary, we have shown that the aMZI can be replaced by sMZI in a Clements
style rectangular network while retaining universality, providing a
deterministic and efficient method of selecting the phases required to
implement a target unitary. We expect that sMZI can replace aMZI in related
designs making use of a rectangular structure, and have used similar logic to
suggest a more compact but equivalent form of the circuit of Fldzhyan et
al.Fldzhyan  _et al._ (2020). Using sMZIs implies that the circuit length
taken up by phase-shifters is approximately halved, compared to aMZIs. This
can be a significant fraction of the overall circuit length in many integrated
photonic platforms - for example in silicon-on-insulator Shen  _et al._
(2017), silica-on-silicon Mennea  _et al._ (2018), and lithium niobate-on-
insulator Boes  _et al._ (2018), the length of a thermo-optic phase shifter is
comparable to that of a beam-splitter, whereas in silicon nitride the phase-
shifters are relatively long and can be the dominant contribution to the
circuit length Taballione  _et al._ (2020); Arrazola  _et al._ (2021). We
expect that this saving in length will allow larger universal linear optical
circuits to fit on a chip, while reducing the propagation loss, helping to
realise large-scale and complex quantum photonic computation and simulation.

We acknowledge funding from the EPSRC UK Quantum Technologies Programme
(EP/T001062/1) via the Quantum Computing and Simulation hub. Bryn Bell is
supported by a European Commission Marie Skłodowska Curie Individual
Fellowship (FrEQuMP, 846073).

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * O’Brien (2007) J. L. O’Brien, Science 318, 1567 (2007).
  * Kok  _et al._ (2007) P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Rev. Mod. Phys. 79, 135 (2007).
  * Carolan  _et al._ (2015) J. Carolan, C. Harrold, C. Sparrow, E. Martín-López, N. J. Russell, J. W. Silverstone, P. J. Shadbolt, N. Matsuda, M. Oguma, M. Itoh, G. D. Marshall, M. G. Thompson, J. C. F. Matthews, T. Hashimoto, J. L. O’Brien, and A. Laing, Science 349, 711 (2015).
  * Metcalf  _et al._ (2014) B. J. Metcalf, J. B. Spring, P. C. Humphreys, N. Thomas-Peter, M. Barbieri, W. S. Kolthammer, X.-M. Jin, N. K. Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R. Smith, and I. A. Walmsley, Nature Photonics 8, 770 (2014).
  * Qiang  _et al._ (2018) X. Qiang, X. Zhou, J. Wang, C. M. Wilkes, T. Loke, S. O’Gara, L. Kling, G. D. Marshall, R. Santagati, T. C. Ralph, J. B. Wang, J. L. O’Brien, M. G. Thompson, and J. C. F. Matthews, Nature Photonics 12, 534 (2018).
  * Spring  _et al._ (2013) J. B. Spring, B. J. Metcalf, P. C. Humphreys, W. S. Kolthammer, X.-M. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N. K. Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R. Smith, and I. A. Walmsley, Science 339, 798 (2013).
  * Carolan  _et al._ (2014) J. Carolan, J. D. A. Meinecke, P. J. Shadbolt, N. J. Russell, N. Ismail, K. Wörhoff, T. Rudolph, M. G. Thompson, J. L. O’Brien, J. C. F. Matthews, and A. Laing, Nature Photonics 8, 621–626 (2014).
  * Crespi  _et al._ (2013) A. Crespi, R. Osellame, R. Ramponi, D. J. Brod, E. F. Galvão, N. Spagnolo, C. Vitelli, E. Maiorino, P. Mataloni, and F. Sciarrino, Nature Photonics 7, 545 (2013).
  * Bell  _et al._ (2019) B. A. Bell, G. S. Thekkadath, R. Ge, X. Cai, and I. A. Walmsley, Opt. Express 27, 35646 (2019).
  * Wang  _et al._ (2018) J. Wang, S. Paesani, Y. Ding, R. Santagati, P. Skrzypczyk, A. Salavrakos, J. Tura, R. Augusiak, L. Mančinska, D. Bacco, D. Bonneau, J. W. Silverstone, Q. Gong, A. Acín, K. Rottwitt, L. K. Oxenløwe, J. L. O’Brien, A. Laing, and M. G. Thompson, Science 360, 285 (2018).
  * Taballione  _et al._ (2019) C. Taballione, T. A. W. Wolterink, J. Lugani, A. Eckstein, B. A. Bell, R. Grootjans, I. Visscher, D. Geskus, C. G. H. Roeloffzen, J. J. Renema, I. A. Walmsley, P. W. H. Pinkse, and K.-J. Boller, Opt. Express 27, 26842 (2019).
  * Sparrow  _et al._ (2018) C. Sparrow, E. Martín-López, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda, J. L. O’Brien, D. P. Tew, and A. Laing, Nature 557, 660–667 (2018).
  * Harris  _et al._ (2017) N. C. Harris, G. R. Steinbrecher, M. Prabhu, Y. Lahini, J. Mower, D. Bunandar, C. Chen, F. N. C. Wong, T. Baehr-Jones, M. Hochberg, S. Lloyd, and D. Englund, Nature Photonics 11, 447–452 (2017).
  * Shen  _et al._ (2017) Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, and M. Soljačić, Nature Photonics 11, 441–446 (2017).
  * Zhuang  _et al._ (2015) L. Zhuang, C. G. H. Roeloffzen, M. Hoekman, K.-J. Boller, and A. J. Lowery, Optica 2, 854 (2015).
  * Pérez  _et al._ (2017) D. Pérez, I. Gasulla, L. Crudgington, D. J. Thomson, A. Z. Khokhar, K. Li, W. Cao, M. G. Z., and J. Capmany, Nature Communications 8 (2017), 10.1038/s41467-017-00714-1.
  * Reck  _et al._ (1994) M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, Phys. Rev. Lett. 73, 58 (1994).
  * Clements  _et al._ (2016) W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optica 3, 1460 (2016).
  * Miller (2015) D. A. B. Miller, Optica 2, 747 (2015).
  * Ribeiro  _et al._ (2016) A. Ribeiro, A. Ruocco, L. Vanacker, and W. Bogaerts, Optica 3, 1348 (2016).
  * Fldzhyan  _et al._ (2020) S. A. Fldzhyan, M. Y. Saygin, and S. P. Kulik, Opt. Lett. 45, 2632 (2020).
  * Burgwal  _et al._ (2017) R. Burgwal, W. R. Clements, D. H. Smith, J. C. Gates, W. S. Kolthammer, J. J. Renema, and I. A. Walmsley, Opt. Express 25, 28236 (2017).
  * Pai  _et al._ (2019) S. Pai, B. Bartlett, O. Solgaard, and D. A. B. Miller, Phys. Rev. Applied 11, 064044 (2019).
  * Mennea  _et al._ (2018) P. L. Mennea, W. R. Clements, D. H. Smith, J. C. Gates, B. J. Metcalf, R. H. S. Bannerman, R. Burgwal, J. J. Renema, W. S. Kolthammer, I. A. Walmsley, and P. G. R. Smith, Optica 5, 1087 (2018).
  * Boes  _et al._ (2018) A. Boes, B. Corcoran, L. Chang, J. Bowers, and A. Mitchell, Laser & Photonics Reviews 12, 1700256 (2018).
  * Taballione  _et al._ (2020) C. Taballione, R. van der Meer, H. Snijders, P. Hooijschuur, J. Epping, M. de Goede, B. Kassenberg, P. Venderbosch, C. Toebes, H. van den Vlekkert, P. Pinkse, and J. Renema, arxiv:2012.05673 (2020).
  * Arrazola  _et al._ (2021) J. Arrazola, V. Bergholm, K. Brádler, T. R. Bromley, M. J. Collins, I. Dhand, A. Fumagalli, T. Gerrits, A. Goussev, L. G. Helt, J. Hundal, T. Isacsson, R. B. Israel, J. Izaac, S. Jahangiri, R. Janik, N. Killoran, S. P. Kumar, J. Lavoie, A. E. Lita, D. H. Mahler, M. Menotti, B. Morrison, S. W. Nam, L. Neuhaus, H. Y. Qi, N. Quesada, A. Repingon, K. K. Sabapathy, M. Schuld, D. Su, J. Swinarton, A. Száva, K. Tan, P. Tan, V. D. Vaidya, Z. Vernon, Z. Zabaneh, and Y. Zhang, Nature 591, 54 (2021).


D<EMAIL_ADDRESS>Sunouchimaff12014667
<EMAIL_ADDRESS>Yoshiokamaff22014667 [aff1]The author is
with the Design Department, Sapporo city university, Geijutsu-no-mori 1,
Minami-ku, Sapporo, Hokkaido 005–0864, Japan [aff2]The author is with the
Graduate School of Information Science and Technology, Hokkaido University,
Kita 14, Nishi 9, Kita–ku, Sapporo, Hokkaido, 060–0814, Japan

# Diversity-Robust Acoustic Feature Signatures
Based on Multiscale Fractal Dimension
for Similarity Search of Environmental Sounds

###### keywords:

environmental sound analysis, fractals, content-based retrieval, feature
extraction

This paper proposes new acoustic feature signatures based on the multiscale
fractal dimension (MFD), which are robust against the diversity of
environmental sounds, for the content-based similarity search. The diversity
of sound sources and acoustic compositions is a typical feature of
environmental sounds. Several acoustic features have been proposed for
environmental sounds. Among them is the widely-used Mel-Frequency Cepstral
Coefficients (MFCCs), which describes frequency-domain features. However, in
addition to these features in the frequency domain, environmental sounds have
other important features in the time domain with various time scales. In our
previous paper, we proposed enhanced multiscale fractal dimension signature
(EMFD) for environmental sounds. This paper extends EMFD by using the kernel
density estimation method, which results in better performance of the
similarity search tasks. Furthermore, it newly proposes another acoustic
feature signature based on MFD, namely very-long-range multiscale fractal
dimension signature (MFD-VL). The MFD-VL signature describes several features
of the time-varying envelope for long periods of time. The MFD-VL signature
has stability and robustness against background noise and small fluctuations
in the parameters of sound sources, which are produced in field recordings. We
discuss the effectiveness of these signatures in the similarity sound search
by comparing with acoustic features proposed in the DCASE 2018 challenges. Due
to the unique descriptiveness of our proposed signatures, we confirmed the
signatures are effective when they are used with other acoustic features.

##  1 Introduction

Acoustic feature extraction is a basic audio signal processing issue. Acoustic
features are important and necessary for various contexts and applications
related to environmental sound recognition (ESR), such as large-scale content-
based retrieval, auditory scene analysis, visualization, and event detection
for surveillance. During the last decade, handy digital sound recorders have
gained popularity, and at present, not only professional creators, but also
amateurs have started recording environmental sounds and sharing them on web
services such as Freesound [1, 2] and SoundCloud [3]. These sound recordings
are not only appreciated as music works, but also sampled for creating sound
effects, new music works, and live performances in music genres such as
ambient, drone, and electronic [4] [5]. These sound recordings are also
utilized for research to analyze and understand the variety of sound
environments that we live in [6].

###  1.1 Applications Using Acoustic Features for Environmental Sounds

In recent years, the research on ESR for understanding a scene and its context
has received considerable attention [7]. The workshop challenges on Detection
and Classification of Acoustic Scenes and Events (DCASE) have demonstrated
performance evaluations of systems for the detection and classification of
sound events [8]. Based on the best result from Task 1B of DCASE2020, Koutini
et al. evaluated their Receptive Field (RF) regularized CNN model with some
parameter reduction methods. [9].

Classification is a basic application that uses acoustic features. In 2003,
Cowling and Sitte [10] presented a comprehensive comparative study of
classification techniques that use various acoustic features for environmental
sounds. They reported that the test patterns using each of the Mel-Frequency
Cepstral Coefficients (MFCCs) and the continuous wavelet transform achieved
the best recognition performance. In 2009 and 2012, Chu et al. [11] and Mogi
et al. [12] reported that recognition systems that use the Matching-Pursuit-
based acoustic feature as a time-domain feature shows better classification
performance than systems that use the popular MFCCs only as a frequency-domain
feature. In 2013, Bauge et al. [13] proposed a new acoustic feature for
environmental sounds based on the scattering transform. This feature is robust
against frequency transposition.

Content-based retrieval is another basic application that uses acoustic
features. Web-based sound archives such as Freesound and SoundCloud are
becoming popular and the amount of their sound content is increasing. The
online users who utilize these sound archives can share and browse sound
content by means of content-based retrieval. In 2008, Xue et al. [14] proposed
a similarity search system, which employs a cluster-based indexing approach
for environmental sounds. In 2010, Roma et al. [15] proposed a method for the
retrieval of environmental sounds using the general sound-events taxonomy
defined based on the principles of ecological acoustics. Chechik et al. [16]
compared the scalability of several classification methods using MFCCs for a
large-scale content-based sound retrieval. In 2013, Sunouchi and Tanaka [17]
proposed a new acoustic feature signature, namely, the enhanced multiscale
fractal dimension signature (EMFD) and demonstrated the effectiveness of EMFD
for content-based similarity search of environmental sounds.

In recent years, the workshop challenges on DCASE have focused on improving
machine learning methods for ESR and produced high-performance results for
their tasks. Acoustic features are still essential as input data for the
machine learning methods for ESR. Hence finding new acoustic features that can
properly describe the features of environmental sounds is fundamental in
improving the performance of these ESR applications. In addition, by studying
how acoustic features can describe the features of environmental sounds and
affect the performance of ESR tasks, we can develop an understanding of how we
are listening to environmental sounds.

###  1.2 Acoustic Feature Extraction for ESR

The environmental sounds outside a recording studio are produced by action and
movement. We can identify things by listening to their acoustic properties,
which are the results of the sound production process. However, environmental
sound signals of the same type cannot be physically identical to each other
due to the difference in their production processes. Furthermore, the
different sound signals generated by simultaneous events are mixed with each
other, which makes the properties of each sound source obscured [18].

Various acoustic features have been proposed for content-based audio
retrieval. The feature selection is an important process for ESR [19], [20].
Cepstral features that include MFCCs and their first and second derivatives
(MFCCs$\Delta$ and MFCCs$\Delta\Delta$) are widely used as frequency-domain
acoustic features. MP-based acoustic feature has been proposed as one of the
useful time-domain features for ESR [11],[12] [21].

Recent researches have focused on the evaluation of time-domain features of
environmental sounds. For ESR, we need acoustic features that describe the
non-stationary characteristics of target sounds as a time-domain feature and
are robust against the diversity of environmental sounds [7]. We have
recognized there may be three main causes of the diversity of environmental
sounds.

  * D1)

Small fluctuations of sound source parameters, such as carrier signal
frequency, due to the individuality of the sound source.

  * D2)

Background noises that the person who recorded the target sound did not expect
to record.

  * D3)

Mixed composition of different types of sound sources.

For the third cause D3, it is necessary to apply, for example, independent
component analysis or non-negative matrix factorization to the sound signal
before the feature extraction process [12],[22]. In this study, we focus on
the extraction of new acoustic feature signatures that are robust against the
diversities caused by D1 and D2.

###  1.3 Problems of EMFD Signature and Their Solutions

In our previous work [17], we proposed an EMFD signature that can describe
both the frequency-domain features and time-domain features of target sounds.
The EMFD signature is a feature vector, which consists of the time-varying
multiscale fractal dimension (MFD) values. We demonstrated that EMFD improves
the performance of similarity search by supplementing MFCCs. Unfortunately, it
is found that EMFD includes error values that depend on the number of analysis
windows and the histogram’s bin size used for computing its histogram.
Furthermore, EMFD seems to be oversensitive while discriminating the features
of environmental sounds, and may lack robustness against the diversity of
environmental sounds.

In this study, we extend the EMFD signature by improving the process of
computing its histogram using the kernel density estimation method. By
optimizing the bandwidth parameter used for kernel density estimation, the
histogram of the enhanced multiscale fractal dimension using kernel density
estimation signature (EMFD-KDE) becomes sufficiently smooth and robust against
the diversity of environmental sounds as an acoustic feature signature. In
Sect. 2, we present the basic theory and characteristics of EMFD. In Sect. 3,
we propose a method to compute the EMFD-KDE signature. In Sect. 5 and 6, we
demonstrate that EMFD-KDE improves the performance of the similarity search
system.

Furthermore, we enhance the idea of EMFD and propose a new acoustic feature
signature, namely very-long-range multiscale fractal dimension signature (MFD-
VL). The environmental sounds have important acoustic features over a long
time period. However, EMFD cannot describe the time-domain feature for time
periods longer than 10 ms. In Sect. 4, we propose a method to compute the MFD-
VL signature. In addition, we demonstrate that MFD-VL can describe the
features of the time varying envelope for long periods of time, and that it
has the robustness against the diversity causes D1 and D2.

In Sect. 7, we conclude that the proposed feature signatures of EMFD-KDE and
MFD-VL solve the problems of EMFD and are effective when they are used with
other acoustic features, including MFCCs and acoustic features proposed in the
DCASE 2018 challenges.

##  2 Basic Theory of Enhanced Multiscale Fractal Dimension Signature

Mandelbrot, who advocated a concept of fractal in 1975 for the first time,
demonstrated that some structures in nature could be modeled well by the
theory of fractals [23]. One of the most important characteristics of fractals
is that they have self-similarity properties at multiple scales. In the field
of acoustics, Voss and Clarke analyzed the power spectrum of fluctuating
physical variables including frequency, loudness and pitch in music and
speech[24]. They obtained the $1/f^{\gamma}~{}(0.5\lesssim\gamma\lesssim 1.5)$
aspects in the power spectrum of each variable against the frequency of a
signal passed through a low-pass filter having a range 0 Hz – 1 Hz. Hsu [25]
compared the fractal geometry of classical music works, and found that there
is a relation, defined by the theory of fractals, between the interval of
successive notes and their frequency of occurrence.

###  2.1 Multiscale Fractal Dimension

A fractal dimension is an index value that can describe the characteristics of
a fractal by quantifying their complexity as a ratio of the change in detail
to the change in scale. Acoustic features based on the fractal dimension have
been proposed and utilized for various practical applications in the fields
such as acoustics, music analysis, image analysis, physics, physiology, and
neuroscience. Maragos et al. [26][27] proposed the short-time fractal
dimension of speech signals as an acoustic feature and used it for speech
segmentation and sound classification. Zlatintsi and Maragos [28][29] proposed
a multiscale fractal dimension (MFD) profile as a short-time descriptor and
found that this descriptor can discriminate several aspects among different
musical instruments.

###  2.2 Steps to Compute the EMFD Signature

In our previous work [17], we developed EMFD as a feature signature of
environmental sounds for a similarity search system. The EMFD is computed as
follows.

####  2.2.1 Preprocessing Target Sounds

The maximum amplitude of each target sound that is to be analyzed must be
first normalized to -0.1 db. They are converted to the standard format with
the following specifications: sampling rate of 44.1 kHz and bit depth of 16
bits.

####  2.2.2 Computing the Area of Minkowski Sausage

The fractal dimension of a sound signal can be computed based on the
Minkowski-Bouligand dimension. A covering area can be drawn by moving a unit
disk of radius $r$ along the curve of the waveform. This covering area is
called a Minkowski Sausage. The center of the unit disk should be at any
position on the curve of waveform and the width of Minkowski Sausage becomes
$2r$. Figure 1 shows the Minkowski Sausage obtained by moving the unit disk
along the waveform. To compute the area of the Minkowski Sausage of a discrete
sound signal, the unit disk vector $C(r)$ is defined as Eq. (1), where $r$
denotes the radius of the unit disk and $i$ denotes the discrete position on
the horizon. Figure 2 shows how the model of the unit disk is built. The
vertical distance from the center to the top of unit disk at each horizontal
position is denoted by the unit disk vector $C(r)$. Let $n$ be the sampling
position, $r$ the radius of the unit disk, $p$ the discrete position of the
unit disk, and $\mathrm{sig}(x)$ the amplitude of sound signal at each
sampling position $x$. The area of the Minkowski Sausage $area(n,r)$ at each
sampling position $n$ is computed as Eq. (2).

Figure 1: A sound waveform and a Minkowski sausage Figure 2: Mesh-
Approximation of a unit disk

$C(r)=\left\\{\mathrm{floor}\left(\sqrt{2ri-i^{2}}\right)\Bigm{|}0\leq i\leq
2r,i\in\mathbb{Z}\right\\}$ (1)

$\mathrm{area}(n,r)\\!=\\!\max\limits_{\begin{subarray}{c}0\leq p\leq 2r\\\
p\in\mathbb{Z}\end{subarray}}\\!\left(\\!\mathrm{sig}\left(n\\!-\\!r\\!+\\!p\right)\\!+\\!\mathrm{floor}\\!\left(\\!\sqrt{2rp\\!-\\!p^{2}}\\!\right)\\!\right)\\\
\\!-\\!\min\limits_{\begin{subarray}{c}0\leq p\leq 2r\\\
p\in\mathbb{Z}\end{subarray}}\\!\left(\\!\mathrm{sig}\left(n\\!-\\!r\\!+\\!p\right)\\!-\\!\mathrm{floor}\\!\left(\\!\sqrt{2rp\\!-\\!p^{2}}\\!\right)\\!\right)$
(2)

####  2.2.3 Definition of Multiscale Fractal Dimension

The MFD values are computed for each analysis window whose period is 50ms. Let
$A(r)$ be the area of the Minkowski Sausage drawn by the unit disk of radius
$r$ in each analysis window. The MFD of each analysis window is defined by Eq.
(3). The minimum radius (r=1) corresponds to the sampling period of the signal
(1/44.1 ms) and the range of $r$ from 1 to 132 corresponds to the range of the
time scales from 1/44.1 to 3 ms.

$\mathop{MFD}\\!=\\!\left\\{2\\!-\\!\frac{\log\left(A(r\\!+\\!1)/A(r)\right)}{\log\left((r\\!+\\!1)/r\right)}\Biggm{|}1\\!\leq\\!r\\!\leq\\!132,~{}r\\!\in\\!\mathbb{Z}\right\\}$
(3)

$\mathop{MFD_{enhanced}(x)}=2-\frac{\log\left(A(r(x+1))/A(r(x))\right)}{\log\left(r(x+1)/r(x)\right)}\mathrm{,where}~{}r(x)=\mathrm{round}(1.4^{x})$
(4)

$\mathop{AW(sound,~{}period)}=\left\\{0,period,...,\mathrm{floor}\left(\frac{\mbox{the
length of}~{}sound}{period}-1\right)\times period\right\\}$ (5)

$\mathop{FAW(dbin,~{}rbin)}=\bigl{\\{}~{}t~{}\bigm{|}~{}t\in{AW(sound,~{}50)},\\\
1+(dbin-1)/32\leq MFD_{enhanced}(rbin)~{}\mbox{of analysis
window}~{}t<1+dbin/32\bigr{\\}}$ (6)

$\mathop{EMFD(sound)}=\Biggl{\\{}\frac{\mathrm{card}\left(FAW(dbin,rbin)\right)}{\mathrm{card}\left(AW(sound,~{}50)\right)}\Biggm{|}1\leq
dbin\leq 32,~{}dbin\in\mathbb{Z},~{}1\leq rbin\leq
16,~{}rbin\in\mathbb{Z}\Biggr{\\}}$ (7)

####  2.2.4 Definition of the EMFD Signature

In our previous work [17], we found that MFD has informative values for unit
disks larger than the disk with a radius of 3 ms $(r=132)$. The maximum radius
of the unit disk was extended to 218, which corresponds to 5 ms (1/10 of the
period of analysis window), and the discrete values of the unit disk were
modified to have exponential values. The enhanced MFD value at the $x$-th
discrete value of the unit disk is defined as Eq. (4). The enhanced MFD values
are computed for each analysis window. The EMFD signature is then defined as
the two-dimensional histogram $(16\times 32)$ of the time-varying enhanced
MFD. Let $period$ be the period (ms) of the analysis window and $sound$ be the
target sound. The set of analysis windows of the target sound is defined as
Eq. (5). Let $rbin$ be the bins that correspond to a series of 16 numbers used
to define the different radius of the unit disk for computing the enhanced
MFD, and $dbin$ be the bins that correspond to a series of the 32 small
intervals into which the range of fractal dimension is divided. The set of
analysis windows whose enhanced MFD values fall into the bin $(dbin,rbin)$ is
defined by Eq. (6). The set of values in each bin of the EMFD histogram is
defined by Eq. (7). Figure 3 shows the histogram that visualizes the EMFD
signature of a cuckoo sound $S_{cuckoo}$. The length of $S_{cuckoo}$ is
21.08s.

Figure 3: Visualization of the EMFD signature of a cuckoo sound $S_{cuckoo}$

###  2.3 Known Characteristics of the EMFD signature

Zlatintsi and Maragos [29] concluded that the MFD profiles are useful for
quantifying the multiscale complexity and fragmentation of the different
states of the instrument sound waveforms. In our previous work [17], we
confirmed that EMFD describes the frequency-domain features and several other
effective features of environmental sounds that MFCCs cannot describe.
Furthermore, we confirmed that the EMFD signature has robustness against
changes in volume levels and phase shifting of sound signals in the analysis
window.

###  2.4 Problems of the EMFD signature

The EMFD includes error values that depend on the number of analysis windows
and the histogram’s bin size, which is defined for computing its histogram, as
shown in Eq. (7). In Sect. 3, we extend the existing EMFD by employing a
kernel density estimation method to solve this problem.

Another problem of EMFD is that it cannot describe time-domain features for a
time period longer than 10 ms, although environmental sounds have important
acoustic features over a long time period. To solve this problem, we introduce
MFD-VL signature in Sect. 4 as a newly developed time-domain acoustic feature.

##  3 Extending EMFD Employing a Kernel Density Estimation Method

As mentioned in Subsect. 2.2, EMFD is computed as the two-dimensional
histogram of time-varying enhanced MFD values. Let $NFAW(bin)$ be the number
of analysis windows whose enhanced MFD values fall into the bin, and
$NAW(sound)$ be the total number of analysis windows of the target sound
$sound$, as defined in Eq. (8). The EMFD value of each bin
$EMFD(sound,~{}bin)$ is computed as Eq. (9). This method has the following two
problems.

$\mathop{NAW(sound)}=\mathrm{card}\bigl{(}AW(sound,~{}50)\bigr{)}$ (8)

$\mathop{EMFD(sound,~{}bin)}=\frac{NFAW(bin)}{NAW(sound)}$ (9)

The first problem is that the value of each bin necessarily includes an error,
because the value can be one of the discrete values given by the density of
analysis windows. In particular, a lower number of analysis windows for the
target sound increases the errors. Ideally, the EMFD histogram should be a
continuous probability distribution of the time-varying enhanced MFD values,
regardless of the number of analysis windows.

The second problem is that EMFD computed by the existing method is
oversensitive to discriminate the features of environmental sounds. The tones
and frequencies of each environmental sound may often vary, depending on the
recording conditions and individual characteristics of the sources that
generate this sound, even if person try to record the same type of
environmental sounds in the same way. Therefore, a feature signature of
environmental sounds should have robustness against the diversity of
environmental sounds caused by D1 defined in Subsect. 1.2.

To solve these problems, we introduce the kernel density estimation method to
compute the EMFD histogram.

###  3.1 Definition of the EMFD-KDE Signature

The kernel density estimation method is employed to compute the probability
distribution of the enhanced MFD values at each radius of the unit disk. The
values of each bin of EMFD-KDE are defined as shown in Eq. (4), Eq. (8), Eq.
(10), Eq. (11), and Eq. (12), where $K(\cdot)$ is the kernel function which is
a Gaussian function, and $h$ in Eq. (11) is the smoothing parameter called
bandwidth.

$\mathop{K(x)}=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^{2}}$ (10)

$\mathop{f_{emfd\mbox{\scriptsize-}kde}(dbin_{val},rbin)}=\frac{1}{NAW~{}h}\\\
\times\sum^{NAW}K\left(\frac{dbin_{val}-MFD_{enhanced}(rbin)}{h}\right)$ (11)

$EMFD\mathchar 45\relax
KDE\\!=\\!\Biggl{\\{}f_{emfd\mbox{\scriptsize-}kde}\left(1\\!+\\!\frac{dbin\\!-\\!0.5}{32},rbin\right)\\\
\Biggm{|}1\\!\leq\\!dbin\\!\leq\\!32,~{}dbin\\!\in\\!\mathbb{Z},1\\!\leq\\!rbin\\!\leq\\!16,~{}rbin\\!\in\\!\mathbb{Z}~{}\Biggr{\\}}$
(12)

###  3.2 Optimization of the Bandwidth for Kernel Density Estimation

The bandwidth $h$ is a smoothing parameter, which is usually determined by the
trade-off between the number of data samples and their standard deviation. Let
$n$ be the number of data samples and $\sigma$ be the standard deviation of
the data samples. The bandwidth $h$ of a Gaussian kernel density estimator is
given by the normal reference rule defined by Eq. (13). The normal reference
rule is most commonly used to determine the bandwidth [30].

$\mathop{h}=\left(\frac{4\sigma^{5}}{3n}\right)^{\frac{1}{5}}\approx
1.06\sigma n^{-\frac{1}{5}}$ (13)

We define the bandwidth $h_{rbin}(\alpha)$, which is optimized for each radius
of the unit disk, as Eq. (14), Eq. (15), and Eq. (16), where $avg$ is the
arithmetic mean of the enhanced MFD values of each analysis window at $rbin$,
and $\sigma_{rbin}$ is the standard deviation of the enhanced MFD values at
$rbin$. The smoothing parameter $\alpha$ in Eq. (16) is a constant. Through
the experiments with different values of $\alpha$, the constant value is
determined so that the best result for the target task is obtained.

$avg=\frac{1}{NAW}\sum^{NAW}{MFD_{enhanced}(rbin)}$ (14)

$\sigma_{rbin}\\!=\\!\sqrt{\frac{1}{NAW}\\!\sum^{NAW}\\!\\!\left(\\!MFD_{enhanced}(rbin)\\!-\\!avg\right)^{2}}$
(15)

$\mathop{h_{rbin}(\alpha)}=1.06\sigma_{rbin}NAW^{-\frac{1}{5}}\alpha$ (16)

Figure 4 shows the 3D histogram visualizing the EMFD signature of the cuckoo
sound $S_{cuckoo}$. Figure 5 shows the 3D histogram visualizing the EMFD-KDE
signature of the same cuckoo sound $S_{cuckoo}$. The 3D histogram of the EMFD-
KDE signature is much smoother than that of the EMFD signature. At each radius
of the unit disk, the larger standard deviation of the enhanced MFD values
$\sigma_{rbin}$ results in the smoother histogram.

Figure 4: The 3D histogram image that visualizes the existing EMFD signature
of the cuckoo sound $S_{cuckoo}$ Figure 5: The 3D histogram image that
visualizes the EMFD-KDE signature of the cuckoo sound $S_{cuckoo}$. The
bandwidth is $h_{rbin}(32)$.

##  4 Very Long Range Multiscale Fractal Dimension Signature

The environmental sounds have important acoustic features with varying time
periods. However, EMFD cannot describe the time-domain features for time
periods longer than 10 ms. To solve this problem, we propose a new acoustic
feature signature for environmental sounds, based on the multiscale fractal
dimension. This feature signature is called very-long-range multiscale fractal
dimension signature (MFD-VL). The basic idea of MFD-VL is to extend the size
range of the unit figure to consider the larger ones, which are used to
compute the area of the Minkowski Sausage. In this section, we define the
method to compute the MFD-VL signature and demonstrate its characteristics.

###  4.1 Definition of the MFD-VL Signature

The multiscale fractal dimension values of MFD-VL are computed for an entire
target sound, and not for each fixed-length analysis window of the target
sound. A unit square, instead of a unit disk, is used to compute the area of
the Minkowski Sausage for MFD-VL. Figure 6 shows the Minkowski Sausage
obtained by moving the unit square along the waveform. The method using the
unit square is much faster than the one using the unit disk. Let $n$ denote
the sampling position, $r$ the half side-length of the unit square, $p$ the
discrete position of the unit square, and $sig(x)$ the amplitude value of the
sound signal at each sampling position $x$. The area of the Minkowski Sausage
$area_{sq}(n,r)$ at each sampling position $n$ is computed as Eq. (17).

Figure 6: A sound waveform and a Minkowski Sausage obtained by moving the unit
square

$area_{sq}(n,r)=\max\limits_{\begin{subarray}{c}0\leq p\leq 2r\\\
p\in\mathbb{Z}\end{subarray}}\left(\mathrm{sig}\left(n-r+p\right)+r\right)\\\
-\min\limits_{\begin{subarray}{c}0\leq p\leq 2r\\\
p\in\mathbb{Z}\end{subarray}}\left(\mathrm{sig}\left(n-r+p\right)-r\right)$
(17)

The half side-length $r$ of the unit square for each scale was defined as Eq.
(18), where $sf$ is the sampling frequency of the target sound. In this study,
$sf$ is 44100 Hz. Let $A_{sq}(r)$ denote the area of the Minkowski Sausage
obtained for an entire target sound by moving the unit square whose side
length is $2r$. The MFD-VL signature is defined as Eq. (19). The MFD-VL
signature is a feature vector that contains 10 elements.

$r(x)=\mathrm{round}\left(sf\times 2^{-{\frac{x+2}{2}}}\right)$ (18)

$MFD\mathchar 45\relax
VL\\!=\\!\biggl{\\{}2\\!-\\!\frac{\log\left(A_{sq}(r(x))/A_{sq}(r(x\\!+\\!1))\right)}{\log\left(r(x)/r(x\\!+\\!1)\right)}\\\
\biggm{|}0\\!\leq\\!x\\!\leq\\!9,~{}x\\!\in\\!\mathbb{Z}\biggr{\\}}$ (19)

###  4.2 Basic Characteristics of the MFD-VL Signature

We found several basic characteristics of MFD-VL through the experiments using
test sounds.

####  4.2.1 MFD-VL’s Descriptiveness of the Beats of Single Sine Waves

The MFD-VL signature is expected to describe the acoustic features over very
long time-periods. We found that MFD-VL can discriminate frequencies of
amplitude envelopes between 22.6 Hz and 1 Hz. The range of the wavelength of
the amplitude envelopes corresponds to the range of the side length of the
unit square between 0.044 s and 1 s. Let $f_{beat}$ denote the frequency of
the beats and $f_{content}$ denote the frequency of single sine waves inside
the amplitude envelopes. The set of test sounds $SS_{t1}$ is defined as Eq.
(20), Eq. (21), and Eq. (22). Each test sound is filtered by the pink noise
filter function $f_{pn}$ Eq. (20). A sound that is artificially synthesized
using pure tones usually has distinct or sparse spectra. This kind of sounds
may cause numerical instabilities while calculating their acoustic features.
To solve this problem, the pink noise filter function $f_{pn}$ is used to add
a background pink noise, which is defined as Eq. (20), where $sig$ is an input
signal and $Noise_{pink}$ is a background pink noise whose maximum amplitude
is normalized to -0.1 db. The signal-to-noise ratio is 24 db. The pink noise,
known as $1/f$ noise, is a signal whose power spectral density is inversely
proportional to the signal frequency. The pink noise signal is known to widely
exist in the natural world. The frequency components below 40 Hz contained in
the pink noise are cut off by using a low cut filter before the amplitude
normalization because the components with lower frequencies cannot be recorded
nor played using common microphones and speakers.

$\mathop{f_{pn}\left(sig\right)}=\frac{15}{16}sig+\frac{1}{16}Noise_{pink}$
(20)

$s_{t1}(f_{beat},f_{content})\\!=\\!f_{pn}\bigl{(}\cos(\pi{f_{beat}}~{}t)\sin(2\pi{f_{content}}~{}t)\bigr{)}$
(21)

$\begin{split}SS_{t1}=&\bigl{\\{}s_{t1}(f_{beat},f_{content})\\\
&\quad\bigm{|}f_{content}=440,f_{beat}\in\\{0.5,1,2,4,8,16,32\\}\bigr{\\}}\end{split}$
(22)

Figure 7 shows the line charts of the MFD-VL values of single sine waves of
frequency 440 Hz, which is filtered by Eq. (20), and those of the test sounds
$SS_{t1}$. In Fig. 7, the frequencies of the beats are indicated by the
troughs of the line chart, in which the side length of the unit square is less
than the wavelengths of the beats. This characteristic can be understood
morphologically as shown in Fig. 8. When the side length of the unit square is
more than the wavelength of the beat, the area of the Minkowski Sausage
becomes almost the same as that of a single sine wave without beats. When the
side length of the unit square is less than the wavelengths of the beats, the
shorter side length of the unit square results in a smaller area of Minkowski
Sausage.

Figure 7: Line charts of the MFD-VL signatures of single sine waves of
frequency 440Hz, with and without a beat. The beat frequencies are 0.5Hz, 1Hz,
2Hz, 4Hz, 8Hz, 16Hz, and 32Hz. Figure 8: Minkowski Sausage obtained by moving
the unit square whose side length is less than the wavelength of the beat.

####  4.2.2 MFD-VL’s Descriptiveness of Amplitude Envelope Shapes

Here we analyze some other characteristics of the MFD-VL that are related to
the descriptiveness of the amplitude envelope shapes. Let $f_{pulse}$ be the
frequency of the rectangular pulse waves and $w_{pulse}$ be the ratio of the
rectangular pulse width to the wavelength of the rectangular pulse waves. The
rectangular pulse function $rect({f_{pulse}},~{}{w_{pulse}},~{}t)$ for
generating the amplitude envelopes is defined as Eq. (23). Let $f_{content}$
be the frequency of a single sine wave inside the amplitude envelopes
generated by the rectangular pulse function. The test sound $s_{t2}$ is
defined as Eq. (23) and Eq. (24). The test sound is filtered by the pink noise
filter function $f_{pn}$ Eq. (20).

$rect({f_{pulse}},~{}{w_{pulse}},~{}t)\\!=\\!\begin{cases}1\\!&\\!\left(t\bmod\frac{1}{f_{pulse}}<=\frac{w_{pulse}}{f_{pulse}}\right)\\\
0\\!&\\!(otherwise)\end{cases}$ (23)

$s_{t2}(f_{pulse},w_{pulse},f_{content})\\\
=f_{pn}\bigl{(}rect({f_{pulse}},~{}{w_{pulse}},~{}t)\sin(2\pi{f_{content}}~{}t)\bigr{)}$
(24)

We define the set of test sounds $SS_{t2}$ as Eq. (25). The line charts of
MFD-VL of the single sine wave of frequency 440 Hz, which is filtered by Eq.
(20), and $SS_{t2}$ are showed in Fig. 9. Here, we compare the line charts of
MFD-VL for $s_{t1}$ and those for $s_{t2}$. This comparison shows that the
bottom of the line chart trough of $s_{t2}$ for $f_{pulse}=2$ is deeper than
that of $s_{t1}$ for $f_{beat}=2$, and that the bottom of the line chart
trough of $s_{t2}$ for $f_{pulse}=4$ is also deeper than that of $s_{t1}$ for
$f_{beat}=4$.

$\begin{split}SS_{t2}\\!=&\bigl{\\{}s_{t1}(f_{beat},f_{content})\\!\bigm{|}\\!f_{content}\\!=\\!440,f_{beat}\\!\in\\!\\{2,4\\}\bigr{\\}}\\\
&\cup\bigl{\\{}s_{t2}(f_{pulse},w_{pulse},f_{content})\\\
&\quad\bigm{|}f_{content}\\!=\\!440,f_{pulse}\\!\in\\!\\{2,4\\},w_{pulse}\\!=\\!0.5\bigr{\\}}\end{split}$
(25)

Figure 9: Line charts of the MFD-VL signatures for single sine waves having a
frequency of 440 Hz; those masked by the cosine function, and those masked by
the rectangular pulse function.

For another comparison, we define a set of test sounds $SS_{t3}$ as Eq. (26).
The set of test sounds $SS_{t3}$ contains single sine waves having a frequency
of 440 Hz masked by the rectangular pulse functions with various widths
${w_{pulse}}$. In Fig. 10, we compare the line charts of MFD-VL of single sine
wave having a frequency of 440 Hz, $s_{t1}(f_{content}=440,f_{beat}=4)$, and
$SS_{t3}$. This comparison shows that the narrower width of the amplitude
envelopes made by the rectangular pulse function results in deeper troughs in
the line chart.

$\begin{split}SS_{t3}=&\bigl{\\{}s_{t2}(f_{pulse},w_{pulse},f_{content})\\\
&\bigm{|}{f_{content}}\\!=\\!440,{f_{pulse}}\\!=\\!4,w_{pulse}\\!\in\\!\\{0.2,0.5,0.8\\}\bigr{\\}}\end{split}$
(26)

Figure 10: Line charts of the MFD-VL signatures of single sine waves having a
frequency of 440 Hz; those masked by the cosine function, and those masked by
the rectangular pulse function having various widths.

###  4.3 MFD-VL’s Descriptiveness of Amplitude Envelopes of Simulated
Environmental Sounds

Environmental sounds such as chirping of insects and birds, and water streams
have amplitude-modulated waveforms. It is well-known that some type of
environmental sounds can be simulated by the granular synthesis technique
[31]. The granular synthesis technique splits a sound signal into small pieces
called grains. Each grain has an envelope that contains an actual sonic
content. If the content signal is a simple sinusoid, the synthesized sound
signal is considered to be the same as a sound signal produced using the
amplitude-modulation technique.

For example, the sound signal of the chirping of a cricket has a carrier wave
with a frequency of around 5 kHz, and the carrier wave is modulated into
trains of syllable chirps whose rate is near 30 Hz [32]. The sound of a
cricket can be roughly simulated by the granular synthesis with an envelope
created by the Hann window, which is defined as Eq. (27). Let $f_{c}$ be the
frequency of the carrier wave and $f_{e}$ be the frequency of the envelope
rate. The model of the chirping of a cricket $sim_{cricket}$ can be defined,
as shown in Eq. (28) and Eq. (29). The sound signal of the chirping of a
Gryllus bimaculatus, which is a species of cricket living in Japan, can be
simulated as $\mathop{sim_{cricket}\left(t,5800,30\right)}$. It is known that
the chirping of a Gryllus bimaculatus sometimes has continuous 3-syllable
chirps that are repeated three times per second. Let $f_{c}$ denote the
frequency of the carrier wave, $sn$ be the number of continuous syllable
chirps, $f_{rep}$ be the number of repetitions per second of the set of
continuous syllable chirps, and $f_{e}$ be the frequency of continuous chirps.
Then, the model of chirping of a cricket $sim_{cricket2}$ is defined using Eq.
(30) and Eq. (31). Figure 11 shows the sound waveform of signal
$sim_{cricket2}(t,5800,3,2.73,30)$.

$\mathop{\omega\left(x\right)=\frac{1}{2}\left(1-\cos\left(2{\pi}x\right)\right)}$
(27)

$\mathop{T_{s}\left(t,f\right)}=\min\left(t\bmod\frac{1}{f},\frac{1}{1.1f}\right)$
(28)

$\mathop{sim_{cricket}\left(t,f_{c},f_{e}\right)}=\sin\left(2{\pi}f_{c}t\right)\mathop{\omega\left(1.1f_{e}T_{s}\left(t,f_{e}\right)\right)}$
(29)

$rect_{cricket}({f_{e}},~{}sn,~{}f_{rep},~{}t)\\!=\\!\begin{cases}1\\!&\\!\left(t\bmod\frac{1}{f_{rep}}\\!\leq\\!\frac{sn}{f_{e}}\right)\\\
0.05\\!&\\!(otherwise)\end{cases}$ (30)

$\mathop{sim_{cricket2}\left(t,f_{c},sn,f_{rep},f_{e}\right)}\\\
\\!=\\!rect_{cricket}({f_{e}},~{}sn,~{}f_{rep},~{}t)\mathop{sim_{cricket}\left(t,f_{c},f_{e}\right)}$
(31)

Figure 11: Sound waveform of the signal $sim_{cricket2}(t,5800,3,2.73,30)$

We define the set of test sounds $SS_{t4}$ as Eq. (32). For preventing sparse
spectra of the synthesized sound, each test sound is filtered by the pink
noise filter function $f_{pn}$ Eq. (20), which is defined in Subsubsect.
4.2.1. Figure 12 shows the line charts of MFD-VL of $SS_{t4}$. Figure 13 shows
the 13-coefficients MFCCs values (MFCC13) of $SS_{t4}$. The SPTK toolkit [33]
was used to compute MFCC13. The MFD-VL line charts of the sound set between
$sim_{cricket}$ (group1) and $sim_{cricket2}$ (group2) are clearly different.
These results indicate that the MFD-VL signature can clearly discriminate the
two amplitude envelope signals.

$\begin{split}SS_{t4}\\!=\\!&\bigl{\\{}f_{pn}\bigl{(}sim_{cricket}(t,f,30)\bigr{)}\bigm{|}f\\!\in\\!\\{5300,5800,6300\\}\bigr{\\}}\\\
&\cup\bigl{\\{}f_{pn}\bigl{(}sim_{cricket2}(t,f,3,2.73,30)\bigr{)}\\\
&\quad\bigm{|}f\\!\in\\!\\{5300,5800,6300\\}\bigr{\\}}\end{split}$ (32)

Figure 12: Line charts of the MFD-VL signatures of simulated cricket sounds.
These show the MFD-VL’s robustness against fluctuations of a carrier signal’s
frequency. Figure 13: MFCC13 values of simulated cricket sounds. These show
the descriptiveness of varying carrier signals’ frequency.

###  4.4 MFD-VL’s Robustness against Fluctuations of a Carrier Signal’s
Frequency

We have assumed that one of the main causes of the diversity of environmental
sounds is small fluctuations of sound source parameters, such as carrier
signal frequency, due to the individuality of the sound source (D1). For
example, the frequency of the human voice varies with the size and form of the
individual’s vocal tract which are the sound source parameters with
fluctuations. However, we can recognize the sound recordings of talking by
children, adults, and aged persons as the human voice in the same manner. The
set of test sounds $SS_{t4}$ contains the synthesized sound signals of the
chirping of a cricket with different frequencies of a carrier wave. We can
evaluate the descriptiveness of acoustic features for the small fluctuations
of sound source parameters by comparing the feature vectors of test sounds
contained in $SS_{t4}$.

Figure 12 shows the robustness of MFD-VL against the fluctuation of the
carrier signal’s frequency contained in the amplitude envelopes. We defined
the discrimination rate of both MFD-VL and MFCC13 for the quantitative
evaluation. Let $se(ts,n)$ denote the $n$-th element value of MFD-VL signature
of the test sound clip $ts$. The MFD-VL discrimination rate $DR_{mfdvl}$ for
comparing two test sounds A and B is calculated using Eq. (33). The
$RANGE_{mfdvl}$ in Eq. (33) denotes the width of the range of $se(ts,x)$,
which may take a value between 1 and 2. Therefore, $RANGE_{mfdvl}$ is
calculated as 1 $(=2-1)$. The MFCC13 discrimination rate $DR_{mfcc13}$ for
comparing two test sounds A and B was defined as Eq. (34) where $se2(ts,n)$
was defined as the $n$-th dimension value of the MFCC13 vector of the test
sound clip $ts$. The constant $RANGE_{mfcc13}$ was set to 56.3. This number
was determined from the difference between the minimum and the maximum values
of MFCC13 of 3000 environmental sounds in a dataset for the experimental
evaluation which is described in Sect. 5. The sounds in the dataset were
collected from the Freesound project [2].

$DR_{mfdvl}=\sqrt{\frac{1}{10}\sum_{n=1}^{10}\left(\frac{se(A,n)-se(B,n)}{RANGE_{mfdvl}}\right)^{2}}$
(33)

$\mathop{DR_{mfcc13}}=\sqrt{\frac{1}{13}\sum_{n=1}^{13}\left(\frac{se2(A,n)-se2(B,n)}{RANGE_{mfcc13}}\right)^{2}}$
(34)

Table 1 lists the comparisons of the two discrimination rates $DR_{mfdvl}$ and
$DR_{mfcc}$ to distinguish two test sounds with various carrier signal
frequencies which are contained in $SS_{t4}$. This comparison shows that to
distinguish any two test sounds with various carrier signal frequencies, the
values of $DR_{mfdvl}$ are not more than 0.22%, and those of $DR_{mfcc}$ are
not less than 2.90% and no more than 16.15%. These results indicate that the
MFD-VL signature shows robustness against the fluctuation of the carrier
signal frequency, which is the diversity cause D1 defined in Subsect. 1.2. The
MFD-VL signature can clearly discriminate the characteristics of the amplitude
envelopes with its robustness against the fluctuation of the carrier signal
frequency.

Table 1: Comparison of discrimination rates of varying carrier signals’
frequencies between the MFD-VL signature and MFCC13

###  4.5 MFD-VL’s Robustness against Noises

To evaluate the robustness of MFD-VL against noises, we use test sounds
generated from the simulated cricket sound and pink noise signals. Let
$S_{cricket2}$ denote $sim_{cricket2}(t,~{}5800,~{}3,~{}2.73,~{}30)$,
$Noise_{pink}$ denote a pink noise signal whose maximum amplitude is
normalized to -0.1 db, and $\beta$ be a constant number to control the signal-
to-noise ratio (SNR). Then, a set of test sounds $SS_{t5}$ is calculated using
Eq. (35). The set of constant numbers $\beta\in\\{0.67,0.8,0.89,0.94,0.97\\}$
corresponds to the set of SNRs $\\{6db,12db,18db,24db,30db\\}$ of $SS_{t5}$.
The frequency components below 40 Hz contained in the pink noise are cut off
by using a low cut filter before the amplitude normalization because the
components with lower frequencies cannot be recorded nor played using common
microphones and speakers.

$SS_{t5}=\bigl{\\{}\beta~{}S_{cricket2}+(1-\beta)~{}Noise_{pink}\\\
\bigm{|}\beta\in\\{0.67,0.8,0.89,0.94,0.97\\}\bigr{\\}}$ (35)

Figure 14 shows the line charts of the discrimination rates $DR_{mfdvl}$ and
$DR_{mfcc}$ for distinguishing the simulated cricket sound $S_{cricket2}$ and
each sound in the set of test sounds $SS_{t5}$. The results indicate that the
MFD-VL signature has the higher robustness against the pink noise than MFCC13,
and that the MFD-VL signature is not relatively affected by the pink noise. It
is confirmed that the MFD-VL signature is certainly robust against the
diversity cause D2 defined in Subsect. 1.2.

Figure 14: Line charts of the discrimination rates $DR_{mfdvl}$ and
$DR_{mfcc}$ to distinguish the simulated cricket sound $S_{cricket2}$ and each
sound in the set of test sounds $SS_{t5}$.

##  5 Experimental Evaluation (EXP1)

###  5.1 Experimental Setups

To evaluate a descriptiveness of the acoustic feature signatures based on the
multiscale fractal dimension, we developed a similarity search system using
the k-nearest neighbors (k-NN) method. We examined the performances of the
similarity search tasks by using the proposed acoustic feature signatures.

To produce a sound dataset of environmental sounds, sound samples were
collected from the Freesound project, which has stored many types of
environmental sounds uploaded by users worldwide. The Freesound allows users
to share their sounds and describe metadata regarding their shared sounds on
the web. Each sound is labeled with a group of tags that are relatively well
maintained as user generated contents [1]. The tags represent the objects the
users were listening to in their listening experiences. Based on the following
rules defined by the authors, the sounds were chosen and imported to our
dataset to be used in the similarity search system. The sounds that were
tagged with “field-recording” and with lengths between 1 and 600 s were
chosen. This means that these sound clips have been recorded outside a
recording studio. Each sound source outside the studio is typically unique and
never the same. Therefore, sound clips that are tagged with ”field-recording”
have the diversity cause D1. We imported the top 3000 sounds in the descending
order of downloaded number by unspecified users counted by the Freesound’s
system for each sound. Each sound was converted to a uniform format (1
channel, 44100 Hz sampling frequency, 16 bits bit depth, and maximal amplitude
normalized to -0.1 db) for normalization before extracting acoustic features
including EMFD, EMFD-KDE, MFD-VL, and MFCCs. The average length of the
imported sounds is 70.4 s.

###  5.2 Acoustic Feature Extraction

One of the most well-known acoustic features used for ESR is MFCCs. Here, we
used MFCCs for comparing the descriptiveness with our newly-proposed feature
signatures. The SPTK toolkit was used to compute 13-coefficients MFCCs
(MFCC13) and MFCC39, which represents the first and the second-order
derivatives of MFCC13. MFCC13 and MFCC39 were computed using a fixed width
analysis window of length 50 ms. The feature sets of MFCC13 and MFCC39 consist
of mean values of their coefficients of the analysis window. EMFD and EMFD-KDE
consist of the 512 elements defined in Sect. 2 and 3, and MFD-VL consists of
the 10 elements defined in Sect. 4. In this experiment (EXP1), we determined
that the smoothing parameter $\alpha$ used for computing the EMFD-KDE
signature in Eq. (16) is 32.

Table 2 lists different feature sets to be compared through experimental
evaluation. L1 represents the total number of features in the concatenated
feature sets and L2 represents the number of features in the feature sets
after dimensionality reduction through principal components analysis (PCA). To
achieve the best possible performance of the similarity search using k-NN
method, PCA was applied to the feature vectors of the most frequently
downloaded 600 sounds in the dataset to extract its eigenvectors for
dimensionality reduction. The “prcomp” function of R language was used for PCA
processing. The corresponding L2s of feature sets 1, 3, 6, and 8 were
determined so that each of their cumulative contribution ratios was 99%. The
L2s of feature sets 2, 4, 5, 7, 9, and 10 were fixed to 114.

Table 2: List of acoustic feature sets for the comparison of their
descriptiveness.

Feature sets 1 and 6 are the standard sets for comparing with other feature
sets. Each feature set from 2 to 5 consists of MFCC13, and each feature
signature is based on the multiscale fractal dimension. We defined that each
feature set from 1 to 5 belongs to the group FS1. In addition, each feature
set from 7 to 10 consists of MFCC39 and each feature signature is based on the
multiscale fractal dimension. Moreover, we defined each feature set from 6 to
10 to belong to group FS2.

The suffix “($\times~{}\gamma$)” of each feature vector denotes a weighting
coefficient $\gamma$. Each value of the feature vectors is multiplied by
$\gamma$ when its feature vector is combined with other feature(s). Through
the experimental evaluation, the weighting coefficient $\gamma$ for each
feature vector was appropriately chosen to perform the best result.

###  5.3 Evaluation Method

The similarity search system using k-NN method returns a search result list of
environmental sounds based on the distance in the space of the selected
feature set through a search-key sound. To evaluate the performance using each
feature vector in table 2, we defined the similarity index $SI$ between the
tag group of the search-key sound $tags_{key}$ and that of the retrieved sound
$tags_{s}$ as Eq. (36). This index is known as the Jaccard similarity
coefficient that measures similarity between finite sample sets.

$SI=\frac{\mathrm{card}\left(tags_{key}\cap
tags_{s}\right)}{\mathrm{card}\left(tags_{key}\cup tags_{s}\right)}$ (36)

To improve an accuracy of similarity index $SI$, we removed the commoner
morphological and inflexional endings from all tags by using Porter Stemmer
[34] in advance. Furthermore, the predefined stop words include sound formats,
such as “mp3” and “stereo,” and tool makers, such as “sony” and “tascam,” were
removed from the tag groups for computing the $SI$s. The tags that contain the
text “fieldrecord” were removed from the tag groups because all sounds in the
dataset have them. For each of the 3000 sounds in the dataset, $SI$s between a
search-key sound and each retrieved sound in the search-result list were
computed. We compared the average values of the $SI$s of 3000 sounds for each
acoustic feature set.

###  5.4 Evaluation Results

Table 3 shows the $SI$s of “top $n$” for each feature set. The $SI$ of “top
$n$” is the average value of the $SI$s between a search-key sound and each
retrieved sound in the top $n$ rank of the search-result list. Table 3 shows
the average values of the $SI$s of “top $n$” computed for each of the 3000
sounds. For reference, the average value of the $SI$ between two randomly
chosen sounds in the dataset is 0.014.

Table 3: Evaluation results of each feature set quantified using $SI$s

By comparing feature sets 4 and 9 with 2 and 7, respectively, it was confirmed
that the descriptiveness of EMFD-KDE is superior to that of EMFD. $SI$ of “top
1” of feature set 4 was 8.1% higher than that of feature set 2. $SI$ of “top
1” of the feature set 9 was 4.9% higher than that of feature set 7.

By comparing feature sets 3 and 8 with 1 and 6, respectively, it was confirmed
that the newly-developed MFD-VL signature improves the performance of
similarity search result. The MFD-VL signature can describe some acoustic
features that MFCCs cannot. $SI$ of “top 1” of the feature set 3 was 8.2%
higher than that of the feature set 1. $SI$ of “top 1” of feature set 8 was
1.6% higher than that of the feature set 6.

The results obtained using the feature sets 5 and 10 achieved the best
similarity search performance in each group of feature vectors FS1 and FS2.
$SI$ of “top 1” of feature set 5 was 17.2% higher than that of the feature set
1. $SI$ of “top 1” of the feature set 10 was 8.7% higher than that of the
feature set 6. It was confirmed that EMFD-KDE and MFD-VL are effective as an
acoustic feature signature for the environmental sounds.

##  6 Experimental Evaluation Using Public Dataset and Other Acoustic Features
(EXP2)

###  6.1 Experimental Setups

We developed another similarity search system to evaluate the proposed
acoustic feature signatures by comparing them with the top-ranked acoustic
features used in Task2 of the DCASE 2018 challenge [35]. As a sound dataset of
this similarity search system, we used the train set of ”Freesound Dataset
Kaggle 2018” (FSDKaggle2018) [36] which contains 9473 sound clips provided as
uncompressed PCM 16 bit, 44.1 kHz, mono audio files. The significant
characteristics of the dataset are as follows:

  * •

The sound clips are unequally distributed in the 41 categories of Google’s
AudioSet Ontology [37]. The minimum number of sound clips per category is 94,
and the maximum is 300.

  * •

Each sound clip is annotated with a single ground-truth label of the
categories.

  * •

The duration of the sound clips ranges from 300ms to 30s.

The maximal amplitude of sound clips was normalized to -0.1db before the
acoustic feature extraction.

###  6.2 Acoustic Feature Extraction

We chose the top-ranked acoustic features proposed in the DCASE 2018 challenge
task2, including log-mel energies, Perceptual weighted power spectrogram, and
Logarithmic-filtered log-spectrogram, to compare their descriptiveness with
that of the EMFD-KDE and MFD-VL signatures. Table 4 shows the acoustic feature
sets to be used in EXP2. The significant processes of the feature extraction
are as follows:

  * •

All acoustic features except MFD-VL were computed using fixed-width analysis
windows of length 50ms every 25ms.

  * •

For computing the MFD-VL signature, the length of the sound signal must be
longer than 1s. In case the target sound clip is shorter than or equal to 1s,
we repeated the sound signal so that the total length of the repeated signals
becomes longer than 1s.

  * •

Let the smoothing parameter $\alpha$ be 1, which is used for computing the
EMFD-KDE signature in Eq. (16).

  * •

Dimensionality reduction was performed through PCA using scikit-learn library.
The numbers of features L2s after dimensionality reduction were determined so
that each of their cumulative contribution ratios was 98%.

Table 4: List of acoustic feature sets for the comparison of their
descriptiveness.

###  6.3 Evaluation Method

For each of the 9473 sounds in the dataset, the relevances between the search-
key sound and each retrieved sound in the search-result list were computed
using their category labels. Then, we computed the mean value of Precision@k
for each query on the entire dataset to measure the top-$k$ ranking quality of
the search result lists.

###  6.4 Evaluation Results

Table 5 shows the Precision@k scores at ranking positions 1, 3, and 10 on the
entire dataset. Bold numbers are the best Precision@k scores for each ranking
position $k$. The results obtained using the feature set 8’ (MFCC39 + MFD-VL)
achieved the best similarity search performance. By comparing feature sets (2’
with 10’, 11’, and 12’), (3’ with 13’, 14’, and 15’), and (4’ with 16’, 17’,
and 18’) respectively, it was confirmed that both the EMFD-KDE and MFD-VL
signatures improve the performance of the similarity search task. By comparing
feature sets 1’ to 6’, we confirmed that the Precision@10 score obtained using
the feature set 5’ (EMFD-KDE) achieved better performance than those obtained
using the feature sets 2’ (log-mel energies), 3’ (Perceptual weighted power
spectrogram), and 4’ (Logarithmic-filtered log-spectrogram), respectively.

Table 5: Precision@k scores at ranking positions 1,3, and 10 on the entire
dataset.

Table 6 shows the Precision@3 scores of the search-result list grouped by the
41 categories for each feature set. Bold numbers are the best Precision@3
scores for each category. The top five numbers of the categories in which each
feature set obtained the best Precision@3 score are 10 categories for feature
set 8’ (MFCC39 + MFD-VL), 9 categories for feature set 1’ (MFCC39), 7
categories for feature set 9’ (MFCC39 + EMFD-KDE + MFD-VL), 5 categories for
feature set 17’ (Logarithmic-filtered log-spectrogram + MFD-VL), and 4
categories for feature set 18’ (Logarithmic-filtered log-spectrogram + EMFD-
KDE + MFD-VL).

Table 6: Precision@3 scores for each category.

##  7 Conclusion

Recent research on ESR focused on the evaluation of time-domain features of
environmental sounds. For ESR, an acoustic feature must describe the
nonstationary characteristics of target sounds as time-domain features and
must be robust against the following three main causes of the diversity of
environmental sounds.

  * D1)

Small fluctuations of sound source parameters, such as carrier signal
frequency, due to the individuality of the sound source.

  * D2)

Background noises that the person who recorded the target sound did not expect
to record.

  * D3)

Mixed composition of different types of sound sources.

In this study, we have focused on the extraction of acoustic feature
signatures that are robust against the diversities caused by D1 and D2.

In a previous study, we proposed the EMFD feature signature to describe both
frequency- and time-domain features of target sounds.

However, we also recognized the following problems in EMFD.

  * p1)

EMFD includes error values.

  * p2)

It is oversensitive to discriminate the features of environmental sounds.

  * p3)

It lacks the robustness against the diversity of environmental sounds.

  * p4)

It cannot describe the time-domain features for periods longer than 10 ms.

To solve these problems, we proposed the EMFD-KDE and MFD-VL feature
signatures.

The newly-proposed EMFD-KDE feature signature is the probability distribution
of the enhanced MFD values at each radius of the unit disk computed using the
kernel density estimation method. In Subsect. 3.2, we studied the method to
optimize bandwidth $h$ used for kernel density estimation. Based on the normal
reference rule, we defined bandwidth $h_{rbin}(\alpha)$ optimized for each
radius of the unit disk by using Eq. (16). Through the experiments with
different values of the smoothing parameter $\alpha$, we determined that the
best result of the similarity search task EXP1 is obtained for $\alpha=32$ and
that of EXP2 is obtained for $\alpha=1$. We assume that this difference of the
optimized value of $\alpha$ for each task is caused by the diversity
difference of datasets for each task. In the FSDKaggle2018 dataset for EXP2, a
number of the ground truth labels have been manually verified, the number of
categories is only 41, and the length of the sound clip is shorter than 30sec.
The dataset for EXP1 is more diverse and complex than the FSDKaggle2018
dataset.

In Sect. 4, we proposed the MFD-VL signature and demonstrated its
characteristics through experiments using the simulated cricket’s sounds as
follows.

  * •

MFD-VL can discriminate the frequencies of the amplitude envelopes between
22.6 and 1 Hz.

  * •

It can discriminate the shapes of the amplitude envelopes.

  * •

It is robust against the fluctuation of the carrier signal frequency (, i.e.,
the diversity cause D1 defined in Subsect. 1.2).

  * •

It is robust against background noises (, i.e., the diversity cause D2 defined
in Subsect. 1.2).

The MFD-VL signature is expected to describe the time-domain features for
periods longer than 10 ms. The MFD-VL signature shows stability and robustness
against background noises and small fluctuations of the carrier signal’s
frequency.

From the experimental evaluation results (EXP1), we confirmed that the
descriptiveness of EMFD-KDE supplementing MFCC13 and MFCC39 is evidently
higher than that of EMFD supplementing MFCC13 and MFCC39. We conclude that the
smoothness of the EMFD-KDE signature can solve problems p1, p2, and p3.
Furthermore, the experimental evaluation results showed that the MFD-VL
signature supplementing EMFD-KDE and MFCCs improves the performance of the
similarity search. The MFD-VL signature functions as an effective time-domain
feature and can solve problems p3 and p4.

In Sect. 6, we conducted another experimental evaluation (EXP2) using the
FSDKaggle2018 dataset and other acoustic features. We confirmed clear evidence
that both the EMFD-KDE and MFD-VL signatures have the unique descriptiveness
of environmental sound. These signatures are effective when they are used with
other acoustic features, including MFCC39 and the top-ranked acoustic features
in the DCASE 2018 challenge task2.

The EMFD-KDE signature has 512 feature elements, which is relatively more than
those of other feature signatures. This implies that the EMFD-KDE signature
requires more computational time than the conventional methods for feature
extraction. For the similarity search task, feature signatures of library
sounds can be computed a priori which implies that we do not care about the
computation time of the EMFD-KDE signature. While the EMFD-KDE also requires
more searching time than conventional method, this time is negligible compared
to the time required for the retrieval of matched sounds from the library. The
computational time of EMFD-KDE and the length of searching time using EMFD-KDE
do not matter.

Environmental sounds have the acoustic features in their frequency domain, as
well as other important features in the time domain with various time scales.
We conclude that both the EMFD-KDE and MFD-VL signatures can describe the
essential acoustic features of environmental sounds with robustness against
the diversity of environmental sounds. Further studies are needed to evaluate
the performance of other applications using these acoustic feature signatures,
such as classification tasks using machine-learning systems.

## Acknowledgments

This work was supported by JSPS KAKENHI Grant Number 18K11378.

## References

  * [1] V. Akkermans, F. Font, J. Funollet, B. de Jong, G. Roma, S. Togias, and X. Serra, “FREESOUND 2.0: An Improved Platform for Sharing Audio Clips,” 12th Int. Soc. Music Inf. Retr. Conf., (2011).
  * [2] Music Technology Group of Universitat Pompeu Fabra, “The Freesound Project.,” https://www.freesound.org/ (Retrieved 2020-12-20).
  * [3] SoundCloud Limited, “SoundCloud – Hear the world’s sounds,” https://soundcloud.com/ (Retrieved 2020-12-20).
  * [4] F. López, “Environmental sound matter,” (1998).
  * [5] D. Michael, “Toward a Dark Nature Recording,” Organised Sound, vol. 16, no. 03, pp. 206–210 (2011).
  * [6] T. H. Park, J. Lee, J. You, M.-J. Yoo, and J. Turner, “Towards Soundscape Information Retrieval (SIR),” Proc. ICMC—SMC —2014, Athens, Greece, 2014, pp. 1218–1225. (2014).
  * [7] S. Chachada and C. C. J. Kuo, “Environmental sound recognition: A survey,” APSIPA Trans. Signal Inf. Process., vol. 3, (2014).
  * [8] D. Stowell, D. Giannoulis, E. Benetos, M. Lagrange, and M. D. Plumbley, “Detection and Classification of Acoustic Scenes and Events,” IEEE Trans. Multimed., vol. 17, no. 10, pp. 1733–1746, (2015).
  * [9] K. Koutini, F. Henkel, H. Eghbal-zadeh, and G. Widmer, “CP-JKU Submissions to DCASE’20: Low-Complexity Cross-Device Acoustic Scene Classification with RF-Regularized CNNs,” DCASE2020 Challenge, Tech. Rep., (2020).
  * [10] M. Cowling and R. Sitte, “Comparison of techniques for environmental sound recognition,” Pattern Recognit. Lett., vol. 24, no. 15, pp. 2895–2907, (2003).
  * [11] S. Chu, S. Narayanan, and C.-C. Kuo, “Environmental Sound Recognition With Time-Frequency Audio Features,” IEEE Trans. Audio. Speech. Lang. Processing, vol. 17, no. 6, pp. 1142–1158, (2009).
  * [12] R. Mogi and H. Kasai, “Noise-Robust environmental sound classification method based on combination of ICA and MP features,” Artif. Intell. Res., vol. 2, no. 1, p. 107, (2012).
  * [13] C. Bauge, M. Lagrange, J. Anden, and S. Mallat, “Representing environmental sounds using the separable scattering transform,” ICASSP, IEEE Int. Conf. Acoust. Speech Signal Process. - Proc., pp. 8667–8671, (2013).
  * [14] J. Xue, G. Wichern, H. Thornbug, and A. Spanias, “Fast query by example of environmental sounds via robust and efficient cluster-based indexing,” ICASSP, IEEE Int. Conf. Acoust. Speech Signal Process.- Proc., pp. 5–8, (2008).
  * [15] G. Roma, J. Janer, S. Kersten, M. Schirosa, P. Herrera, and X. Serra, “Ecological Acoustics Perspective for Content-Based Retrieval of Environmental Sounds,” EURASIP J. Audio, Speech, Music Process., vol. 2010, pp. 1–11, (2010).
  * [16] G. Chechik, E. Ie, M. Rehn, S. Bengio, and D. Lyon, “Large-scale content-based audio retrieval from text queries,” Vancouver, British Columbia, Canada, pp. 105–112, (2008).
  * [17] M. Sunouchi and Y. Tanaka, “Similarity Search of Freesound Environmental Sound Based on Their Enhanced Multiscale Fractal Dimension,” Sound Music Comput. Conf. 2013, SMC 2013, pp. 715–721, (2013).
  * [18] S. Handel, “Timbre perception and auditory object identification,” Hearing. Academic Press, p. 468, (1995).
  * [19] D. Mitrović, M. Zeppelzauer, and H. Eidenberger, “On Feature Selection in Environmental Sound Recognition,” 51st Int. Symp. ELMAR, no. September, pp. 28–30, (2009).
  * [20] D. Mitrović, M. Zeppelzauer, and C. Breiteneder, “Features for Content-Based Audio Retrieval,” Adv. Comput., 2010, vol. 78, ch. 3, pp. 71–150, (2010).
  * [21] S. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. Signal Process., vol. 41, no. 12, pp. 3397–3415, (1993).
  * [22] S. Innami and H. Kasai, “NMF-based environmental sound source separation using time-variant gain features,” Comput. Math. with Appl., vol. 64, no. 5, pp. 1333–1342, (2012).
  * [23] B. Mandelbrot, “The Fractal Geometry of Nature,” W. H. Freeman and Company, (1982).
  * [24] R. F. Voss and J. Clarke, “ ‘1/f noise’ in music and speech,” Nature, vol. 258, no. 5533, pp. 317–318, (1975).
  * [25] K. J. Hsu and A. J. Hsu, “Fractal Geometry of Music,” Proc. Natl. Acad. Sci., vol. 87, no. 3, pp. 938–941, (1990).
  * [26] P. Maragos, “Fractal aspects of speech signals: dimension and interpolation,” Acoust. Speech, Signal Process. 1991. ICASSP-91. - Proc., pp. 417–420, (1991).
  * [27] P. Maragos and A. Potamianos, “Fractal dimensions of speech sounds: computation and application to automatic speech recognition.” J. Acoust. Soc. Am., vol. 105, no. 3, pp. 1925–1932, (1999).
  * [28] A. Zlatintsi and P. Maragos, “Musical Instruments Signal Analysis and Recognition Using Fractal Features,” 19th Eur. Signal Process. Conf. (EUSIPCO 2011), Barcelona, Spain, no. Eusipco, pp. 684–688, (2011).
  * [29] A. Zlatintsi and P. Maragos, “Multiscale Fractal Analysis of Musical Instrument Signals With Application to Recognition,” IEEE Trans. Audio. Speech. Lang. Processing, vol. 21, no. 4, pp. 737–748, (2013).
  * [30] D. W. Scott, “Multivariate Density Estimation: Theory, Practice, and Visualization,” 1st ed. Wiley, (1992).
  * [31] D. Keller and B. Truax, “Ecologically-based granular synthesis,” Proc. Int. Comput. Music Conf., Ann Arbor, USA, 1998, pp. 117–120, (1998).
  * [32] J. Thorson, T. Weber, and F. Huber, “Auditory behavior of the cricket,” J. Comp. Physiol. A. Neuroethol. Sens. Neural. Behav. Physiol., vol. 146, no. 3, pp. 361–378, (1982).
  * [33] SPTK working group, “Speech Signal Processing Toolkit (SPTK),” http://sp-tk.sourceforge.net/, (Retrieved 2020-12-20).
  * [34] M. Porter, “An algorithm for suffix stripping,” Progr. Electron. Libr. Inf. Syst., vol. 14, no. 3, pp. 130–137, (1980).
  * [35] E. Fonseca et al., “General-purpose tagging of Freesound audio with AudioSet labels: task description, dataset, and baseline,” in Proceedings of the Detection and Classification of Acoustic Scenes and Events 2018 Workshop (DCASE2018), 2018, pp. 69–73.
  * [36] E. Fonseca et al., “FSDKaggle2018,” 29-Jan-2019. [Online]. Available: https://zenodo.org/record/2552860. [Accessed: 23-Apr-2021].
  * [37] J. F. Gemmeke et al., “Audio Set: An ontology and human-labeled dataset for audio events,” 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2017, pp. 776-780, doi: 10.1109/ICASSP.2017.7952261.

[sunou]Motohiro Sunouchireceived masters of environmental studies in human and
engineered environmental studies from the University of Tokyo in 2004. He is
currently pursuing the Ph.D. in information science at the Hokkaido
University. He has been a research assistant since 2007 and a senior lecturer
since 2016 with the Department of Design, Sapporo City University, Japan. His
research interests lie in the areas of audio signal processing and auditory
culture. [yoshioka]Masaharu YoshiokaMasaharu Yoshioka, Ph. D is a Professor of
Faculty of Information Science and Technology, Global Station of Big Data and
Cybersecurity, and Institute for Chemical Reaction Design and Discovery of
Hokkaido University. He received the B.E. and M.E. degrees of precision
engineering and the Ph.D. degree of precision machinery engineering from
University of Tokyo, Japan, in 1991, 1993, and 1996, respectively. From April
1996 to March 2000, he was a Research Associate of National Center for Science
and Information Systems, Japan. From April 2000 to May 2001, he was a Research
Associate of National Institute of Informatics, Japan. From June 2001, he
joined the Graduate School of Engineering as a Associate Professor and this
school is reorganized as Graduate School of Information Science and Technology
in 2004. From January 2019, he became a Professor of Faculty of Information
Science and Technology. He also joined Institute for Chemical Reaction Design
and Discovery from January 2020. From October 2017, he also serve as a
visiting researcher at RIKEN Center for Advanced Intelligence Project. His
research interests includes application of knowledge engineering technology
for information access and knowledge management, Linked Open Data, and
application of knowledge engineering technology for a particular research
domain (e.g., cheminformatics and nanoinformatics).