# The Multiwavelength Picture of GRB 221009A’s Afterglow

Marc Klinger1, Andrew M. Taylor1, Tyler Parsotan2, Andrew Beardmore3,
Sebastian Heinz4 and Sylvia J. Zhu1
1Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen,
Germany
2Astrophysics Science Division, NASA Goddard Space Flight Center,Greenbelt, MD
20771, USA.
3School of Physics and Astronomy, University of Leicester, Leicester LE1 7RH,
UK
4Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter
Street, Madison, WI 53706, USA E-mail<EMAIL_ADDRESS>0000-0002-4697-1465 0000-0001-9473-4758 0000-0002-4299-2517
0000-0002-8433-8652 0000-0002-6468-8292

( Accepted 2023 November 30. Received 2023 November 30; in original form 2023
August 17)

###### Abstract

We present counts-level fits to the multi-instrument (keV–GeV) data of the
early afterglow (4 ks, 22 ks) of the brightest gamma-ray burst detected to
date, GRB 221009A. The complexity of the data reduction, due to the
unprecedented brightness and the location in the Galactic plane, is critically
addressed. The energy spectrum is found to be well described by a smoothly
broken power law with a break energy at a few keV. Three interpretations
(slow/fast cooling or the transition between these) within the framework of
forward shock synchrotron emission, from accelerated and subsequently cooled
electrons, are found. The physical implications for each of these scenarios
are discussed.

###### keywords:

gamma-ray bursts, radiation mechanisms: non-thermal, acceleration of
particles, methods: data analysis

††pubyear: 2023††pagerange: The Multiwavelength Picture of GRB 221009A’s
Afterglow–The Multiwavelength Picture of GRB 221009A’s Afterglow

##  1 Introduction

The non-thermal afterglow emission of gamma-ray bursts (GRBs) constitutes a
two-dimensional flux surface as a function of energy and time. In the energy
dimension the cleanest observational window ranges from hard X-ray (keV) up to
very-high-energy (VHE; >100 GeV) gamma-ray energies. In the VHE band and
above, the absorption due to the extragalactic background light (EBL) limits
spectral inference. From the soft X-ray down to infra-red bands, both
photoelectric absorption and dust scattering also prohibit decisive
conclusions, such that optical and UV data serve only as lower limits to the
intrinsic source flux. We emphasise that the systematic uncertainty introduced
by correcting for these exponential absorption effects is difficult to
constrain. At even lower energies, i.e. from infra-red down to the radio band,
the propagation becomes unhindered once again; however, an additional emission
component, usually attributed to the reverse shock, appears to begin to
dominate within this range, limiting the insight on the forward shock
component. In the dimension of time, the afterglow’s most informative data
window lies at early times, starting with the clear onset of the afterglow.
Towards later times the significance of the GRB signal over the background
quickly drops due to the decaying light curves in combination with the limited
sensitivity of GRB detection instruments.

GRB 221009A is the brightest gamma-ray burst detected to date. Its extensive
multi-wavelength coverage promises a correspondingly large amount of new
insights. Contrary to these expectations, the extraordinarily high flux and
its location in the Galactic plane (latitude b = 4.3 degrees) led to
additional systematic errors in several instruments, hindering the measurement
of the GRB’s flux. This is not only limited to the prompt emission phase
($\lesssim T_{0\mathrm{,GBM}}+700$ s, Lesage et al. 2023); for example its
complex dust echos significantly challenge the determination of the afterglow
X-ray spectrum up to more than 100 ks (e.g. Williams et al., 2023).

We focus thus on the early afterglow with contemporaneous data from keV up to
GeV energies, in particular two time intervals around 4 ks and 22 ks after the
trigger time of the Fermi Gamma-Ray Burst Monitor (Fermi-GBM). In both
intervals the Swift X-ray Telescope (Swift-XRT; Burrows et al. 2005), Swift
Burst Alert Telescope (Swift-BAT; Barthelmy et al. 2005) and Fermi Large Area
Telescope (Fermi-LAT; Atwood et al. 2009) observations overlap temporally.

Several earlier works have already inferred spectral information on GRB
221009A in this energy window in temporal proximity. The Fermi-GBM data shows
an afterglow component with an observed photon spectral index of
$p_{\gamma}=-\differential\log N_{\gamma}/\differential\log E_{\gamma}\approx
2$ (40 keV–8 MeV), already shining through during low phases of the flaring
prompt phase and turning more clearly into the afterglow after $\sim 700$ s
(Lesage et al., 2023; Zhang et al., 2023). At 1.5 ks GRB 221009A gets occulted
by the Earth and at 4 ks only an energy flux level of
$EF_{E}\approx~{}10^{-8}\,{\rm erg}/({\rm cm}^{2}~{}{\rm s})$ at 10 keV is
estimated via the Earth occultation technique (Lesage et al., 2023). Williams
et al. (2023) find preference for a break around 7 keV with photon spectral
indices $p_{\gamma}$ of 1.7 (before the break) and 2.2 (after the break) from
a combined fit of Swift-XRT, Swift-BAT and the heavily absorbed Swift-UVOT
data around 4.2 ks, when linking the spectral indices of a broken power law by
0.5. The residuals show a significant offset between Swift-XRT and Swift-BAT
fluxes, highlighting the challenges of the Swift-XRT data in windowed timing
(WT) mode (see Section S1 of the supplementary materials for details). The
photon spectral index $p_{\gamma}$ above the break seems to be in agreement
with the rough value of $2$ determined from 3 energy bins (20-1220 keV,
5.1-25.7 ks) by Frederiks et al. (2023). Earlier observations (1.35-1.86 ks)
by GECAM-C (20 keV-6 MeV) show a photon spectral index of $p_{\gamma}\approx
2.1$, in agreement with their Earth occultation analysis (20-200 keV),
although suggesting a harder value (1.6) above 200 keV (An et al., 2023). The
combined analysis of Swift-XRT and Fermi-LAT in Laskar et al. (2023) suggests
a break between the two instruments based on the time averaged photon spectral
indices ($p_{\gamma}\approx 1.6-1.9$ for Swift-XRT, $p_{\gamma}\approx 2.1$
for Fermi-LAT), without specific information on the time interval or the data
reduction in WT mode. Other works like Liu et al. (2023) and Stern & Tkachev
(2023) also report photon spectral indices for Fermi-LAT of approximately
$p_{\gamma}\approx 2.1$.

GRB 221009A gained particular attention due to its detection by LHAASO at VHE
energies (LHAASO Collaboration et al., 2023). Significant detection by the
Water Cherenkov Detector Array (WCDA) is however limited to less than 3 ks and
upper limits are not significantly constraining for the time windows we
consider here. A power-law extrapolation of the light curve to 4 ks indicates
a flux at the level of a few $10^{-9}{\rm erg}/({\rm cm}^{2}~{}{\rm s})$ with
big uncertainties of an order of magnitude arising from low statistics and EBL
absorption. Due to poor observational conditions, in particular bright
moonlight, no other VHE telescopes detected GRB 221009A.

Commonly, these observations are interpreted in a one-zone radiation model
coupled to fireball blast wave dynamics (for a review see e.g. Zhang, 2018,
and references therein). Applying this to GRB 221009A required departures from
this simple picture, in various directions (Ren et al., 2023; Sato et al.,
2023; Laskar et al., 2023; Gill & Granot, 2023; O’Connor et al., 2023; Kann et
al., 2023).

In order to clarify this emerging early afterglow spectrum — namely, a hard
spectrum (1.6) up to a break at a few keV followed by a consistent softer
spectral index ($\gtrsim$2) up to 10 GeV — and its compatibility with this
one-zone standard model, we provide in this article a more comprehensive
combined fit on the counts-level, using data from 0.6 keV up to 100 GeV. In
particular, we give a short summary of our reduced afterglow model, the data
reduction and the joint-fit spectral analysis procedure in section 2, present
and discuss our results in sections 3 and 4 and conclude in section 5.

##  2 Methods

###  2.1 Afterglow model

The ultra-relativistic outflow of a GRB jet decelerates as it sweeps up
material from the surrounding environment. A fraction of this swept up
material is shock-accelerated to non-thermal energies and subsequently
radiates away a part of its energy into the X-ray to gamma-ray energy bands.
Our model is described in detail in Klinger et al. (2023) and we give a short
summary of the relevant parts here.

Following simple energy conservation arguments (Blandford & McKee, 1976), we
estimate the self-similarly decelerating bulk Lorentz factor $\Gamma$ of the
outflow to be 34 (18) around 4 ks (22 ks), assuming an isotropic energy of
$E_{\mathrm{iso}}=10^{55}$ erg (Lesage et al., 2023) converted completely into
radiation and a constant circum-burst density
$n_{\text{up}}=1/\mathrm{cm}^{3}$. We emphasise the weak dependence on these
parameters $\Gamma\propto(E_{\mathrm{iso}}/n_{\text{up}})^{1/8}$, which does
not change dramatically for a wind-like density profile. We also assume that a
fraction of the upstream ram pressure is converted to turbulent magnetic
fields downstream (parameterised by $\varepsilon_{B}$, see Klinger et al.
(2023)).

During the particle acceleration process, the two competing energy loss
processes for the electrons, namely adiabatic and radiative energy losses,
collectively dictate the shape of the produced electron spectrum. We
parameterise the latter as a smoothly broken power law with smoothness
parameter $s_{e}$, and spectral indices based on the injected electron
spectral index $p_{e}$ corresponding to the slow/fast cooling regime (see
Klinger et al. (2023) for the definitions of $E_{e}^{\mathrm{min}}$,
$E_{e}^{\mathrm{cb}}$, $E_{e}^{\mathrm{max}}$). We note that we use the
injection energy scale, $E_{e}^{\mathrm{min}}$, as a free parameter, rather
than the often adopted constant fraction parameters (e.g. Sari et al., 1996):
the non-thermal particle number ($\zeta_{e}$) and energy injection
($\varepsilon_{e}$) rates (Although see van der Horst et al., 2014; Misra et
al., 2021; Warren et al., 2015, 2018; Ressler & Laskar, 2017; Asano et al.,
2020). In the results section, however, the corresponding constraint on
$\zeta_{e}$ and $\varepsilon_{e}$ are provided. Furthermore, the parameter
$\eta$ controls the electron acceleration time — in units of the Bohm
acceleration time — which, in turn, dictates the maximum electron and thereby
eventually also synchrotron and inverse Compton photon energy.

The quasi-steady-state energy distribution of these electrons subsequently
dictates the spectral energy distribution (SED) of the non-thermal emission
from the GRB. For the treatment of the radiative emission, we use the reduced
synchrotron-self-Compton (SSC) model of Klinger et al. (2023), within which
both a single-component synchrotron model (reduced syn) as well as a two-
component SSC model (reduced SSC) can be treated as subclasses. It reduces all
the degenerate complexity of the normalisation of the two components of the
photon spectrum to two parameters: 1) the normalisation of the synchrotron
component $F_{\text{syn}}$ at an observed energy of 100 keV and 2) the
relative factor of the inverse Compton component $N_{\text{IC}}$,
phenomenologically defined as the ratio of the maximum value of the energy
flux spectra of both components. This source spectrum is then transformed to
the observer’s frame by assuming a point-like emission zone, an observation
angle of $0^{\circ}$ and a redshift of $z=0.15$.

In summary, we explore with our reduced syn model the parameters $p_{e}$,
$s_{e}$, $\varepsilon_{B}$, $E_{e}^{\mathrm{min}}$, $\eta$, $F_{\text{syn}}$
and in case of the reduced SSC model the additional parameter $N_{\text{IC}}$
(for priors see Section S4 of the supplementary materials). We fix
$E_{\mathrm{iso}}=10^{55}$ erg, $n_{\text{up}}=1~{}\mathrm{cm}^{-3}$ and
$t_{\text{obs}}\approx 4.8~{}(21.8)$ ks.

###  2.2 Data reduction

####  2.2.1 Swift-XRT

The Swift-XRT started observing GRB 221009A 3.3 ks after the Fermi-GBM
trigger. As reported in Williams et al. (2023), the first 12 snapshots of XRT
data, out to 89 ks post-trigger, were obtained in Windowed Timing (WT) mode,
which is a fast (1.78 ms time resolution), 1D, readout mode, used to minimise
the effects of pile-up on the detected photons when a source is bright.

1D spatial profiles from the early WT observations indicated the presence of a
strong dust scattering component, which was later seen as a set of expanding
rings (i.e. dust scattering echos) in the 2D Photon Counting (PC) mode images.
Compared to a regular point source, these rings add a significant, but
intrinsically inseparable component to the 1D region used to extract the GRB
signal for the two time periods covered here ($3.9-4.5$ ks and $21.6-22.1$ ks;
see Figure S1). By modelling the dust echo profile in WT mode, we have
attempted to minimise its effect on the GRB spectral extractions and
identified an appropriate background region for the later time interval (see
Section S1 of the supplementary materials for details).

In addition to the effect of the dust rings, we model the observed GRB flux
suppression due to photoelectric absorption, whose exponential nature
introduces particularly large uncertainties at lower energies ($\lesssim$ few
keV). We use the conventionally adapted combination of an absorber close to
the GRB and a small contribution of the Milky ways absorption
($N_{\mathrm{H}}=5.38\times 10^{21}\,\mathrm{cm}^{-2}$), as in Williams et al.
(2023).

####  2.2.2 Swift-BAT

The Swift-BAT produces survey data when it is not triggering on a GRB. As
outlined in Williams et al. (2023), we used the BatAnalysis python
package111This package is open source and is available on github at:
https://github.com/parsotat/BatAnalysis  (Parsotan et al., 2023) to download
individual Swift-BAT survey data in the time windows 3907-4534s (4 ks) and
21.6-22.1 ks (22 ks) and process them to produce spectra in 8 energy bins from
14-195 keV (compare Section S2 of the supplementary materials).

####  2.2.3 Fermi-LAT

GRB 221009A became visible to Fermi-LAT during orbital windows of around 2 ks,
being occulted by the Earth for a bit more than 3 ks in-between. The exact
overlapping windows with Swift-XRT and Swift-BAT result in very low statistic
for Fermi-LAT, such that we extend the temporal selection range for Fermi-LAT
to the maximal visible windows, 4-6 ks and 21.2-22.5 ks.

Due to the systematic effects from energy dispersion, we select only photons
with energies above 200 MeV (Fermi-LAT Collaboration, 2018a). We limit our
energy range to 100 GeV, whereas the highest detected photon energies are 7.8
GeV (3 GeV) for 4 ks (22 ks).

Motivated by the large point spread function at low energies, we include
photons within a cone of 10° around GRB 221009A. We highlight its distance of
less than 5° from the Galactic plane means that the Galactic diffuse
background emission is non-negligible on these timescales. This limits our
ability to reliably infer the spectral index in the Fermi-LAT energy range. We
find that other background sources play no significant role (see Section S3 of
the supplementary materials for details).

###  2.3 Joint-fit spectral analysis procedure

The joint-fit spectral analysis is based on the method described in Klinger et
al. (2023). We define the likelihood $\mathcal{L}$ for each instrument on the
counts level based on the corresponding suitable statistic: Poisson (4 ks) and
C-stat222Poisson data with Poisson background (Cash, 1979) (22 ks) for Swift-
XRT, $\chi^{2}$ for Swift-BAT, and PG-stat333Poisson data with Gaussian
background for Fermi-LAT. The instrument response function is incorporated via
forward folding. To account for inter-instrumental cross-calibration
uncertainties we add a constant floating norm of up to $\pm 15\%$ to each
instrument (Madsen et al., 2017; Krimm et al., 2013; Fermi-LAT Collaboration,
2018b). To account for the expected large background contribution of the dust
features of 10-20% at 4 ks, we increase the range in this time window to a
factor between 0.95 and 1.3 (see Section S1 of the supplementary materials).

Following a Bayesian approach, we furthermore assume flat priors (see Section
S4 of the supplementary materials) to sample the parameters’ posterior
distributions and integrate the total evidence for model comparison via the
Bayes factor.

##  3 Results

We expect from single-instrument and combined-subset fits (see Section S5 of
the supplementary materials for details) that the SED should roughly resemble
a broken power law. This is based on a statistical preference for a break
around 10 keV at 4 ks, while this is less clear for the 22 ks interval with
lower statistics. At higher energies Swift-BAT and Fermi-LAT can be fitted
together by a power law without statistical preference for a break, with
photon spectral indices of $\approx 2.2$ (4 ks) and $\approx 2.1$ (22 ks). We
note that the large uncovered MeV gap from 0.2-200 MeV and the large
uncertainties of Fermi-LAT’s spectral index complicate the exclusion of
additional features in-between.

We proceed with the results obtained from a combined fit including all three
instruments, in order to extract the spectral properties of this rough broken-
power-law behaviour in the framework of our reduced syn model: The break
energy, the indices below and above the break and a potential additional
spectral feature in between Swift-BAT and Fermi-LAT.

We start by performing a combined fit using our reduced syn model with a
single component ($N_{\mathrm{IC}}=0$), our syn-only case. Motivated by the
non-standard Swift-XRT background treatment, we perturb syn-only case by
releasing our assumption on the Swift-XRT effective area correction (more than
30%), the XRT-floating case (subsection 3.1). Finally, we allow for an
additional SSC component ($N_{\mathrm{IC}}>0$ in the reduced SSC) in the SSC
case (subsection 3.2).

Figure 1 shows for each time bin (4 ks in blue and 22 ks in green) the 3 best-
fit models with a $1\sigma$-envelope and for the syn-only case the
corresponding residuals at the counts-level for all three detectors. The syn-
only case’s posterior distributions are given in supplementary materials
(Section S6). The single instrument power law fits are added to the SED plot
in grey to guide the eye.

Figure 1: Top: $1\sigma$ envelope of intrinsic SED at 4 ks (blue) and 22 ks
(green) for the cases: syn-only ($N_{IC}=0$, darkest), XRT-floating
($N_{IC}=0$ with extended floating norm factor for Swift-XRT, lowest at 10
keV) and SSC ($N_{IC}>0$, resulting in large uncertainties above 10 GeV). Grey
butterflies correspond to single-instrument power law fits (compare Table S5).
The blue line corresponds to the best fit SED and the blue dashed one shows
the observed, absorbed best fit spectrum (dotted line) for 4 ks interval.
Bottom: Counts-level comparison for the syn-only-case (blue dashed line from
top panel) for each detector (Swift-XRT, Swift-BAT and Fermi-LAT) with
corresponding residuals below. Floating norm factors for detectors ("const")
and Fermi-LAT’s galactic diffuse background template ("galdiff") correspond to
default priors from syn-only-case.

The best fit shape of the syn-only case resembles the broken power law
expectation with a break energy of a few keV and a photon spectral index of
$\approx 2.2$.

###  3.1 Perturbation 1: XRT-floating case ($N_{\mathrm{IC}}=0$)

Inspecting the residuals of Swift-XRT and Swift-BAT (Figure 1) as well as the
posterior distributions of their cross-correlation uncertainty factors (Figure
S2) indicates an increased effective area correction between them. The non-
standard background treatment due to the dust ring pollution motivates us to
explore in our XRT-floating case the effect on our results when we extend the
Swift-XRT floating norm prior distribution to a shift of more than 30%
downwards. Such a preference for large effective area corrections of up to a
factor of 2 is not uncommon for bright GRB afterglows (see e.g. Klinger et al.
(2023) for GRB 190114C).

We observe for 4 ks that the Swift-XRT data gets shifted down by a factor of
$2\pm 0.3$, while the correction factors for Swift-BAT and Fermi-LAT remain
negligible ($\approx 1$). This results in a change of the spectral index above
the break from $p_{\gamma}\approx 2.3-2.2$, corresponding to a change of the
electron injection index from $p_{e}\approx 2.6-2.5$. We emphasise that the
uncertainty to this values is likely underestimated due to the systematic
effects in the Swift-XRT analysis as well as the absorption model. For 22 ks,
Swift-XRT requires a shift down by a factor of $2.2\pm 0.5$ and the evidence
for a break is less clear than at 4 ks. Additionally, the posterior
distributions become multimodal and very broad, such that the systematic
uncertainties from the Swift-XRT analysis are limiting the conclusions drawn
from the syn-only case at 22 ks.

We conclude that an extended Swift-XRT floating norm enlarges mainly the
uncertainty on the inferred spectral index of the electrons.

###  3.2 Perturbation 2: SSC case ($N_{\mathrm{IC}}>0$)

An additional spectral component is expected for the case $N_{\mathrm{IC}}>0$,
becoming relevant at the highest energies. The ad-hoc assumption of
acceleration at the Bohm rate ($\eta=1$) results for both time intervals in a
cut-off (the observed synchrotron burn-off) around a few GeV, where the Fermi-
LAT counts become sparse. We investigate with this SSC perturbation the
possibility that an IC component from an additional free parameter
($N_{\mathrm{IC}}>0$) might dominate the Fermi-LAT observation by displacing
the synchrotron cut-off to lower energies.

We find no preference for this behaviour. Instead, the only difference to the
syn-only case’s fit result is the additional uncertainty above a few GeV, see
Figure 1. This visualises the limited power of the Fermi-LAT data to constrain
an IC component.

We conclude that our combined fit is not sensitive to the onset of an
additional IC component.

##  4 Discussion

We discuss now the posterior distribution of the best fit scenarios, resulting
at 4 ks in an apparent small uncertainty on the best fit values of
$\varepsilon_{B}$, $p_{e}$ and $E_{e}^{\mathrm{min}}$. We examine the validity
of these uncertainties in the following. In the reduced syn model, both, the
break energy and the spectral index below ($E_{\gamma}<10~{}$ keV) are driven
by the combination of (i) the photoelectric absorption factor
($N_{\mathrm{H}}$ in proximity to the GRB) and (ii) the intrinsic break
determined by the interplay between magnetic field ($\varepsilon_{B}\to$
cooling break energy) and minimum injected electron energy
($E_{e}^{\mathrm{min}}$). Note that we directly treat $E_{e}^{\mathrm{min}}$
as a fit parameter.

Focusing on the intrinsic break energy within the reduced syn model, the
observed energy $E_{\gamma}^{\mathrm{b,obs}}\approx$ few keV can be matched by
a break in the synchrotron component in three conceptually different
solutions. They correspond to three different features in the electron
spectrum: 1) slow cooling ($E_{e}^{\mathrm{min}}\ll E_{e}^{\mathrm{cb}}$),
i.e., a cooling break from the change from the dominance of adiabatic to
synchrotron cooling; 2) fast cooling ($E_{e}^{\mathrm{min}}\gg
E_{e}^{\mathrm{cb}}$), i.e., a cooling tail in the electron spectrum or 3)
their transition ($E_{e}^{\mathrm{min}}\approx E_{e}^{\mathrm{cb}}$), i.e., a
cut-off in the electron spectrum resulting in the low energy tail of the
synchrotron spectrum of the lowest energy electrons transitioning to the
synchrotron spectrum of the cooled electrons.

Figure 2: Left: posterior distributions in
$E_{e}^{\mathrm{min}}-\varepsilon_{B}$ plane for the 4 ks (small blue region)
and 22 ks (large green boomerang) with the two arms fixing the observed break
energy $E_{\gamma}^{\mathrm{b,obs}}\sim 10$ keV, corresponding to the slow
(red) and fast (purple) cooling solution, meeting in the transition region
(grey). Note that we treat $E_{e}^{\mathrm{min}}$ directly as a parameter
instead of fixing it via $\varepsilon_{e}$ and our prior is
$E_{e}^{\mathrm{min}}>10^{8}$ eV. For three representative solutions S (slow),
F (fast) and T (transition) the differences in electron and observed photon
spectral energy distributions are visualised.

Figure 2 illustrates three representative cases for each regime, in terms of
their electron (bottom) and photon (top) SEDs and shows the syn-only-case’s
posterior distribution of the corresponding parameters ($E_{e}^{\mathrm{min}}$
vs. $\varepsilon_{B}$, on the left). Placing $E_{\gamma}^{\mathrm{b,obs}}\sim$
few keV restricts the parameter space to the two arms (boomerang), as is found
for both time intervals (blue/green for 4 ks/22 ks). Note that we also provide
a second x-axis and y-axis, to indicate the corresponding magnetic field
strength and ratio between number and energy of upstream electrons converted
into non-thermal particles, respectively
444$B=\Gamma\sqrt{32\pi\varepsilon_{B}n_{\text{up}}m_{p}c^{2}}$ and
$\frac{\varepsilon_{e}}{\zeta_{e}}=\frac{(p_{e}-1)}{(p_{e}-2)}\frac{E_{e}^{\mathrm{min}}}{\Gamma
m_{p}c^{2}}\approx\frac{E_{e}^{\mathrm{min}}}{10\;\mathrm{GeV}}$. We note that
in particular the later mapping is based on multiple uncertain assumptions.

For all three solutions the photon spectrum above the break energy corresponds
to the cooled electrons and is the same ($p_{\gamma}=-\differential\log
N_{\gamma}/\differential\log E_{\gamma}=(p_{e}+1)/2$), whereas the spectral
indices below differ: 4/3 for the transition solution, 1.5 for the fast
cooling solution or $(p_{e}+1)/2\approx 1.7$ (using an injected electron
spectral index of $p_{e}=2.4$). Despite the power-law spectral index of
$p_{\gamma}\approx 1.6$ for Swift-XRT alone, we find from our combined fits a
strong preference in the data for a very hard spectrum below a few keV
($p_{\gamma}\approx 0$ from a broken-power law fit), which limits the best fit
parameters on the boomerang further into the lower-right corner, corresponding
to the transition regime.

It is important to note in this context the strong effect (one order of
magnitude at 1 keV) of the exponential absorption factor in this spectral
regime below the break energy, as visible in the dotted line in Figure 1
showing the observed absorbed spectrum.

The interpretation of this hardening to discriminate between the three
solutions is thus limited by the uncertainties of the absorption model,
introduced by e.g. the chemical composition of the photoelectric model (see
e.g. Dalton & Morris 2020; Valan et al. 2023). For heavily absorbed sources
with large statistics like GRB 221009A, it becomes relevant to propagate this
uncertainty into the likelihood evaluation. A quantitative evaluation of these
uncertainties is additionally complicated by the non-standard background
estimation as well as uncertainties in the instrument response from e.g.
instrumental edges and aging effects and goes beyond the scope of this paper.

We thus conclude that an inferred limitation of the solution type to the
transition regime in the 4 ks interval is not robust due to the
underestimation of the above systematic uncertainties.

The inference of physical properties of the radiation zone from the observed
break energy $E_{\gamma}^{\mathrm{b,obs}}$, however, depends strongly on the
solution type (qualitatively similar for wind-like density profiles):

$E_{\gamma}^{\mathrm{b,obs}}\propto\begin{cases}\varepsilon_{B}^{-3/2}E_{\text{iso}}^{-1/2}n_{\text{up}}^{-1}\hskip
42.67912pt&\text{slow cooling}\\\
\varepsilon_{B}^{1/2}E_{\text{iso}}^{1/4}n_{\text{up}}^{-1/4}\left(E_{e}^{\mathrm{min}}\right)^{2},&\text{fast
cooling}\end{cases}$ (1)

Within the slow cooling solution, fixing the observed break position to the
keV regime would effectively fix the comoving magnetic field to $B\approx
0.03\;\mathrm{G}$, $\varepsilon_{B}\to 10^{-5}$ (for $n_{\text{up}}\approx 1$
cm-3). It is important to note that the steady state approximation is only
accurate to a factor of order unity in this case.

On the other hand, the fast cooling solution is compatible with a large range
of magnetic field values, up to a few Gauss. Here the dependence on the
minimum injected electron energy $E_{e}^{\mathrm{min}}$ becomes dominant,
placing it from 10 to 100 GeV. It is important to notice that in the fast
cooling solution, in particular when approaching the transition solution,
$E_{\gamma}^{\mathrm{b,obs}}$ relates directly to the electron spectrum at the
minimum injection scale, where the spectral shape is poorly understood. In
particular the nature of the typically assumed step-like cut-on of the non-
thermal electron spectrum is expected to be affected by the thermal population
(see e.g. Warren et al. 2022), but commonly neglected.

In terms of the micro-physical parameters, the left side of the slow cooling
arm with $\varepsilon_{e}\ll\zeta_{e}$ requires the injection scale of the
non-thermal electrons to be lower than the average proton energy (i.e.
$E_{e}^{\mathrm{min}}<(p-2)/(p-1)\Gamma m_{p}c^{2}\approx 10$ GeV). This is in
contrast to the typical values obtained from particle-in-cell simulations,
$\varepsilon_{e}\approx 3\zeta_{e}$ (Marcowith et al., 2016, ch. 4), a result
which better aligns with the fast cooling solution we find
($\varepsilon_{e}/\zeta_{e}\approx 3-10$). It is also worth mentioning that
for the right side of the figure ($E_{e}^{\mathrm{min}}\sim 100$ GeV), typical
values of $\varepsilon_{e}\sim 0.1$ would require a value of $\zeta\sim 0.01$
(a percent of the electrons). However, effects like pair-loading of the
upstream medium can significantly influence this ratio, preventing us from
drawing strong conclusions on the microphysics (Grošelj et al., 2022).

Comparing our finding of $p_{\gamma}\approx 2.2-2.3$ to Williams et al. (2023)
and the observational picture sketched in section 1, we find consistency.
Modeling results in other works find consistent values of the spectral index
of the electrons, spanning $p_{e}\approx 2.2-2.8$, whereas the estimates for
$\varepsilon_{B}\approx 10^{-1}-10^{-7}$ vary greatly (Ren et al., 2023; Sato
et al., 2023; Laskar et al., 2023; Kann et al., 2023; Gill & Granot, 2023;
LHAASO Collaboration et al., 2023). In particular our results for the electron
spectral index $p_{e}$ help to constrain its value to 2.4–2.5.

Additionally, one could constrain a possible inverse Compton component using
the LHAASO non-observation at 4 ks to $N_{\text{IC}}<1$. We highlight, that
our findings for GRB 221009A are thus also consistent with, but not more
conclusive than the picture emerging from GRB 190114C (Klinger et al., 2023)
and GRB 190829A (H.E.S.S. Collaboration et al., 2021), such that the physical
origin of the VHE afterglow observations can not be determined more clearly
from this brightest GRB detected so far.

##  5 Conclusion

GRB 221009A was extraordinarily bright, which led to complications in the
analysis of the data in several energy bands. Furthermore, its position close
to the Galactic plane led to increased background levels which further
broadened the systematic uncertainties.

Its SED from keV to GeV energies is in both time intervals (4 ks and 22 ks
after $T_{0\mathrm{,GBM}}$) well described by a smoothly broken power-law,
with break energy $E_{\gamma}^{\mathrm{b,obs}}\approx$ few keV. The photon
spectral index above the break is ($p_{\gamma}=-\differential\log
N_{\gamma}/\differential\log E_{\gamma}=2.2-2.3$) with dominant systematic
uncertainties driven by the Swift-XRT background estimation. For the same
reasons the photon spectral index below the break is less clear
($p_{\gamma}\approx 1.5-2$).

Within our GRB spectral model, the inferred break energy
$E_{\gamma}^{\mathrm{b,obs}}\approx$ few keV can have three possible
interpretations. Either (1) as a slow cooling break, resulting in a weak
magnetic field with only an upper limit on the minimum electron energy. Or (2)
a fast cooling break, requiring a stronger magnetic field and a higher minimum
electron energy. Or (3) the transition between (1) and (2), requiring a weak
magnetic field and a large minimum electron energy. We summarise these results
in Table 1. We note that the preference at 4 ks for the transition regime is
not robust, and that both time intervals only limit $B$ and
$E_{e}^{\mathrm{min}}$ to the boomerang shape region shown in Figure 2.

|  |  $B$ [G] | $\varepsilon_{B}$ |  |  $E_{e}^{\mathrm{min}}$ [eV] | $\varepsilon_{e}/\zeta_{e}$
---|---|---|---|---|---|---
1) slow |  | few | $10^{-5}$ |  | $<10^{11}$ | $<30$
2) fast |  | 0.1-few | $10^{-4}-1$ |  | $10^{10}-10^{11}$ | $3-30$
3) transition |  | 0.1 | $10^{-4}$ |  | $10^{11}$ | $30$
Table 1: A summary of the constraints on key parameters for each of the three
solutions of the energy break origin.

## Acknowledgements

MK, AT and SZ acknowledge support from DESY (Zeuthen, Germany), a member of
the Helmholtz Association HGF. The authors would also like to thank Michelle
Tsirou for helpful input on the Fermi-LAT analysis. This work was supported by
the International Helmholtz-Weizmann Research School for Multimessenger
Astronomy, largely funded through the Initiative and Networking Fund of the
Helmholtz Association. SH acknowledges support from NSF grant AST2205917. All
figures by the authors under a CC BY 4.0 license.

## Data Availability

The X-ray data underlying this paper are publicly available from the UK Swift
archive at https://www.swift.ac.uk/swift_live/ or the HEASARC archive at
https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/swift.pl. Fermi-LAT data is
available from https://fermi.gsfc.nasa.gov/ssc/data/access/. The fits have
been performed with the publicly available Multi-Mission Maximum Likelihood
framework (3ML; Vianello et al., 2015) and the UltraNest package (Buchner,
2021). The code of the authors will be made available upon request.

## References

  * An et al. (2023) An Z.-H., et al., 2023, arXiv e-prints, p. arXiv:2303.01203
  * Asano et al. (2020) Asano K., Murase K., Toma K., 2020, ApJ, 905, 105
  * Atwood et al. (2009) Atwood W. B., et al., 2009, ApJ, 697, 1071
  * Barthelmy et al. (2005) Barthelmy S. D., et al., 2005, Space Sci. Rev., 120, 143
  * Blandford & McKee (1976) Blandford R. D., McKee C. F., 1976, Physics of Fluids, 19, 1130
  * Buchner (2021) Buchner J., 2021, J. Open Source Softw., 6, 3001
  * Burrows et al. (2005) Burrows D. N., et al., 2005, Space Sci. Rev., 120, 165
  * Cash (1979) Cash W., 1979, ApJ, 228, 939
  * Dalton & Morris (2020) Dalton T., Morris S. L., 2020, MNRAS, 495, 2342
  * Fermi-LAT Collaboration (2018b) Fermi-LAT Collaboration 2018b, Caveats About Analyzing LAT Pass 8 (P8R3_V3) Data, https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html, Accessed: 2023-11-30
  * Fermi-LAT Collaboration (2018a) Fermi-LAT Collaboration 2018a, Pass 8 Analysis and Energy Dispersion, https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Pass8_edisp_usage.html, Accessed: 2023-11-30
  * Frederiks et al. (2023) Frederiks D., et al., 2023, ApJ, 949, L7
  * Gill & Granot (2023) Gill R., Granot J., 2023, MNRAS, 524, L78
  * Grošelj et al. (2022) Grošelj D., Sironi L., Beloborodov A. M., 2022, ApJ, 933, 74
  * H.E.S.S. Collaboration et al. (2021) H.E.S.S. Collaboration Abdalla H., Aharonian F., Ait Benkhali F., Angüner E. O., Arcaro C., et al., 2021, Science, 372, 1081
  * Kann et al. (2023) Kann D. A., et al., 2023, ApJ, 948, L12
  * Klinger et al. (2023) Klinger M., Tak D., Taylor A. M., Zhu S. J., 2023, MNRAS, 520, 839
  * Krimm et al. (2013) Krimm H. A., et al., 2013, ApJS, 209, 14
  * LHAASO Collaboration et al. (2023) LHAASO Collaboration et al., 2023, Science, 0, eadg9328
  * Laskar et al. (2023) Laskar T., et al., 2023, ApJ, 946, L23
  * Lesage et al. (2023) Lesage S., et al., 2023, ApJ, 952, L42
  * Liu et al. (2023) Liu R.-Y., Zhang H.-M., Wang X.-Y., 2023, ApJ, 943, L2
  * Madsen et al. (2017) Madsen K. K., Beardmore A. P., Forster K., Guainazzi M., Marshall H. L., Miller E. D., Page K. L., Stuhlinger M., 2017, AJ, 153, 2
  * Marcowith et al. (2016) Marcowith A., et al., 2016, Reports on Progress in Physics, 79, 046901
  * Misra et al. (2021) Misra K., et al., 2021, MNRAS, 504, 5685
  * O’Connor et al. (2023) O’Connor B., et al., 2023, Science Advances, 9, eadi1405
  * Parsotan et al. (2023) Parsotan T., Laha S., Palmer D. M., Lien A., Cenko S. B., Krimm H., Markwardt C., 2023, ApJ, 953, 155
  * Ren et al. (2023) Ren J., Wang Y., Zhang L.-L., Dai Z.-G., 2023, ApJ, 947, 53
  * Ressler & Laskar (2017) Ressler S. M., Laskar T., 2017, ApJ, 845, 150
  * Sari et al. (1996) Sari R., Narayan R., Piran T., 1996, ApJ, 473, 204
  * Sato et al. (2023) Sato Y., Murase K., Ohira Y., Yamazaki R., 2023, MNRAS, 522, L56
  * Stern & Tkachev (2023) Stern B., Tkachev I., 2023, Soviet Journal of Experimental and Theoretical Physics Letters,
  * Valan et al. (2023) Valan V., Larsson J., Ahlgren B., 2023, ApJ, 944, 73
  * Vianello et al. (2015) Vianello G., et al., 2015, arXiv e-prints, p. arXiv:1507.08343
  * Warren et al. (2015) Warren D. C., Ellison D. C., Bykov A. M., Lee S.-H., 2015, MNRAS, 452, 431
  * Warren et al. (2018) Warren D. C., Barkov M. V., Ito H., Nagataki S., Laskar T., 2018, MNRAS, 480, 4060
  * Warren et al. (2022) Warren D. C., Dainotti M., Barkov M. V., Ahlgren B., Ito H., Nagataki S., 2022, ApJ, 924, 40
  * Williams et al. (2023) Williams M. A., et al., 2023, ApJ, 946, L24
  * Zhang (2018) Zhang B., 2018, The Physics of Gamma-Ray Bursts, doi:10.1017/9781139226530.
  * Zhang et al. (2023) Zhang H.-M., Huang Y.-Y., Liu R.-Y., Wang X.-Y., 2023, ApJ, 956, L21
  * van der Horst et al. (2014) van der Horst A. J., et al., 2014, MNRAS, 444, 3151


# Distortion measure of spectrograms for
classification of respiratory diseases

Jeremy Levy, Alexander Naitsat and Yehoshua Y. Zeevi Jeremy Levy is with the
Faculty of Electrical Engineering and the Faculty of Biomedical Engineering,
Technion—Israel Institute of Technology, Haifa 3200003 IsraelAlexander Naitsat
and Yehoshua Zeevi are with the Faculty of Electrical Engineering,
Technion—Israel Institute of Technology, Haifa 3200003, Israel

###### Abstract

A new method for classification of respiratory diseases is presented. The
method is based on a novel class of features, extracted from pulmonary sounds,
by parameterizing their spectrograms that are represented as surfaces, and by
utilizing geometrical distortions defined with reference to these surfaces.
This method yields a set of highly descriptive features for analysis of
pulmonary sound recordings. Furthermore, by combining these features with Mel-
frequency cepstral coefficients, we introduce a powerful model for automatic
diagnosis of common respiratory pathologies. Compared with baseline methods,
our model achieves superior results for binary and multi-class classifications
of common respiratory deceases. Our new approach to classification of one-
dimensional signals is applicable to other signals in the context of their
representations in combined spaces or manifolds.

##  I Introduction

Machine learning has begun to have impact on various facets of medical signals
and image analysis [1, 2]. Computational means for lung sounds analysis have
been the subject of various studies using machine learning [3] and other
related approaches [4], aiming at detecting respiratory pathologies.

Nevertheless, the shortage in well-trained physicians who specialize in
osculation pulmonary diseases is alarming. This is especially the case in
certain parts of the word like Africa, where the spread of obstructive lung
diseases is already beyond control [5]. Since advanced means of mass screening
by medical imaging is not an option in such places, let alone that expert
radiologists are not available there, there is an urgent need for inexpensive,
yet reliable methods for automatic diagnostics of lung diseases by means of an
electronic stethoscope.

In view of the above outstanding medical problem, it is of utmost importance
to develop new algorithms for classification of lung sounds and to advance a
fully automatic pipeline for diagnosis of common respiratory pathologies.

The physics of sounds and the natural characteristic of nonstationarity of
lung sounds, lend themselves to their analysis in the combined space of time-
frequency, i.e. the spectrogram. We carry this approach one step further, by
considering the spectrogram as a surface and by exploiting its geometric
properties in generating novel features that are powerful in classification of
lung diseases by means of machine learning.

##  II Related work

There exist several deep learning and model-based methods for automatic
classification of lung pathologies based on their fingerprints that are hidden
in the pulmonary sounds. For instance, the recent method of [6] implements a
deep transfer learning-based multi-class classifier for diagnosis of COVID–19,
using cough recordings.

Chanbres et al. [7] employ the algorithm of the Essentia library [8] for
extracting sound features from cough recordings. This system was trained on
the dataset of the ICBHI 2017 challenge [9] and it uses a boosted-decisional-
tree algorithm to classify sounds like crackles and wheezes.

Kandaswamy et al. [10] developed a method of analyzing lung sounds signals
using wavelet transform, combined with Artificial Neural Network (ANN) for the
classification. They obtained an accuracy of 94.56% on the validation set and
91.33% on the test set. The dataset was composed of 126 recordings,
categorized into inspiratory wheezes, fine crackles, stridor, squawk, and
rhonchus, apart from normal vesicular sounds.

A similar method with ANN layer was used in [11] for the analysis of cough
sounds recorded from pediatric patients.

Sankur et al. [12] proposed an auto-regressive (AR) model combined with KNN
classifier for the task of diagnosis of the lung sounds. They obtained Correct
Classification (CC) of 93.75% on the test set. The dataset was composed of
asthma, chronic bronchitis, emphysema, pneumonia, pleurisy and bronchiectasis
pathologies.

##  III Methods

###  III-A Database analysis

The Respiratory Sound Database from
(https://www.kaggle.com/vbookshelf/respiratory-sound-database) was used for
the analysis. A total of 918 lung sounds recordings from 126 patients were
used. This database incorporates 7 different pathologies: URTI, Asthma, COPD,
LRTI, Bronchectasis, Pneumonia, Bronchiolitis, and healthy recordings.

The histogram depicted in Figure 1 presents the repartition of the labels
among the cases included in the database. Due to the very poor occurrence of
the Asthma and LRTI pathologies, the corresponding recordings were excluded.

Figure 1: Repartition of labels among the cases included in this database. The
database is imbalanced, about 50% of the patients have COPD disease.

###  III-B Preprocessing

One of the major problems that one must overcome in the process of analysis of
lung sounds is $S/N$ level. Sound generated by instruments and other ambient
activities affect significantly the quality of the lung sounds signal. It is
therefore crucial to improve the level of the $S/N$ without distorting the
stethoscope’s signal.

Our algorithm employs the Saviztky-Golay filter [13] for denoising lung
sounds. The purpose of this filter is to smooth the signal, and increase the
SNR without altering the signal. This filter has been widely used in the field
of time series analysis [14], especially for lung sound analysis [15].

The filter aims to fit a specific polynomial suitable for a signal frame,
using least squares method. The central point of the window is replaced with
that of the polynomial, producing a smoother signal. Denote a polynomial of
the degree $N$ by

$p(n)=\,\sum\limits_{k=1}^{N}\,a_{k}n^{k}\,,$

then, the aim of Saviztky-Golay filter is to minimize the following error:

${\mathcal{E}}_{N}=\,\sum\limits_{n=-M}^{M}\,\Bigl{(}p(n)-x[n])^{2}\Bigr{)}\,,$

where $2M+1$ is the width of the window.

A large value of $M$ will produce a smoother signal, but may neglect some
important variations in the signal. A low value of $M$ may ”over fit” the
data. This is in a way the principle of uncertainty. Secondly, $N$, which
specify the degree of the polynomial can produce a smooth signal for low
value. On the other hand, high value of $N$ may “over fit” the data. By
experimenting with various combinations of these filter parameters, we
converged on the values of $N=3,M=11$ that yielded the best results. Figure 2
shows an example of a filtered signal, and the corresponding raw data.

Figure 2: Effect of the Saviztky-Golay filter on lung sounds recording. Over
the segment of the time period 530–580 sec., the signal is very noisy. The
Saviztky-Golay filter eliminates the sharp changes.

###  III-C Feature Engineering

Our algorithm for the feature extraction can be decomposed into the two main
parts: the signal processing and the geometry processing parts.

In the signal processing part, we preprocess data and compute spectrogram
images of the denoised audio signals.

In the geometry processing part, the original problem of signal processing is
converted to another problem in which we analyze geometric representations of
time-varying signals. The first goal of the geometry processing part is to
represent spectrogram images by discrete surfaces, and the second goal is to
reduce the problem of signal processing to the problem of analyzing shapes of
the obtained surfaces. In this process, each sound signal is, first,
identified with its spectrogam surface, and then the measure of dissimilarity
between two signals is estimated as an effort that takes to deform one
spectrogam surface onto the another.

By using various metrics for estimating this effort, we are capable of
computing different geometrical quantities associated with spectrogram
surfaces. We refer to these quantities as to the geometrical distortions or
distortion measures.

In our method, geometrical distortions are used as descriptive features for
analyzing the original sound recordings. Based on these features, we train an
automatic system for classification of respiratory diseases from sound
recording.

Summarizing all the above, our algorithm for extracting features consists of
the following steps:

  1. 1.

First, we use the FFT algorithm to compute the spectrogram of each lung sounds
recording. Each recorded signal is represented as a spectrogram surface
$z(x,y)$, where $x$ is the time in seconds, $y$ is the normalized frequency
and $z$ is the corresponding power-frequency (see Figure 3.

  2. 2.

We next triangulate the spectrogram surfaces $z(x,y)$ using the Delaunay
triangulation algorithm. That is, we take $n$ samples
$\operatorname{\mathcal{V}}=\big{\\{}(x_{1},y_{1},z_{1}),\dots,(x_{n},y_{n},z_{n})\big{\\}}$
of each spectrogram surface and compute a set of triangles
$\operatorname{\Im}$ that connect vertices over $\operatorname{\mathcal{V}}$.
The mesh of triangles obtained by that triangulation is denoted by
$(\operatorname{\mathcal{V}},\operatorname{\Im})$. It provides a discrete
representation of the original sound signal (see supplemental material for
details).

  3. 3.

Each mesh $M=(\operatorname{\mathcal{V}},\operatorname{\Im})$ is then mapped
to the plane by a piecewise affine map (simplicial map). Here by piecewise
affine or simplicial map we refer to a map $f$ of a mesh $M$, defined in such
a way that: (i) the restriction of $f$ to each triangle
$t\in\operatorname{\Im}$ is an affine transformation, denoted by $f_{t}$; (ii)
components of $f$ coincide on common edges of neighboring triangles. [See
examples of such mappings in Figure 4.]

  4. 4.

Finally, we compute geometrical distortions induced by each simplicial map $f$
and use these distortions as descriptive features of the original signal (see
examples of color encoded distortions in Figure 4).

Apparently, the first two stages are aimed at computing a discrete geometric
representation of audio signals, whereas the goal of the last two stages is to
compare these discrete representations. The details of the above algorithm
stages are presented in the next sections.

###  III-D Sampling and Triangulation

Assume that $S$ is a manifold surface, embedded in $\mathbb{R}^{3}$, and that
$\operatorname{\mathcal{V}}$ is a finite set of points (vertices) sampled on
$S$. Then, a common way to discretize $S$ is to divide that surface into a
finite set of triangles $\operatorname{\Im}$ such that: (i) vertices of the
obtained triangles belong to $\operatorname{\mathcal{V}}$; $(\textbf{ii})$ for
any pair of non-disjoint triangles $t_{1},t_{2}\in\operatorname{\Im}$ the
intersection $t_{1}\cap t_{2}$ is either a common edge of $t_{1}$ and $t_{2}$
or a common vertex of these triangles. We will refer to the pair
$(\operatorname{\mathcal{V}},\operatorname{\Im})$ as to the triangle mesh of a
surface $S$.

In our case, each input signal $I$ is represented by a spectrogram surface
$S=S(I)$ that can be written in following parametric form:

$S(I)=\big{\\{}\big{(}x,y,z(x,y)\big{)}|x\in X,y\in Y\big{\\}}\,,$ (1)

where $X$ and $Y$ are the time and frequency ranges of the audio signal $I$.
We divide $X$ and $Y$ into a number of uniformly distributed points
$x_{1},...,x_{N}$ and $y_{1},...,y_{N}$, and the vertex set
$\operatorname{\mathcal{V}}$ of $S$ is defined by

$\operatorname{\mathcal{V}}=\big{\\{}(x_{i},y_{j},z(x_{i},y_{j}))|\,i,j=1,2,\dots,N\big{\\}}\,.$

We tested the above sampling method and found empirically that it works well
for $N=150$. However, in some scenarios, using adaptive sampling can
potentially yield even better results. See Section VI for the discussion on
more advanced sampling schemes.

The triangle set $\operatorname{\Im}$ of $S$ is constructed by the standard
algorithm for Delaunay triangulation [16] that minimizes the minimum angle of
all the angles of triangles in $\operatorname{\Im}$. This triangulation
algorithm avoids generation of slim triangles which appearance may lead to
numerical issues in the stage of the feature extraction.

After representing data by triangular meshes, we proceed to the next step of
analyzing geometric properties of these meshes.

###  III-E Shape descriptors

Given two meshes of spectrogram surfaces, we wish to define a metric suitable
for of quantifying geometrical (de)similarities between these meshes. In
computer vision, such metrics are often referred to as shape descriptors.
There exist many approaches to computing shape descriptors for a collection of
3D objects. These approaches can be divided qualitatively into the following
categories:

  1. 1.

Spectral methods. In spectral approaches, shape descriptors are derived from
discrete representations of the Laplace-Beltrami operator, defined on surfaces
[17]. Cotangent weights are most commonly used for approximating Laplace-
Beltrami operators over meshes. By using cotangent weights, the Laplace
operator action on a mesh $M$ can be represented by a sparse Laplacian matrix
$L=L(M)$. In such a case, the spectral descriptors of
$M=(\operatorname{\mathcal{V}},\operatorname{\Im})$ are often defined as
$n$-largest eigenvalues of $L$, for a constant number
$n<|\operatorname{\mathcal{V}}|$ [18].

  2. 2.

Metric methods. These methods represent each mesh $M$ by a matrix $G$ of
pairwise distances between vertices of $M$. Usually, these are the Euclidean
or geodesic distances. A desimilarity measure between two meshes $M_{1}$ and
$M_{2}$ is defined in the metric approaches as a function of the distance
matrices $G_{1}$ and $G_{2}$ of these two meshes. For example, metric
descriptors of triangulated surfaces can be obtained by solving the problem of
the General Multi-Dimensional Scaling (GMDS) [19], or by solving other related
problems that involve computations of geodesic distances [20, 21].

  3. 3.

Deformation-based methods. In deformation-based methods, a distance between
two shapes $S_{1}$ and $S_{2}$ is estimated by computing an optimal
deformation $f_{12}$ of $S_{1}$ onto $S_{2}$ and by measuring changes in
various geometric features induced by $f_{12}$. There exist many criteria for
definition of map’s optimality. Most of these criteria are targeted at
preserving the map injectivity and avoiding visual distortions, as much as
possible.

Note that for a large collection $\\{S_{1},...,S_{m}\\}$ of shapes it may be
very expensive to compute optimal deformations $f_{ij}$, for each $1\leq
i<j\leq m$. Therefore, instead of matching all the pairs of shapes, a more
practical approach would be to compute an optimal mapping $f_{i}$ of each
shape $S_{i}$ into a simple target domain. Depending on the dimensionality of
an object and its topology, a simple target domain could be the plane, a
sphere [22], the unit circle [23], a cube [24], or a ball [25].

Our model employs deformation-based descriptors for measuring similarities
between triangle meshes. Note that all the meshes that constitute a peak
surface of spectrograms have the topology of a planar disc. Therefore, a
natural candidate for the optimal deformation of a such mesh $M$ is a length-
minimizing mapping of $M$ into the plane. We refer to this mapping process as
to the surface flattening, for short. In our model, surface flattening
algorithms are used for computing deformation-based descriptors of
spectrographic shapes.

If $f$ is a flattening of a mesh $M$, then we select the shape descriptors of
$M$ to be the geometrical distortions that measure how Euclidean lengths are
deformed under $f$. In such a way, each mesh $M$ can be associated with its
signature vector $(E_{1},...E_{2})$, where numbers $E_{i}$ are various
estimates of the metric deviations induced by flattening $M$ into the plane.

In the the sequel, we address in the detail the surface flattening and the
distortion estimation processes.

###  III-F Surface flattening

The problem of flattening triangular meshes into the plane, also referred to
as the parametrization problem, constitutes one of the central issues in
geometry processing. Consequently, there exist many algorithms for flattening
triangulated surfaces111See [26, 27] for the survey of parametrization
algorithms and for other related methods used in geometry processing.. These
algorithms are aimed at computing a locally-injective parametrization that
minimizes distortions of fundamental geometric quantities, such as angles and
lengths. Surface parametrization tasks can be reduced to the following
optimization problem:

$\displaystyle f^{*}={}$
$\displaystyle\underset{f}{\operatorname{argmin}}~{}E(f);$ (2) s.t.
$\displaystyle\det df_{t}>0,t\in\Im\,,$

where $f^{*}$ is a piecewise affine mapping of a mesh
$(\operatorname{\mathcal{V}},\operatorname{\Im})$ that minimizes the chosen
distortion criteria $E$ under the following constraints: for each mesh
triangle $t$, the component of $f^{*}$ on $t$ is an orientation preserving
map. These constraints are expressed by the determinant signs of Jacobian
matrices $df_{t}$, $t\in\operatorname{\Im}$. Negative determinants of the
Jacobians yield inverted triangles in the image of $f$. Therefore, satisfying
the orientation constraints is the necessary condition for inducing one-to-one
parametrization of surface meshes.

We adopt the recently-proposed Adaptive Block Coordinate Descent (ABCD)
algorithm [28], combined with the Tutte embedding method [29], to solve the
optimization problem (2) and thereby the parametrization problem. In
particular, we initialize the parametrization problem (2) by mapping
triangular meshes onto a circle via the method of [29]. We then employ the
ABCD algorithm to induce locally injective parametrization characterized by
minimal length distortions.

Note that, since (2) is a non-convex problem, solving it with different
initial maps may lead to distinct local minima. Therefore, choosing an
appropriate initialization method is crucial for adequate approximation of the
global minimizer $f^{*}$.

We tested a number of different initialization schemes, and found that using a
convex combination mapping of meshes [29] onto a planar disc, yields the best
results. Note that the algorithm of [29] is actually a variant of the
classical Tutte embedding algorithm that is widely used in shape processing
applications. This method guarantees a bijective mapping onto convex planar
domains and it has a low computational cost. Figure 3 demonstrates this
initialization scheme and the related process of distortion minimization.

Figure 3: Visualization of the process of measuring distortion measures for a
given lung sound segment. The figure shows from left to right: a sound
spectrogram, a triangular mesh
$(\operatorname{\mathcal{V}},\operatorname{\Im})$ representation of the
spectrogram surface, an initial mapping of
$(\operatorname{\mathcal{V}},\operatorname{\Im})$ onto the plane and the final
flattening of the mesh, computed by means of the ABCD algorithm.

We proceed to discuss the process of feature extraction. It includes the local
sub-step of extracting features of individual triangles and the global sub-
step in which local features are summed over large subsets of mesh triangles.

###  III-G Measuring local distortions

If $M=(\operatorname{\mathcal{V}},\operatorname{\Im})$ is a triangle mesh and
$f$ is a simplicial mapping of $M$, then a local distortion induced by $f$, on
a triangle $t$, is defined to be a function $E(\sigma_{1},\sigma_{2})$ of the
singular values $\sigma_{1}(df_{t})$ and $\sigma_{2}(df_{t})$ of the Jacobian
$df_{t}$.

The Jacobian singular values uniquely define the shape of a triangulated
surface, up to rotation and sliding of mesh triangles. Generally speaking,
local distortions estimates how extensively is the shape of $t$ distorted
under $f$.

These measures are instrumental in many applications in computer vision,
including shape classification and shape analysis [25, 30]. In our algorithm,
geometric distortions are used as measures of dissimilarity of triangulated
surfaces222In the context of the lung sound representation on the surface of
the local spectrum, the geometric distortions assess how much the local
spectrum of the sound is affected by the simplicial mapping $f$..

Note that for a dense triangulation, feeding singular values
$\big{\\{}\sigma_{i}(df_{t})|t\in\operatorname{\Im},\,i=1,2\big{\\}}$ to a
deep learning model preserves all the information contained in the pixels of
the spectrogram.

Our algorithm employs several distortion measures. These distortions belong to
the following major classes of geometric measures:

Isometric distortions. These measures estimate distortions of the Euclidean
length. We use the following isometric distortions:

  * •

ARAP (As-Rigid-As-Possible) energy [31]

$E_{\mathrm{ARAP}}(\sigma_{1},\sigma_{2})=(\sigma_{1}^{2}-1)^{2}+(\sigma_{2}^{2}-1)^{2}\,;$

  * •

Symmetric Dirichlet energy [32]

$E_{\mathrm{SD}}(\sigma_{1},\sigma_{2})=\displaystyle\frac{1}{4}\,(\sigma_{1}^{2}+\sigma_{1}^{-2}+\sigma_{2}^{2}+\sigma_{2}^{-2})\,;$
(3)

  * •

Quasi-isometric $(qi)$ dilatation [33],[34]

$E_{\mathrm{QI}}(\sigma_{1},\sigma_{2})=\max\limits\,\\{\sigma_{1},\sigma_{2}^{-1}\\}\,;$

Conformal distortions. These distortions estimate how far $f$ is from being an
angle-preserving mapping. Our algorithm uses the following estimates of
conformal distortions:

  * •

Quasi-conformal (qc) dilatation [35]

$E_{\mathrm{QC}}(\sigma_{1},\sigma_{2})=\max\limits\,\left\\{\displaystyle\frac{\sigma_{1}}{\sigma_{2}},\frac{\sigma_{2}}{\sigma_{1}}\right\\}\,;$
(4)

  * •

MIPS energy [36],[37]

$E_{\mathrm{MIPS}}(\sigma_{1},\sigma_{2})=\displaystyle\frac{\sigma_{1}}{\sigma_{2}}+\displaystyle\frac{\sigma_{2}}{\sigma_{1}}=\displaystyle\frac{\sigma_{1}^{2}+\sigma_{2}^{2}}{\sigma_{1}\sigma_{2}}\,;$
(5)

MIPS (Most Isometric Parametrizations) is a quadratic function, widely used
for optimizing conformal distortions over triangular domains [37].

Area distortions. These distortions estimate dilatation and compression of
triangle areas induced by $f$. We use the following measure of the area
distortion:

  * •

Unsigned area distortion [38]

$E_{\mathit{AD}}(\sigma_{1},\sigma_{2})=\max\limits\,\Bigl{\\{}|\sigma_{1}\sigma_{2}|,|\,\sigma_{1}\sigma_{2}|^{-1}\Bigr{\\}}\,;$
(6)

Scale distortions. These distortions assess the degree to which mesh triangles
are scaled by $f$. Scale distortions are closely related to discrete harmonic
mappings [39] and to stretch-minimization mappings [40]. We use the following
scale distortions:

  * •

Dirichlet energy [36]

$E_{\text{Dirichlet}}\,(\sigma_{1},\sigma_{2})=\frac{1}{2}\,\Bigl{(}\sigma_{1}^{2}+\sigma_{2}^{2}\Bigr{)}\,;$
(7)

  * •

Conformal factor [30]

$E_{\text{CF}}(\sigma_{1},\sigma_{2})=\,\displaystyle\frac{\sigma_{1}+\sigma_{2}}{2}\,.$
(8)

Note that conformal factors are closely related to conformal distortions such
as quasi-conformal dilatation and MIPS energy. Indeed, according to the
uniformization theorem [41] any disc-topology surface $S$ can be mapped into
the plane by a conformal map $f_{S}$. The map $f_{S}$ can be described by its
conformal factors, up to a composition of $f_{S}$ with a rigid transformation.
For this reason, the conformal factor has been used by [30] as a geometric
signature for a collection of 3D surfaces.

All those distortion measures are rotation invariant, since they are functions
of signed singular values of the Jacobian. This work aims to show that the
dimensionality of the data can be considerably reduced by employing weighted
sums of local distortions over different subsets of $\operatorname{\Im}$. The
obtained quantities will be referred to as a global distortions.

###  III-H Measuring global distortions

Let $f$ be a simplicial map of the mesh
$(\operatorname{\mathcal{V}},\operatorname{\Im})$, $E$ be a local distortion,
$\operatorname{\Im}_{0}$ be a subset of $\operatorname{\Im}$. The global
distortion of $f$, computed with respect to $E$ over $\operatorname{\Im}_{0}$,
is then defined as follows:

$D_{\operatorname{\Im}_{0}}(f,E)=\displaystyle\frac{\sum\limits_{t\in\operatorname{\Im}_{0}}\,E\big{(}\sigma_{1}(df_{t}),\sigma_{2}(df_{t})\big{)}\text{area}(t)}{\sum\limits_{t\in\operatorname{\Im}_{0}}\,\text{area}(t)}\,,$
(9)

where $df_{t}$ is the Jacobian of $f$ on $t$, $\sigma_{1}(df_{t})$ and
$\sigma_{2}(df_{t})$ are the Jacobian singular values and $\text{area}(t)$
denotes the area of a triangle $t$.

Figure 4: Flattening surfaces and measuring the resultant geometrical
distortions, attained as shape descriptors. The process is visualized for the
MIPS distortion energy, defined by (5). We first use ABCD algorithm [28],
initialized with Tutte embedding, to map the triangulated surfaces into the
plane. We then compute Jacobian singular values $\sigma_{1}(t)$ and
$\sigma_{2}(t)$ over each mesh triangle $t$. Finally, for each distortion
measure $E(\sigma_{1},\sigma_{2})$ and mapping $f$ we compute the two
quantities $E_{1}(f)$ and $E_{2}(f)$, defined according to (9) and (10).
Figure 5: Visualization of a simplicial map $f$ with color encoded local
distortions $E_{\text{SD}}$ , defined by (3). At the top, we show how an
individual triangle $t$ is mapped under an affine component $f_{t}$ of the map
$f$.

In many cases, values of local distortions are distributed non-uniformly over
mesh triangles. As demonstrated by Figures 4 and 5, a small number of highly
distorted triangles may have more impact on the global distortion
$D_{\operatorname{\Im}}(f,E)$ than the rest of the mesh triangles. Therefore,
in order to extract more information from each distortion measure, one can
divide the triangle set $\operatorname{\Im}$ into a number of disjoint
subsets. We employ this approach to extract more features for each distortion
measure $E(\sigma_{1},\sigma_{2})$ and to compensate for the adverse effects
of a non-uniform distribution of distortions. In particular, we divide
triangles into the two subsets according to triangle frequency.

The frequency of a triangle $t$ is defined as the frequency of the center of
gravity of $t$. The median over the frequencies of all the triangles of the
surface is computed. Then, the global distortion
$D_{\operatorname{\Im}_{1}}(f,E)$ is computed for all the triangles with a
frequency lower than the median ($t\in\operatorname{\Im}_{1}$), and a second
global distortion $D_{\operatorname{\Im}_{2}}(f,E)$ is computed for the
triangles with a frequency greater that the median
($t\in\operatorname{\Im}_{2}$). In such a way, each distortion measure $E$
contributes the following two features: $D_{\operatorname{\Im}_{1}}(f,E)$ and
$D_{\operatorname{\Im}_{2}}(f,E)$. We will denote these features by $E_{1}(f)$
and $E_{2}(f)$, for short. That is,

$E_{i}(f)=D_{\operatorname{\Im}_{i}}(f,E),~{}i=1,2\,,$ (10)

where $D_{\operatorname{\Im}_{i}}$ is defined according to (9).

To summarize, we measure global distortions over the two subsets of triangles
and use the obtained quantities as a shape descriptors of spectogram surfaces.
This approach has the following advantages over the previously available
distortion-based models for shape analysis [30, 25]:

  1. 1.

A wider set of distortion measures is used.

  2. 2.

The overall number of features is further increased by dividing distortions
into the low and high frequencies.

  3. 3.

The method operates on triangular meshes instead of tetrahedral meshes.
Compared with the volumetric method of [25], extracting features in our
algorithm has a lower computation cost333We use triangular meshes because our
data is represented by disc-topology surfaces. However, it is possible to
represent this data by tetrahedral meshes and to estimate volumetric
distortions of these meshes. See the Discussion section for more details..

A high-level overview of the proposed model is summarized in Figure 6.

Figure 6: High-level architecture of the proposed model. The model is composed
of three main steps: signal processing, geometric learning and classification
by means of machine learning.

###  III-I Baseline

A baseline (i.e. a reference mode) has to be created for comparison the model
created with. A different approach has been selected for this purpose, based
on a set of features has been handcrafted. Twelve Mel Frequency Cepstrum
Coefficients (MFCCs) were extracted from the audio files. MFCC is the most
widely used feature extraction method in automatic speech recognition [42]. In
the feature extraction phase, 6 statistical parameters have been extracted
from each of the 12 MFCC coefficients as follows: mean, standard deviation,
min, max, mean of the absolute difference, standard deviation of the absolute
difference. Altogether, 72 features.

The reference model has been applied to all classifiers, with the same
training cost in the case of the proposed model. The models that we use for
comparison have been trained and tested on the same train/test subdivision of
the data. The proposed model was first applied with each distortion measure
separately, to assess the efficiency of each of them. Combination of the
distortion measures were then tested.

Finally, the MFCC-based model has been combined with the proposed model (based
on distortion measures). For each recording, 88 features have been computed:
the 16 features based on distortions measures (2 for each distortion), and the
72 features based on MFCC coefficients. As the number of features increased
significantly, a feature selection step has been applied, based on the ranking
of features, determined by implementation in the Random Forest classifier.
Altogether, 45 features have been selected.

##  IV Implementation

###  IV-A Training

A total of two classification tasks have been conducted: a multi-class
classification, with the 5 pathologies and the healthy recordings, and a
binary classification, for each of the 5 pathologies against the class of
healthy recordings.

The dataset was subdivided into 80% training set and 20% test set. For each
task, several classifiers were experimented: Logistic regression (LR), support
vector machine (SVM), Random forest (RF), $K$ nearest neighbors (KNN), and
AdaBoost (AB). For all of these models we used the 16 engineered features. For
each model, hyper-parameters such as the number of estimators or number of
neighbors were optimized using 5–fold cross validation. A large random grid of
hyper-parameters was searched for (see Supplements). In the case of the multi-
class classification, the performance measure used for optimization was the
accuracy, whereas for the binary tasks the area under the receiver operating
characteristics curve (AUROC) was used. A weight has been assigned to each
class, inversely proportional to the class frequencies in the training set.

For each iteration of the cross-fold, training examples were divided into
training and validation set by stratifying among patients, which means several
recordings from the same patient are always in the same set.

###  IV-B Evaluation

All the models were trained on the same test set. That is, for all the models,
the database was split into the same training and test subsets.

The following metrics were used for the performance evaluation:

$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}\,,$ $Recall=\frac{TP}{TP+FN}\,,$
$Jaccard=\frac{TP+TN}{2*(P+N)-(TP+TN)}\,,$

where $TP,TN,FP,FN$ are respectively the the True Positives, True Negatives,
False Positives and False Negatives. $P$ is the number of positives samples,
$N$ is the number of negatives samples.
The area under the ROC curve, $AUROC$, is computed.

##  V Results

###  V-A MFCC baseline model

The results of the MFCC baseline model are summarized in Table I,

TABLE I: Results of the binary classification utilizing the MFCC-based model. | URTI | COPD | Bronchiectasis | Pneumonia | Bronchiolitis
---|---|---|---|---|---
$AUROC$ | 0.82 | 0.98 | 0.95 | 0.90 | 0.81
Accuracy | 0.83 | 0.99 | 0.96 | 0.91 | 0.70
Recall | 0.83 | 0.99 | 0.96 | 0.91 | 0.70

###  V-B Classification

The results obtained by each classifier are summarized in Table II and the
results of the binary classification appear in Table III. For clarity, we
report only the best results of the binary classifications in Table III.

TABLE II: Results obtained the 5 classifiers on the test set. The best score for each metric is underlined. Classifier | Accuracy | Jaccard score | Recall
---|---|---|---
KNN | 0.86 | 0.76 | 0.86
LR | 0.88 | 0.80 | 0.89
AB | 0.88 | 0.80 | 0.89
RF | 0.88 | 0.79 | 0.88
SVM | 0.88 | 0.78 | 0.88
TABLE III: Score for binary classifiers, for the 5 different pathologies. | URTI | COPD | Bronchiectasis | Pneumonia | Bronchiolitis
---|---|---|---|---|---
$AUROC$ | 0.89 | 0.99 | 1.00 | 0.87 | 1.00
Accuracy | 0.90 | 0.97 | 1.00 | 0.87 | 1.00
Recall | 0.90 | 0.97 | 1.00 | 0.87 | 1.00
TABLE IV: Results of the combined model, for classifying 5 pathologies. | URTI | COPD | Bronchiectasis | Pneumonia | Bronchiolitis
---|---|---|---|---|---
$AUROC$ | 0.93 | 1.00 | 1.00 | 1.00 | 1.00
Accuracy | 0.92 | 1.00 | 1.00 | 0.90 | 1.00
Recall | 0.92 | 1.00 | 1.00 | 0.90 | 1.00

Ranking of the features according to their importance, as determined by Random
Forest classifier, is depicted in Figure 7.

Figure 7: Feature importance, determined by using a Random Forest classifier.
The ranking of the features allows us to understand what are the most useful
features in the classification, and yields further insight into the model.

The proposed model obtains a better $AUROC$ than the baseline models for
almost all the binary tasks. For differentiation of pneumonia pathology from
the rest of diseases, the proposed model yields a lower $AUROC$ value than the
baseline model (0.87 vs 0.90).

###  V-C Combination of models

After comparing the proposed models with the MFCC-based model, a combination
of the models has been implemented as follows:

For each recording the feature vectors from the 2 models have been
concatenated, leading to a total of 88 features. Figure 8 presents the ranking
of the features, for each of the 5 binary tasks and the multiclass task of
identifying the 5 pathologies. Although there are 16 distortion measures
features and 72 MFCC features, for most of the pathologies the occurrence of
the distortion measure features is relatively high. In particular, there are
six distortion measures out of ten most highly-ranked features for the
Bronchiectasis and URTI pathologies. Likewise, distortion measures appear
among the four most highly-ranked features used in classification of the
Bronchiolitis and COPD diseases. Indeed, for this pathology the MFCC-based
model outperformed the proposed model. However, if our task is to identify a
Pneumonia lung sound, then only a single distortion measure appears in the
feature ranking list. Indeed, for this pathology the MFCC based model
outperformed the proposed hybrid model.

Figure 8: Ranking of the 10 best features of the combined model, determined by
means of the performance of the Random Forest classifier. On panel (a), the
ranking is presented for the binary task of identifying Bronchiectasis. On
panel (b), for the Bronchiolitis. On panel (c), for COPD. On panel (d), for
Pneumonia. On panel (e), for URTI. Finally on panel (f), the ranking is
presented for the multi class task of classifying the 5 pathologies.

Secondly, Table IV presents the results of the combined model for the binary
tasks. This hybrid model receives a $AUROC$ of 1.00 for Bronchiectasis,
Bronchiolitis and COPD diseases. To summarize, the combination of both
approaches with the feature selection step improved the results.

Finally, Figure 9 presents the accuracy on the test set of the combined model
for the multiclass classification, over the number of features selected. The
results are computed by using the test set. The most accurate classification
results are achieved with 40 features. The accuracy decreases with higher
number of features, as the model begin to overfit the data.

Figure 9:  Accuracy obtained in classification of the pulmonary sounds on test
set versus the number of features, used in the combined model. Note the
decrease in the accuracy, due to over-fitting.

##  VI Discussion

The success penetration of innovative machine learning methodologies into
certain disciplines of medical diagnostics, such as radiology, has increased
the interest in achieving an automatic process for analysis of lung sounds.

This a process can be used for second opinion in the diagnosis of pulmonary
diseases, it can be useful for monitoring patients in critical care units and
may even substitute for medical experts in mass screening of malaria, as
example, in parts of the world where there is severe deficiency in medical
expertise and far from sufficient manpower.

With these goals in mind, a new approach has been proposed in this work, based
on geometric properties of the spectrographic representation of lung sounds.
This model outperformed MFCC-based model in four out of five pathologies for
which we had access to sufficient data. Further, as we have shown, the two
approaches can be combined by concatenating feature vector, or by linear
regression between the models output probabilities. These approaches can even
be incorporated into a deep learning framework which has emerged as the most
widely used approach in the field of machine learning. Such a hybrid paradigm
has recently attracted the attention of the community of AI [43].

It should be noted that our model matches the performance of the existing deep
learning methods, while requiring much fewer data samples for training.
Furthermore, compared with purely data-driven approaches, our handcrafted
features allow a better theoretical understanding of the sound classification
problem. In particular, by analyzing properties of simplicial maps we can
identify distortion-preserving transformations of audio signals in the time-
frequency domain (see Figure 10 for the illustration).

Figure 10:  This figure illustrates that, in certain scenarios, distortion
measures are invariant under translations of spectrograms in the time-
frequency domain. We show from left to right: a surface representation of the
given spectrogram, the same surface shifted along the time and the frequency
axes. Colors encode the MIPS distortion (5) induced by the projection
$f_{\text{proj}}:(x,y,z)\mapsto(x,y,0)$. It is easy to show that
$f_{\text{proj}}$ induces equal global distortions over the depicted surfaces,
for any local distortion $E(\sigma_{1},\sigma_{2})$. Indeed, triangles that
belong to the flat regions have zero distortions under $f_{\text{proj}}$,
whereas triangles of the peaks produce the same amount of a local distortion
because $E(\sigma_{1},\sigma_{2})$ is invariant under translations.

With the available, rather limited volume of dataset, we cannot predict what
will be the optimal subset of features to select in the combined method that
incorporates MFCC and distortion measures. Indeed, the subset of features
varies from disease to disease. A larger dataset is needed to assess what is
the global optimal subset of features.

Our novel approach to classification, based on distortion measures, highlights
an interesting direction for future research would be to consider a higher
dimensional distortion measures. In particular, a spectrogram surface can be
represented by the tetrahedral volume enclosed by that surface and the plane
$z=0$. In such a case, sound spectrograms can be characterized by 3D
distortions induced by mapping tetrahedral meshes into canonical domains.
Although the tetrahedral approach has a higher computational cost, it can
yield potentially more accurate results because volumetric distortions can
detect both the changes made to the boundary surface and the changes made to
the interior volume.

It should be also interesting to combine our model with various types of shape
descriptors, such as the metric and spectral geometrical features, listed in
Section III-E.

We also plan to examine more methods for discretizing the spectrogram images.
In particular, a curvature based method [44] can be used for a more accurate
sampling of spectrogram images and for building triangle meshes with an
optimal number of vertices.

Finally, we stress that our approach to the classification of one-dimensional
signals is also applicable to higher-dimensional signals. Distortion measures
can be extended in a straightforward manner to $\mathbb{R}^{n}$ and to
piecewise linear manifolds embedded in $\mathbb{R}^{n}$, for any $n\geq 2$
444See [27], for $n$-dimensional variants of distortion measures introduced in
Section III-G. Indeed, if
$\boldsymbol{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a simplicial map
and $s$ is an $n$-dimensional simplex, then a local distortion of
$\boldsymbol{f}$ over $s$ can be expressed by a function
$E(\sigma_{1},\sigma_{2},\dots,\sigma_{n})$, where $\sigma_{i}$ denotes the
$i^{\text{th}}$ singular value of the Jacobian matrix
$d\boldsymbol{f}_{s}\in\mathbb{R}^{n\times n}$. So that our distortion-based
analysis of surfaces has been extended to $m$-manifolds embedded in
$\mathbb{R}^{n}$ and to their discrete representations, for any $2\leq m\leq
n$. For instance, consider a simultaneous recording of different time-varying
signals such as pulmonary sounds, heart rate, oxygen saturation, body
plethysmography, etc. Instead of computing spectrograms for each signal
separately, one can represent a $n$-channel data stream by a 2-manifold
embedded in $\mathbb{R}^{n}$. The obtained manifold can be discretized using
the sampling method of [44] and a Delaunay-based algorithm for triangulation.
An extension of our approach to manifolds thereby allows a more general
analysis of multichannel biomedical data, collected from various devices. We
therefore consider other applications of the proposed distortion-based model
in related fields of biomedical signal processing, medical imaging and voice
recognition.

## References

  * [1] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” Medical image analysis, vol. 42, pp. 60–88, 2017.
  * [2] M. Topalovic, N. Das, P.-R. Burgel, M. Daenen, E. Derom, C. Haenebalcke, R. Janssen, H. A. Kerstjens, G. Liistro, R. Louis, et al., “Artificial intelligence outperforms pulmonologists in the interpretation of pulmonary function tests,” European Respiratory Journal, vol. 53, no. 4, 2019.
  * [3] D. Chamberlain, R. Kodgule, D. Ganelin, V. Miglani, and R. R. Fletcher, “Application of semi-supervised deep learning to lung sound analysis,” in 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 804–807, 2016.
  * [4] D. Emmanouilidou, K. Patil, J. West, and M. Elhilali, “A multiresolution analysis for detection of abnormal lung sounds,” in 2012 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pp. 3139–3142, 2012.
  * [5] E. W. Pefura-Yone, A. P. Kengne, A. D. Balkissou, C. G. Magne-Fotso, M. Ngo-Yonga, J. R. Boulleys-Nana, N. R. Efe-de Melingui, P. I. Ndjeutcheu-Moualeu, C. L. Mbele-Onana, E. C. Kenmegne-Noumsi, et al., “Prevalence of obstructive lung disease in an african country using definitions from different international guidelines: a community based cross-sectional survey,” BMC research notes, vol. 9, no. 1, pp. 1–10, 2016\.
  * [6] A. Imran, I. Posokhova, H. N. Qureshi, U. Masood, M. S. Riaz, K. Ali, C. N. John, M. I. Hussain, and M. Nabeel, “Ai4covid-19: Ai enabled preliminary diagnosis for covid-19 from cough samples via an app,” Informatics in Medicine Unlocked, vol. 20, p. 100378, 2020.
  * [7] G. Chambres, P. Hanna, and M. Desainte-Catherine, “Automatic detection of patient with respiratory diseases using lung sound analysis,” in 2018 International Conference on Content-Based Multimedia Indexing (CBMI), pp. 1–6, IEEE, 2018.
  * [8] D. Bogdanov, N. Wack, E. Gómez Gutiérrez, S. Gulati, H. Boyer, O. Mayor, G. Roma Trepat, J. Salamon, J. R. Zapata González, X. Serra, et al., “Essentia: An audio analysis library for music information retrieval,” in Britto A, Gouyon F, Dixon S, editors. 14th Conference of the International Society for Music Information Retrieval (ISMIR); 2013 Nov 4-8; Curitiba, Brazil.[place unknown]: ISMIR; 2013. p. 493-8., International Society for Music Information Retrieval (ISMIR), 2013.
  * [9] B. Rocha, D. Filos, L. Mendes, I. Vogiatzis, E. Perantoni, E. Kaimakamis, P. Natsiavas, A. Oliveira, C. Jácome, A. Marques, et al., “A respiratory sound database for the development of automated classification,” in International Conference on Biomedical and Health Informatics, pp. 33–37, Springer, 2017.
  * [10] A. Kandaswamy, C. Kumar, R. Ramanathan, S. Jayaraman, and N. Malmurugan, “Neural classification of lung sounds using wavelet coefficients,” Computers in Biology and Medicine, vol. 34, no. 6, pp. 523 – 537, 2004.
  * [11] Y. Amrulloh, U. Abeyratne, V. Swarnkar, and R. Triasih, “Cough sound analysis for pneumonia and asthma classification in pediatric population,” in 2015 6th International Conference on Intelligent Systems, Modelling and Simulation, pp. 127–131, IEEE, 2015.
  * [12] B. Sankur, Y. P. Kahya, E. Çağatay Güler, and T. Engin, “Comparison of ar-based algorithms for respiratory sounds classification,” Computers in Biology and Medicine, vol. 24, no. 1, pp. 67 – 76, 1994.
  * [13] A. Savitzky and M. J. E. Golay, “Smoothing and differentiation of data by simplified least squares procedures,” Analytical Chemistry, vol. 36, pp. 1627–1639, Jan. 1964.
  * [14] J. Chen, P. Jönsson, M. Tamura, Z. Gu, B. Matsushita, and L. Eklundh, “A simple method for reconstructing a high-quality ndvi time-series data set based on the savitzky–golay filter,” Remote Sensing of Environment, vol. 91, no. 3, pp. 332 – 344, 2004.
  * [15] N. S. Haider, R. Periyasamy, D. Joshi, and B. K. Singh, “Savitzky-Golay Filter for Denoising Lung Sound,” Brazilian Archives of Biology and Technology, vol. 61, 00 2018.
  * [16] H. Edelsbrunner et al., Geometry and topology for mesh generation.  Cambridge University Press, 2001.
  * [17] M. Reuter, F.-E. Wolter, and N. Peinecke, “Laplace–beltrami spectra as ‘shape-dna’of surfaces and solids,” Computer-Aided Design, vol. 38, no. 4, pp. 342–366, 2006.
  * [18] R. M. Rustamov, “Laplace-beltrami eigenfunctions for deformation invariant shape representation,” in Proceedings of the fifth Eurographics symposium on Geometry processing, pp. 225–233, 2007.
  * [19] A. M. Bronstein, M. M. Bronstein, and R. Kimmel, “Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching,” Proceedings of the National Academy of Sciences, vol. 103, no. 5, pp. 1168–1172, 2006.
  * [20] A. B. Hamza and H. Krim, “Geodesic object representation and recognition,” in International conference on discrete geometry for computer imagery, pp. 378–387, Springer, 2003.
  * [21] A. Elad and R. Kimmel, “On bending invariant signatures for surfaces,” IEEE Transactions on pattern analysis and machine intelligence, vol. 25, no. 10, pp. 1285–1295, 2003.
  * [22] R. M. Rustamov, M. Ovsjanikov, O. Azencot, M. Ben-Chen, F. Chazal, and L. Guibas, “Map-based exploration of intrinsic shape differences and variability,” ACM Transactions on Graphics (TOG), vol. 32, no. 4, pp. 1–12, 2013.
  * [23] D. M. Boyer, Y. Lipman, E. S. Clair, J. Puente, B. A. Patel, T. Funkhouser, J. Jernvall, and I. Daubechies, “Algorithms to automatically quantify the geometric similarity of anatomical surfaces,” Proceedings of the National Academy of Sciences, vol. 108, no. 45, pp. 18221–18226, 2011.
  * [24] R. Sandhu, A. Dominitz, Y. Gao, and A. Tannenbaum, “Volumetric mapping of genus zero objects via mass preservation,” arXiv preprint arXiv:1205.1225, 2012.
  * [25] A. Naitsat, S. Cheng, X. Qu, X. Fan, E. Saucan, and Y. Y. Zeevi, “Geometric approach to detecting volumetric changes in medical images,” Journal of Computational and Applied Mathematics, vol. 329, pp. 37–50, 2018.
  * [26] K. Hormann, B. Lévy, and A. Sheffer, “Mesh parameterization: Theory and practice,” 2007.
  * [27] A. Naitsat, G. Naitzat, and Y. Y. Zeevi, “On inversion-free mapping and distortion minimization,” 2020.
  * [28] A. Naitsat, Y. Zhu, and Y. Y. Zeevi, “Adaptive block coordinate descent for distortion optimization,” Computer Graphics Forum, vol. 39, no. 6, pp. 360–376, 2020.
  * [29] M. S. Floater, “One-to-one piecewise linear mappings over triangulations,” Math. Comp., vol. 72, pp. 685–696, Oct. 2002.
  * [30] M. Ben-Chen and C. Gotsman, “Characterizing shape using conformal factors,” in Proceedings of the 1st Eurographics Conference on 3D Object Retrieval, 3DOR ’08, (Goslar, DEU), p. 1–8, Eurographics Association, 2008\.
  * [31] O. Sorkine and M. Alexa, “As-rigid-as-possible surface modeling,” in Proceedings of the Fifth Eurographics Symposium on Geometry Processing, SGP ’07, (Goslar, DEU), p. 109–116, Eurographics Association, 2007.
  * [32] J. Smith and S. Schaefer, “Bijective parameterization with free boundaries,” ACM Trans. Graph., vol. 34, July 2015.
  * [33] A. Naitsat, E. Saucan, and Y. Y. Zeevi, “Geometric approach to estimation of volumetric distortions,” in Proceedings of the 11th Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications: Volume 1: GRAPP, GRAPP 2016, (Setubal, PRT), p. 105–112, SCITEPRESS - Science and Technology Publications, Lda, 2016.
  * [34] O. Sorkine, D. Cohen-Or, R. Goldenthal, and D. Lischinski, “Bounded-distortion piecewise mesh parameterization,” in IEEE Visualization, 2002. VIS 2002., pp. 355–362, 2002.
  * [35] A. Naitsat, E. Saucan, and Y. Y. Zeevi, “Computing quasi-conformal maps in 3d with applications to geometric modeling and imaging,” in 2014 IEEE 28th Convention of Electrical & Electronics Engineers in Israel (IEEEI), pp. 1–5, IEEE, 2014.
  * [36] K. HORMANN, “Mips : An efficient global parametrization method,” Curve and Surface Design: Saint-Malo, pp. 153–162, 1999.
  * [37] X.-M. Fu, Y. Liu, and B. Guo, “Computing locally injective mappings by advanced mips,” ACM Trans. Graph., vol. 34, July 2015.
  * [38] J. M. P. Degener and R. Klein, “An adaptable surface parameterization method,” 2003.
  * [39] D. Ezuz, J. Solomon, and M. Ben-Chen, “Reversible harmonic maps between discrete surfaces,” ACM Trans. Graph., vol. 38, Mar. 2019.
  * [40] S. J. G. J. S. Sander, Pedro V. and H. Hoppe, “Signal-specialized parameterization,” 2002.
  * [41] W. Abikoff, “The uniformization theorem,” The American Mathematical Monthly, vol. 88, no. 8, pp. 574–592, 1981.
  * [42] D. O’Shaughnessy, “Invited paper: Automatic speech recognition: History, methods and challenges,” Pattern Recognition, vol. 41, no. 10, pp. 2965 – 2979, 2008.
  * [43] D. Rueckert and J. A. Schnabel, “Model-based and data-driven strategies in medical image computing,” Proceedings of the IEEE, vol. 108, no. 1, pp. 110–124, 2019.
  * [44] E. Saucan, E. Appleboim, and Y. Y. Zeevi, “Sampling and reconstruction of surfaces and higher dimensional manifolds,” Journal of Mathematical Imaging and Vision, vol. 30, no. 1, pp. 105–123, 2008.

## Supplements

### Supplementary Note 1

For Random Forests classifier, the grid focused on:

  * •

Number of estimators (100, 110, 120, 150, 200, 250, 300).

  * •

Number of features to consider at every split (could be all features, or just
the square of overall features).

  * •

Maximum number of levels in tree (from 10 to 110, with pace of 10).

  * •

Minimum number of samples required to split a node (2,5,10).

  * •

Minimum number of samples required at each leaf node (1,2,4).

  * •

Enable/Disable bootstrap.

For SVM classifier, the grid contained:

  * •

Regularization parameter, C. The strength of the regularization is inversely
proportional to C. An exponential search was used, from 10-9 to 1013.

  * •

Radial Basis Function (RBF) Kernel was used.

  * •

Degree of the kernel function (relevant only for polynomial kernel): 1, 2, 3,
4, 5, 10.

  * •

Value of gamma. An exponential search was used, from 10-9 to 1013.

  * •

Whether to use the shrinking heuristic.

Regarding the KNN classifier, the random grid has the following features:

  * •

Number of neighbors to use for k-neighbors queries: small numbers were tried,
like 3, 5, 7, 11 and then numbers from 41 to 131, with pace of 10.

  * •

Weight function used in prediction. Whether uniform, or distance (meaning
weight points by the inverse of their distance).

  * •

The algorithm used. It could be ball tree, kd tree, or brute force.

  * •

The power parameter for Minkowski metric, p: 1 for manhattan distance, 2 for
Euclidean distance.

In the case of Logistic Regression model:

  * •

The norm used in the penalization. It can be L1, L2, or elasticnet.

  * •

Tolerance for stopping criteria. An exponential search was used, from
$10^{-9}$ to $10^{13}$.

  * •

Inverse of regularization strength, C: an exponential search was used, from
$10^{-9}$ to $10^{13}$.

  * •

The algorithm to use in the optimization problem. There are liblinear, lbfgs,
sag, saga and newton-cg.

Adaboost combines a series of weak learners, with the aim of creating an
improved classifier. The weak learners vote for the final prediction label,
and a majority voting is performed. The grid search for the Adaboost
classifier is constitued of:

  * •

Number of estimators: A linear search was used, from 50 to 200 with step of
10.

  * •

Learning rate: An exponential search was used, from $10^{-9}$ to 10.

  * •

Base estimator: RandomForest, DecisionTree, SVM, LogisticRegression.

### Supplementary: Delaunay triangulation

A Delaunay triangulation of a finite set $P\subset{\mathbb{R}}^{3}$ is defined
as a collection $I$ of triangles, such that:

  * •

$\mathit{conv}(\operatorname{\mathcal{V}})=\cup_{T\epsilon\Im}T$

  * •

$P=\cup_{T\epsilon\Im}V(T)$

  * •

For every distinct pair $T,U\in\Im$, the intersection $T\cap U$ is either a
vertex, an edge, or empty.

For a finite set $P\subset{\mathbb{R}}^{3}$, the Delaunay triangulation is
unique if the points are in a general position, which means no 5 points are
cospherical. The assumption of unique Delaunay triangulation has been made in
this work.