# Post-merger gravitational-wave signal from neutron-star binaries: a new look
at an old problem

Konrad Topolski Institut für Theoretische Physik, Goethe Universität, Max-von-
Laue-Str. 1, 60438 Frankfurt am Main, Germany  Samuel D. Tootle Institut für
Theoretische Physik, Goethe Universität, Max-von-Laue-Str. 1, 60438 Frankfurt
am Main, Germany  Luciano Rezzolla Institut für Theoretische Physik, Goethe
Universität, Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany  Frankfurt
Institute for Advanced Studies, Ruth-Moufang-Str. 1, 60438 Frankfurt am Main,
Germany  School of Mathematics, Trinity College, Dublin 2, Ireland

###### Abstract

The spectral properties of the post-merger gravitational-wave signal from a
binary of neutron stars encodes a variety of information about the features of
the system and of the equation of state describing matter around and above
nuclear saturation density. Characterising the properties of such a signal is
an “old” problem, which first emerged when a number of frequencies were shown
to be related to the properties of the binary through “quasi-universal”
relations. Here we take a new look at this old problem by computing the
properties of the signal in terms of the Weyl scalar $\psi_{4}$. In this way,
and using a database of more than 100 simulations, we provide the first
evidence for a new instantaneous frequency, $f^{\psi_{4}}_{0}$, associated
with the instant of quasi time-symmetry in the postmerger dynamics, and which
also follows a quasi-universal relation. We also derive a new quasi-universal
relation for the merger frequency $f^{h}_{\rm mer}$, which provides a
description of the data that is four times more accurate than previous
expressions while requiring fewer fitting coefficients. Finally, consistently
with the findings of numerous studies before ours, and using an enlarged
ensamble of binary systems we point out that the $\ell=2,m=1$ gravitational-
wave mode could become comparable with the traditional $\ell=2,m=2$ mode on
sufficiently long timescales, with strain amplitudes in a ratio
$|h^{21}|/|h^{22}|\sim 0.1-1$ under generic orientations of the binary, which
could be measured by present detectors for signals with large signal-to-noise
ratio or by third-generation detectors for generic signals should no collapse
occur.

Neutron stars (1108), Gravitational waves (678), Compact binary stars (283)

††software: Einstein Toolkit (Haas & et al., 2020), Carpet (Schnetter et al.,
2004), FIL (Most et al., 2019a), FUKA (Papenfort et al., 2021), Kadath
(Grandclement, 2010), Watpy (The CoRe Collaboration, 2022), Kuibit (Bozzola,
2021), Mathematica (Wolfram Research, 2020)

##  1 Introduction

The observation of a gravitational-wave (GW) signal from the binary neutron-
star (BNS) merger event GW170817 (The LIGO Scientific Collaboration & The
Virgo Collaboration, 2017), and the detection of an electromagnetic (EM)
counterpart, has testified the enormous potential of GW astronomy. Starting
from early works with simplified equations of state (EOSs) (see, e.g., Shibata
et al., 2005; Anderson et al., 2008; Liu et al., 2008; Baiotti et al., 2008;
Hotokezaka et al., 2011), increasingly more comprehensive simulations of these
events, which involve an ever more detailed description of the microphysics
(Bauswein et al., 2019; De Pietri et al., 2019; Gieg et al., 2019; Tootle et
al., 2022; Most et al., 2022; Camilletti et al., 2022; Ujevic et al., 2023),
of the magnetic-field evolution (Rezzolla et al., 2011; Dionysopoulou et al.,
2013; Ciolfi et al., 2019; Sun et al., 2022; Zappa et al., 2023), and its
amplification (Kiuchi et al., 2015; Palenzuela et al., 2022; Chabanov et al.,
2023), and of transport of neutrinos (Foucart et al., 2022; Zappa et al.,
2023), allow one to make predictions from the early inspiral up to the long-
term evolution of the postmerger remnant (De Pietri et al., 2020; Kiuchi et
al., 2022). During each stage in the evolution of the binary, the features of
the GW and EM signals change in a characteristic manner, encoding information
on the properties of the constituent neutron stars and of the hypermassive
neutron star (HMNS) produced after the merger and, hence, on the governing
EOS.

Characterising the properties of the post-merger GW signal is a rather “old”
problem, which has first emerged when a number of peculiar frequencies were
shown to be related with the properties of the binary through quasi-universal
relations, i.e., relations that are almost independent of the specific EOS.
These relations have been suggested for the GW frequency at merger $f_{\rm
mer}$ (Read et al., 2013; Bernuzzi et al., 2014; Takami et al., 2015; Rezzolla
& Takami, 2016; Most et al., 2019; Bauswein et al., 2019; Weih et al., 2020;
Gonzalez et al., 2022), the dominant frequency in the postmerger spectrum
$f_{2}$ (see, e.g., Oechslin & Janka, 2007; Bauswein & Janka, 2012; Read et
al., 2013; Rezzolla & Takami, 2016; Gonzalez et al., 2022), and other
frequencies identifiable in the transient period right after the merger
(Bauswein & Stergioulas, 2015; Takami et al., 2015; Rezzolla & Takami, 2016).
Fits to these quasi-universal relations have been employed in a number of
studies (see, e.g., (Bauswein et al., 2016; Baiotti & Rezzolla, 2017) for some
reviews). These EOS-insensitive relations can help enormously in constraining
the EOS of matter at nuclear densities, marking the possible appearance of
phase transitions (Most et al., 2019; Weih et al., 2020; Liebling et al.,
2021; Prakash et al., 2021; Fujimoto et al., 2023; Tootle et al., 2022; Espino
et al., 2023), inform waveform models (Bose et al., 2018; Breschi et al.,
2019); however, see Raithel & Most (2022) for possible violations of these
relations.

Another relatively “old” problem in the characterisation of the GW signal from
BNS is the one about the relative weight of the lower-order multipole
$\ell=2,m=1$. Numerical simulations have highlighted that the HMNS can be
subject to a nonaxisymmetric instability that powers the growth of
$\ell=2,m=1$ mode of the rest-mass density distribution and, hence, of the
corresponding GW signal (see, e.g., East et al., 2015; Lehner et al., 2016a,
b; Radice et al., 2016; East et al., 2019; Papenfort et al., 2022). This mode
can already be seeded by the initial asymmetry of the system in the unequal-
mass case, or develop by the shearing of the contact layers of the binary
constituents upon merger. While the $\ell=2,m=2$ GW mode is the primary
contributor to the GW signal, it is damped faster than the other modes leading
to interesting secular behaviours.

Figure 1: 3D representation of the complete information in the GW signal from
a representative binary (Tootle et al., 2022). The left panel shows the
$\ell=2,m=2$ mode of the GW signal $\psi_{4}(t)$ (light red) and its amplitude
$|\psi_{4}(t)|$ (dark red), the instantaneous frequency $f_{{}_{\rm GW}}(t)$
(black), and the power spectral density (PSD) $\sqrt{2}f\,\tilde{\psi}_{4}$
(blue) as a function of the frequency $f$. Also indicated are the frequency at
merger $f_{\rm mer}^{\psi_{4}}$, the frequency at quasi time-symmetry
$f^{\psi_{4}}_{0}$, the dominant frequency of the HMNS emission $f_{2}$, and
the corresponding times where these frequencies appear. The right panel shows
the same quantities but when computed in terms of the GW strain.

We here take a new look at both of these old problems by considering the
spectral properties of the GW signal when computed in terms of the Weyl scalar
$\psi_{4}$. In this way, we are able to find three novel features that can
enrich our understanding of the GW signal from BNS mergers. In particular, we
first highlight the presence of a new instantaneous frequency, which we dub
$f^{\psi_{4}}_{0}$, that can be associated with the instant of quasi time-
symmetry in the postmerger dynamics. Interestingly, we find that a quasi-
universal relation exists for $f^{\psi_{4}}_{0}$ as a function of the tidal
deformability $\kappa_{2}^{{}_{T}}$ and of the binary mass ratio $q$. Second,
by employing a large number of BNS simulations, some of which are taken from
the CoRe database (Gonzalez et al., 2022), we obtain a new quasi-universal
relation for $f_{\rm mer}$ as a function of $\kappa_{2}^{{}_{T}}$ and $q$ that
not only requires a smaller number of coefficients, but also provides a more
accurate description of the data. Finally, as already suggested in (Papenfort
et al., 2022), we provide evidence that the $\ell=2,m=1$ GW mode could become
the most powerful mode on secular timescales after the merger.

##  2 Numerical and physical framework

Our analysis is based on the GW signal computed via numerical simulations of
BNS mergers in full general relativity computed with the codes described in
(Radice et al., 2014a, b; Most et al., 2019a, b; Papenfort et al., 2021;
Tootle et al., 2021) and using a number of different EOSs (see below). In
addition, we employ part of the data contained in the CoRe database (Gonzalez
et al., 2022), from where we select only simulations with the highest-
resolution. The combined data of $118$ irrotational binaries covers the range
$q:=M_{2}/M_{1}\in[0.485,1]$ in the mass ratio,
$M:=M_{1}+M_{2}\in[2.4,3.33]~{}M_{\odot}$ in the total ADM mass at infinite
separation, and $\kappa_{2}^{{}_{T}}\in[33,458]$ in the tidal deformability.
The dataset comprises a variety of EOSs including some with quark matter
(Prakash et al., 2021; Logoteta, Domenico et al., 2021; Alford et al., 2005;
Demircik et al., 2022; Tootle et al., 2022).

A crucial role in our analysis is played by the use of the Weyl scalar
$\psi_{4}$ in place of the standard dimensionless strain polarisations
$h_{+,\times}$. The two quantities are mathematically equivalent and related
by two time derivatives (i.e., $\psi_{4}=\partial^{2}_{t}(h_{+}-ih_{\times})$;
see (Bishop & Rezzolla, 2016) for a review). However, while $\psi_{4}$ is
computed from the simulations, $h_{+,\times}$ are obtained after a nontrivial
double time integration (the transformation from $h_{+,\times}$ to is
$\psi_{4}$ trivial as it involves derivatives and not integrals; see (Calderon
Bustillo et al., 2022) for a data-analysis framework based on $\psi_{4}$,
which can obviously be employed for all types of compact-object binaries).
More importantly, the evolution of the GW frequency from $\psi_{4}$ is less
rapid than from the strain, i.e., $\partial_{t}\ln f^{\psi_{4}}_{{}_{\rm
GW}}(t)\ll\partial_{t}\ln f^{h}_{{}_{\rm GW}}(t)$, thus making it easier and
more robust to characterise the features of the $\psi_{4}$ GW signal. In this
sense, while $\psi_{4}$ and $h_{+,\times}$ are related by simple time
derivatives, the analysis carried out with the former does provide additional
information as it allows for the determination of properties that are harder
to capture with the latter.

##  3 Old and new frequencies

Figure 1 reports the complete information of the GW signal from a
representative binary in our sample. Using a 3D representation, we report on
the left the $\ell=2,m=2$ mode of the GW signal $\psi_{4}(t)$ (light red) and
its amplitude $|\psi_{4}(t)|$ (dark red), the instantaneous frequency
$f^{\psi_{4}}_{{}_{\rm GW}}(t)$ (black), and the power spectral density (PSD)
$\sqrt{2}f\,\tilde{\psi}_{4}$ (blue) as a function of the frequency $f$ (see
Rezzolla & Takami, 2016, for details on the definition). Also indicated are
the three main frequencies in our analysis: the frequency at merger $f_{\rm
mer}^{\psi_{4}}$, i.e., the GW frequency at the first maximum of $|\psi_{4}|$,
the frequency at quasi time-symmetry $f^{\psi_{4}}_{0}$, i.e., the GW
frequency at the first minimum of $|\psi_{4}|$, and the dominant frequency of
the HMNS emission $f_{2}^{\psi_{4}}$. To help the eye, we also mark with lines
the corresponding times $t^{\psi_{4}}_{\rm mer}$ (dashed), $t^{\psi_{4}}_{0}$
(dotted), and frequencies (dashed, dotted and dot-dashed respectively). The
right panel of Fig. 1 shows the same quantities but when computed from the
strain. By comparing the black lines in the left and right panels it is
straightforward to realise that the variation of $f^{h}_{{}_{\rm GW}}(t)$ is
much larger than that in $f^{\psi_{4}}_{{}_{\rm GW}}(t)$ over the same
interval of $\sim 1\,$ms after the merger. It is this very rapid change in
$f^{h}_{{}_{\rm GW}}(t)$ that makes the identification of $f^{h}_{0}$
extremely difficult, if not impossible. Note also that while in both
representations $f_{\rm mer}<f_{2}<f_{0}$, the numerical values of the various
quantities are similar but not identical. However, $f_{2}^{\psi_{4}}\simeq
f_{2}^{h}$ to very good precision (the largest differences are $\lesssim 4\%$)
simply because this frequency is relative to a mostly monochromatic GW signal;
hence, hereafter we simply assume $f_{2}^{\psi_{4}}=f_{2}^{h}=:f_{2}$.
Finally, because in all representations $f_{0}$ is largest frequency measured,
even a crude measure of largest frequency in the signal will serve as a first
estimate of the $f_{0}$ frequency.111From a numerical point of view, we note
that the $f_{0}$ frequency is always below $\sim 4\,{\rm kHz}$ and is
therefore much smaller than the typical sampling frequency of the $\psi_{4}$
scalar, that is $\simeq 80\\!-\\!100\,{\rm kHz}$.

Besides marking the time of the first amplitude minimum, from a physical point
of view $t^{\psi_{4}}_{0}$ corresponds to the time when the two stellar cores
have reached the minimum separation and are about to bounce-off each other. At
this instant, the corresponding amplitude of $\psi_{4}$ shows a clear minimum,
while the instantaneous GW frequency a local maximum [the discussion in the
Supplemental Material (SM) illustrates this behaviour very clearly by
employing the toy model introduced in (Takami et al., 2015)].

Figure 2: 3D representation of the GW frequencies $f^{\psi_{4}}_{\rm mer}$ and
$f^{\psi_{4}}_{0}$ as measured from the data (coloured circles) and presented
as a function of $\kappa^{{}_{T}}_{2}$ and $q$. Also reported are the best-fit
surfaces, while shown below are the relative errors of the fit in the two
principal directions. Stars mark binaries modelled with the V-QCD EOS and thus
having a strong first-order phase transition (Demircik et al., 2022; Tootle et
al., 2022).

##  4 Quasi-universal relations

We next proceed to the derivation of quasi-universal relations that can be
employed to deduce the physical properties of the binary. Following the
approach started already in (Takami et al., 2014, 2015; Rezzolla & Takami,
2016), which captures the logarithmic variation of a properly rescaled mass
and frequency, we express the relevant frequencies in terms of a power
expansion of the mass ratio $q$, i.e.,
$\log_{10}\left[({M}/M_{\odot})({f}/\rm{Hz})\right]=a_{0}+(b_{0}+b_{1}q+b_{2}q^{2})\left({\kappa_{2}^{{}_{T}}}\right)^{n}$,
where $f$ is any of the frequencies we consider (i.e., $f^{\psi_{4}}_{\rm
mer},f^{h}_{\rm mer},f^{\psi_{4}}_{0},f_{2}$), $a_{0},b_{0},b_{1},b_{2},n$ are
fitting coefficients. Hereafter, we will refer to this generic fitting
functions as $\mathcal{F}_{1}$.

Figure 2 provides a 3D representation of the measured GW frequencies
$f^{\psi_{4}}_{\rm mer}$ and $f^{\psi_{4}}_{0}$ as a function of
$\kappa^{{}_{T}}_{2}$ and $q$ (see also the SM for fits to $f^{h}_{\rm mer}$
and $f_{2}$). Also reported is the fitting surface described by
$\mathcal{F}_{1}$, with the best-fit parameters listed in Table $2$ of the SM
for all the frequencies considered. Furthermore, for each frequency we report
below the relative error of the fit in the two principal directions of the
fit, $\kappa^{{}_{T}}_{2}$ and $q$. Despite their simple form, our fits for
$f^{\psi_{4}}_{\rm mer}$ and $f^{\psi_{4}}_{0}$ capture the data very well,
showing average relative errors that are $\lesssim 1\%$ and maximal relative
errors $\lesssim 2\%$ for other than equal-mass binaries.

It is interesting to compare our functional fitting form $\mathcal{F}_{1}$ for
$f^{h}_{\rm mer}$, which needs only five fitting coefficients, with the one
proposed in (Breschi et al., 2022) for irrotational binaries, which we will
refer to as $\mathcal{F}_{2}$, and that requires twice as many coefficients.
In order to compare $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ it is first
necessary to distinguish the “pipeline”, that is, the technical procedure
employed to extract the frequencies from the data. We thus indicate with
$\mathcal{P}_{1}$ the pipeline discussed above and with $\mathcal{P}_{2}$ that
released in (Gonzalez et al., 2022). Naturally, each fitting function can be
applied to either pipeline, so that $\mathcal{F}_{1}(\mathcal{P}_{1})$
indicates the use of our fitting form to data computed with our pipeline. In
Fig. 3, we present the relative differences between the measured frequencies
for the 118 binaries considered and the corresponding values from the fit,
with different rows referring to the four possibilities.

Overall, the comparison in Fig. 3 shows that $\mathcal{F}_{1}$ leads to
smaller relative errors with a maximum residual error of $\sim 2\%$ and an an
average residual error that is between two and four times smaller than for
$\mathcal{F}_{2}$. As a cautionary note we should remark that we have
specialised the fitting $\mathcal{F}_{2}$, which is more general and can
include spinning and eccentric binaries, to the case relevant for this
comparison, namely, irrotational binaries. Hence, our conclusions apply only
to such binaries.

Figure 3: Relative differences between the measured $f_{\rm mer}^{h}$
frequencies and the corresponding values from the fit. Different rows refer to
the four different possibilities of applying the fitting function
$\mathcal{F}_{i}$ to the data pipeline $\mathcal{P}_{i}$ with $i=1,2$; arrows
indicate relative differences above $10\%$. Note how the new fitting function
and pipeline, i.e., $\mathcal{F}_{1}(\mathcal{P}_{1})$, provide errors that
are about four times smaller.

Figure 4: Evolution of the ratio of the two main GW modes
$|\psi_{4}^{21}|/|\psi_{4}^{22}|$, with the left panel showing a selected set
of binaries to highlight the behaviour for different mass ratios (colour code)
and the right one reporting all of the binaries. Black dotted horizontal lines
are used to mark the position when the two modes have equal amplitudes, while
straight coloured lines report a linear extrapolation after averaging the last
$5\,{\rm ms}$ of the evolution; note that the large majority of binaries
exhibits a growing trend that is maintained throughout the computed evolution.

##  5 Secular GW emission

The last point we cover regards the relative strengths of the $\ell=2,m=2$ and
$\ell=2,m=1$ GW modes. The importance of the latter was first pointed out in
(Paschalidis et al., 2015; Lehner et al., 2016a) and is produced by
corresponding asymmetries in the rest-mass density. The emergence of an $m=1$
deformation is well-known to occur in isolated stars (Chandrasekhar, 1969;
Shibata et al., 2000; Baiotti et al., 2007; Franci et al., 2013; Löffler et
al., 2015) that have a sufficiently large amount of rotational kinetic energy
$T$ and emerges when the ratio $T/|W|$, where $W$ is the gravitational binding
energy, exceeds a certain threshold (in a systematic analysis, Baiotti et al.,
2007, have shown that this happens for $T/|W|\gtrsim 0.25$). In such stars,
the $m=1$ mode in the rest-mass density would grow exponentially reaching
equipartition with the $m=2$ mode, and, subsequently, represent the largest
deformation. A similar phenomenology seems to be present also for the
postmerger remnant, as already hinted in (Papenfort et al., 2022), but as
shown more clearly by the numerous binaries considered here. Figure 4 reports
the evolution of the ratio of the GW amplitudes in the two modes
$|\psi_{4}^{21}|/|\psi_{4}^{22}|$, with the left panel showing a selected set
of binaries and with the right panel reporting all binaries (the raw
timeseries are smoothed over a window $\Delta t=0.5\,{\rm ms}$). While only a
few binaries in the sample reach $|\psi_{4}^{21}|/|\psi_{4}^{22}|=1$ within
the simulated time, the large majority exhibits a trend that we try to capture
by extrapolating linearly in time after averaging the last $5\,{\rm ms}$ of
the evolution. Using a colour code to distinguish binaries with different $q$,
it becomes clear that the initial strength of the $m=1$ mode is inversely
proportional to the mass ratio, so that for an extremely asymmetric binary,
i.e., $q\lesssim 0.6$, $|\psi_{4}^{21}|/|\psi_{4}^{22}|$ can be more than two
orders of magnitude larger than for equal-mass binaries. At the same time, the
initial mode-amplitude ratio does not depend on the merger dimensionless spin
(Papenfort et al., 2022). Unsurprisingly, a similar behaviour can be observed
when computing the signal-to-noise (SNR) ratio in the $m=2$ and $m=1$ modes.
Specifically, by computing a time-windowed SNR ratio we can quantify the
growing contribution of the subdominant mode in a similar fashion to estimates
given by (Lehner et al., 2016b) and find that its increases to
$\simeq\mathcal{O}(1)$ if the $m=1$ is not suppressed (see SM for details).

If confirmed by systematic, long-term evolutions, this finding would change
the standard picture in which the largest signal in the BNS postmerger is to
be expected at the $f_{2}$ frequency. Rather, the trend reported here suggests
that the most powerful feature in the PSD for long-lived HMNSs may actually
appear at a frequency $\simeq\tfrac{1}{2}f_{2}$. Because this falls in a more
favourable region of the detectors sensitivities, and assuming a detection
angle not favouring either of the modes, the corresponding signal-to-noise
ratio will grow proportionally to the ratio between detectors noise at
$\tfrac{1}{2}f_{2}$ and $f_{2}$, hence by a factor of $\sim 2$ for LIGO or
Virgo.

While promising, this prospect should be accompanied by some caveats. First,
it is possible that the growth rate may be weaker than the one estimated here.
Second, the extrapolation assumes that the HMNS will not collapse to a black
hole before reaching $|\psi_{4}^{21}|/|\psi_{4}^{22}|\sim 1$ and while this is
likely for soft EOSs and low-mass binaries, it may not happen if the EOS is
stiff and the binary massive. Third, all binaries in our sample have zero
deformation. A robust conclusion that can be inferred from the results shown
in Fig. 4 is that remnants with a long lifetime, as it was likely the case for
GW170817 (Rezzolla et al., 2018; Gill et al., 2019; Murguia-Berthier et al.,
2021), will reasonably have the $\ell=2,m=1$ as the least-damped mode. Hence,
considerable spectral power should be present at frequencies
$\frac{1}{2}f_{2}$ and $f_{2}$, with the main strain amplitudes in a ratio
$|h^{21}|/|h^{22}|\sim 0.1-1$ for generic orientations (e.g., for an
inclination of $2\arctan(1/2)\sim 53^{\circ}$ two modes have the same spin-
weighted spherical-harmonics coefficients).

##  6 Conclusion

Leveraging on a rich literature developed over the last ten years on this
subject, we have considered again the spectral properties of the signal when
computed in terms of the Weyl scalar $\psi_{4}$ rather than in terms of the GW
strain $h_{+,\times}$. Exploiting the better behaviour of $\psi_{4}$, we were
able to highlight three novel features that can be used to better infer
physical information from the detected signal.

First, by employing a large number of simulations spanning a considerable set
of EOSs and mass ratios, we have shown the existence of a new instantaneous
frequency, $f^{\psi_{4}}_{0}$, that can be associated with the instant of
quasi time-symmetry in the postmerger dynamics. This corresponds to when the
stellar cores in the merger remnant have reached their minimum separation and
are about to bounce-off each other. Just like other spectral frequencies of
the BNS GW signal, $f^{\psi_{4}}_{0}$ also follows a quasi-universal behaviour
as a function of the tidal deformability $\kappa_{2}^{{}_{T}}$ and of the
binary mass ratio $q$, for which we provide a simple and yet accurate
analytical expression. Second, we have obtained a new quasi-universal relation
for the merger frequency $f^{h}_{\rm mer}$ as a function of
$\kappa_{2}^{{}_{T}}$ and $q$. The new expression not only requires a smaller
number of fitting coefficients than alternative expressions in the literature,
but it also provides a more accurate description of the data, with a residual
error that is four times smaller on average. Finally, we have pointed out the
evidence that the $\ell=2,m=1$ could become the most powerful GW mode on
sufficiently long timescales, with strain amplitudes for the dominant modes
that are in a ratio $|h^{21}|/|h^{22}|\sim 0.1-1$. Should this mode not be
suppressed by the collapse of the HMNS to a black hole or by other dissipative
effects such as magnetic fields, considerable spectral power should be present
at frequencies $\frac{1}{2}f_{2}$, where it could be detected in conditions of
smaller signal-to-noise ratios or by third-generation detectors.

The results presented here can be improved by enlarging the number of BNS
simulations considered, by increasing the variance in the microphysical
description (e.g., including simulations with magnetic fields and neutrino
transport), by performing additional long-term evolutions, and by extending
the fitting approach to binaries with spins and eccentricity. We will explore
these extensions in future work.

Data policy. The relevant data that supports the findings of this paper is
available from the first author and can be shared upon a reasonable request.

We thank K. Chakravarti, K. Takami, C. Ecker, and C. Musolino for useful input
and discussions. Support comes from the State of Hesse within the Research
Cluster ELEMENTS (Project ID 500/10.006). LR acknowledges funding by the ERC
Advanced Grant “JETSET: Launching, propagation and emission of relativistic
jets from binary mergers and across mass scales” (Grant No. 884631). The
simulations from which parts of the used data are derived were performed on
HPE Apollo HAWK at the High Performance Computing Center Stuttgart (HLRS)
under the grants BNSMIC and BBHDISKS, and on SuperMUC at the Leibniz
Supercomputing Centre.

## References

  * Alford et al. (2005) Alford, M., Braby, M., Paris, M., & Reddy, S. 2005, Astrophys. J., 629, 969, doi: 10.1086/430902
  * Anderson et al. (2008) Anderson, M., Hirschmann, E. W., Lehner, L., et al. 2008, Phys. Rev. D, 77, 024006, doi: 10.1103/PhysRevD.77.024006
  * Baiotti et al. (2007) Baiotti, L., De Pietri, R., Manca, G. M., & Rezzolla, L. 2007, Phys. Rev. D, 75, 044023, doi: 10.1103/PhysRevD.75.044023
  * Baiotti et al. (2008) Baiotti, L., Giacomazzo, B., & Rezzolla, L. 2008, Phys. Rev. D, 78, 084033, doi: 10.1103/PhysRevD.78.084033
  * Baiotti & Rezzolla (2017) Baiotti, L., & Rezzolla, L. 2017, Rept. Prog. Phys., 80, 096901, doi: 10.1088/1361-6633/aa67bb
  * Bauswein et al. (2019) Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., et al. 2019, Phys. Rev. Lett., 122, 061102, doi: 10.1103/PhysRevLett.122.061102
  * Bauswein et al. (2019) Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., et al. 2019, Physical Review Letters, 122, 061102, doi: 10.1103/PhysRevLett.122.061102
  * Bauswein & Janka (2012) Bauswein, A., & Janka, H.-T. 2012, Phys. Rev. Lett., 108, 011101, doi: 10.1103/PhysRevLett.108.011101
  * Bauswein & Stergioulas (2015) Bauswein, A., & Stergioulas, N. 2015, Phys. Rev. D, 91, 124056, doi: 10.1103/PhysRevD.91.124056
  * Bauswein et al. (2016) Bauswein, A., Stergioulas, N., & Janka, H.-T. 2016, European Physical Journal A, 52, 56, doi: 10.1140/epja/i2016-16056-7
  * Bernuzzi et al. (2014) Bernuzzi, S., Nagar, A., Balmelli, S., Dietrich, T., & Ujevic, M. 2014, Phys. Rev. Lett., 112, 201101, doi: 10.1103/PhysRevLett.112.201101
  * Bishop & Rezzolla (2016) Bishop, N. T., & Rezzolla, L. 2016, Living Reviews in Relativity, 19, 2, doi: 10.1007/s41114-016-0001-9
  * Bose et al. (2018) Bose, S., Chakravarti, K., Rezzolla, L., Sathyaprakash, B. S., & Takami, K. 2018, Phys. Rev. Lett., 120, 031102, doi: 10.1103/PhysRevLett.120.031102
  * Bozzola (2021) Bozzola, G. 2021, The Journal of Open Source Software, 6, 3099, doi: 10.21105/joss.03099
  * Breschi et al. (2022) Breschi, M., Bernuzzi, S., Chakravarti, K., et al. 2022, arXiv e-prints, arXiv:2205.09112, doi: 10.48550/arXiv.2205.09112
  * Breschi et al. (2019) Breschi, M., Bernuzzi, S., Zappa, F., et al. 2019, Phys. Rev. D, 100, 104029, doi: 10.1103/PhysRevD.100.104029
  * Calderon Bustillo et al. (2022) Calderon Bustillo, J., Wong, I. C. F., Sanchis-Gual, N., et al. 2022, arXiv e-prints, arXiv:2205.15029, doi: 10.48550/arXiv.2205.15029
  * Camilletti et al. (2022) Camilletti, A., Chiesa, L., Ricigliano, G., et al. 2022, Mon. Not. R. Astron. Soc., 516, 4760, doi: 10.1093/mnras/stac2333
  * Chabanov et al. (2023) Chabanov, M., Tootle, S. D., Most, E. R., & Rezzolla, L. 2023, Astrophys. J. Lett., 945, L14, doi: 10.3847/2041-8213/acbbc5
  * Chandrasekhar (1969) Chandrasekhar, S. 1969, Ellipsoidal figures of equilibrium
  * Ciolfi et al. (2019) Ciolfi, R., Kastaun, W., Kalinani, J. V., & Giacomazzo, B. 2019, Phys. Rev. D, 100, 023005, doi: 10.1103/PhysRevD.100.023005
  * De Pietri et al. (2019) De Pietri, R., Drago, A., Feo, A., et al. 2019, Astrophys. J., 881, 122, doi: 10.3847/1538-4357/ab2fd0
  * De Pietri et al. (2020) De Pietri, R., Feo, A., Font, J. A., et al. 2020, Phys. Rev. D, 101, 064052, doi: 10.1103/PhysRevD.101.064052
  * Demircik et al. (2022) Demircik, T., Ecker, C., & Järvinen, M. 2022, Phys. Rev. X, 12, 041012, doi: 10.1103/PhysRevX.12.041012
  * Dionysopoulou et al. (2013) Dionysopoulou, K., Alic, D., Palenzuela, C., Rezzolla, L., & Giacomazzo, B. 2013, Phys. Rev. D, 88, 044020, doi: 10.1103/PhysRevD.88.044020
  * East et al. (2015) East, W. E., Paschalidis, V., & Pretorius, F. 2015, Astrophys. J. Lett., 807, L3, doi: 10.1088/2041-8205/807/1/L3
  * East et al. (2019) East, W. E., Paschalidis, V., Pretorius, F., & Tsokaros, A. 2019, Phys. Rev. D, 100, 124042, doi: 10.1103/PhysRevD.100.124042
  * Ellis et al. (2018) Ellis, J., Hütsi, G., Kannike, K., et al. 2018, Phys. Rev. D, 97, 123007, doi: 10.1103/PhysRevD.97.123007
  * Espino et al. (2023) Espino, P. L., Prakash, A., Radice, D., & Logoteta, D. 2023, arXiv e-prints, arXiv:2301.03619, doi: 10.48550/arXiv.2301.03619
  * Foucart et al. (2022) Foucart, F., Duez, M. D., Haas, R., et al. 2022, arXiv e-prints, arXiv:2210.05670, doi: 10.48550/arXiv.2210.05670
  * Franci et al. (2013) Franci, L., De Pietri, R., Dionysopoulou, K., & Rezzolla, L. 2013, Phys. Rev. D, 88, 104028, doi: 10.1103/PhysRevD.88.104028
  * Fujimoto et al. (2023) Fujimoto, Y., Fukushima, K., Hotokezaka, K., & Kyutoku, K. 2023, Phys. Rev. Lett., 130, 091404, doi: 10.1103/PhysRevLett.130.091404
  * Gieg et al. (2019) Gieg, H., Dietrich, T., & Ujevic, M. 2019, arXiv e-prints, arXiv:1908.03135.  https://arxiv.org/abs/1908.03135
  * Gill et al. (2019) Gill, R., Nathanail, A., & Rezzolla, L. 2019, Astrophys. J., 876, 139, doi: 10.3847/1538-4357/ab16da
  * Gonzalez et al. (2022) Gonzalez, A., Zappa, F., Breschi, M., et al. 2022, arXiv e-prints, arXiv:2210.16366, doi: 10.48550/arXiv.2210.16366
  * Grandclement (2010) Grandclement, P. 2010, J. Comput. Phys., 229, 3334, doi: 10.1016/j.jcp.2010.01.005
  * Haas & et al. (2020) Haas, R., & et al. 2020, The Einstein Toolkit, The ”DeWitt-Morette” release, ET_2020_11, Zenodo, doi: 10.5281/zenodo.4298887
  * Hotokezaka et al. (2011) Hotokezaka, K., Kyutoku, K., Okawa, H., Shibata, M., & Kiuchi, K. 2011, Phys. Rev. D, 83, 124008, doi: 10.1103/PhysRevD.83.124008
  * Kiuchi et al. (2015) Kiuchi, K., Cerdá-Durán, P., Kyutoku, K., Sekiguchi, Y., & Shibata, M. 2015, Phys. Rev. D, 92, 124034, doi: 10.1103/PhysRevD.92.124034
  * Kiuchi et al. (2022) Kiuchi, K., Fujibayashi, S., Hayashi, K., et al. 2022, arXiv e-prints, arXiv:2211.07637, doi: 10.48550/arXiv.2211.07637
  * Lehner et al. (2016a) Lehner, L., Liebling, S. L., Palenzuela, C., et al. 2016a, Classical and Quantum Gravity, 33, 184002, doi: 10.1088/0264-9381/33/18/184002
  * Lehner et al. (2016b) Lehner, L., Liebling, S. L., Palenzuela, C., & Motl, P. M. 2016b, Phys. Rev. D, 94, 043003, doi: 10.1103/PhysRevD.94.043003
  * Liebling et al. (2021) Liebling, S. L., Palenzuela, C., & Lehner, L. 2021, Classical and Quantum Gravity, 38, 115007, doi: 10.1088/1361-6382/abf898
  * Liu et al. (2008) Liu, Y. T., Shapiro, S. L., Etienne, Z. B., & Taniguchi, K. 2008, Phys. Rev. D, 78, 024012, doi: 10.1103/PhysRevD.78.024012
  * Löffler et al. (2015) Löffler, F., De Pietri, R., Feo, A., Maione, F., & Franci, L. 2015, Phys. Rev. D, 91, 064057, doi: 10.1103/PhysRevD.91.064057
  * Logoteta, Domenico et al. (2021) Logoteta, Domenico, Perego, Albino, & Bombaci, Ignazio. 2021, A&A, 646, A55, doi: 10.1051/0004-6361/202039457
  * Lucca et al. (2021) Lucca, M., Sagunski, L., Guercilena, F., & Fromm, C. M. 2021, Journal of High Energy Astrophysics, 29, 19, doi: 10.1016/j.jheap.2020.10.002
  * Most et al. (2022) Most, E. R., Haber, A., Harris, S. P., et al. 2022, arXiv e-prints, arXiv:2207.00442.  https://arxiv.org/abs/2207.00442
  * Most et al. (2019) Most, E. R., Papenfort, L. J., Dexheimer, V., et al. 2019, Physical Review Letters, 122, 061101, doi: 10.1103/PhysRevLett.122.061101
  * Most et al. (2019) Most, E. R., Papenfort, L. J., Dexheimer, V., et al. 2019, Phys. Rev. Lett., 122, 061101, doi: 10.1103/PhysRevLett.122.061101
  * Most et al. (2019a) Most, E. R., Papenfort, L. J., & Rezzolla, L. 2019a, Mon. Not. R. Astron. Soc., 490, 3588, doi: 10.1093/mnras/stz2809
  * Most et al. (2019b) Most, E. R., Papenfort, L. J., Tsokaros, A., & Rezzolla, L. 2019b, Astrophys. J., 884, 40, doi: 10.3847/1538-4357/ab3ebb
  * Murguia-Berthier et al. (2021) Murguia-Berthier, A., Ramirez-Ruiz, E., De Colle, F., et al. 2021, Astrophys. J., 908, 152, doi: 10.3847/1538-4357/abd08e
  * Oechslin & Janka (2007) Oechslin, R., & Janka, H.-T. 2007, Phys. Rev. Lett., 99, 121102, doi: 10.1103/PhysRevLett.99.121102
  * Palenzuela et al. (2022) Palenzuela, C., Liebling, S., & Miñano, B. 2022, Phys. Rev. D, 105, 103020, doi: 10.1103/PhysRevD.105.103020
  * Papenfort et al. (2022) Papenfort, L. J., Most, E. R., Tootle, S., & Rezzolla, L. 2022, Mon. Not. Roy. Astron. Soc., 513, 3646, doi: 10.1093/mnras/stac964
  * Papenfort et al. (2021) Papenfort, L. J., Tootle, S. D., Grandclément, P., Most, E. R., & Rezzolla, L. 2021, Phys. Rev. D, 104, 024057, doi: 10.1103/PhysRevD.104.024057
  * Paschalidis et al. (2015) Paschalidis, V., East, W. E., Pretorius, F., & Shapiro, S. L. 2015, Phys. Rev. D, 92, 121502, doi: 10.1103/PhysRevD.92.121502
  * Prakash et al. (2021) Prakash, A., Radice, D., Logoteta, D., et al. 2021, Phys. Rev. D, 104, 083029, doi: 10.1103/PhysRevD.104.083029
  * Radice et al. (2016) Radice, D., Bernuzzi, S., & Ott, C. D. 2016, Phys. Rev. D, 94, 064011, doi: 10.1103/PhysRevD.94.064011
  * Radice et al. (2014a) Radice, D., Rezzolla, L., & Galeazzi, F. 2014a, Mon. Not. R. Astron. Soc. L., 437, L46, doi: 10.1093/mnrasl/slt137
  * Radice et al. (2014b) —. 2014b, Class. Quantum Grav., 31, 075012, doi: 10.1088/0264-9381/31/7/075012
  * Raithel & Most (2022) Raithel, C. A., & Most, E. R. 2022, Astrophys. J. Lett., 933, L39, doi: 10.3847/2041-8213/ac7c75
  * Read et al. (2013) Read, J. S., Baiotti, L., Creighton, J. D. E., et al. 2013, Phys. Rev. D, 88, 044042, doi: 10.1103/PhysRevD.88.044042
  * Rezzolla et al. (2011) Rezzolla, L., Giacomazzo, B., Baiotti, L., et al. 2011, Astrophys. J. Letters, 732, L6, doi: 10.1088/2041-8205/732/1/L6
  * Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, Astrophys. J. Lett., 852, L25, doi: 10.3847/2041-8213/aaa401
  * Rezzolla & Takami (2016) Rezzolla, L., & Takami, K. 2016, Phys. Rev. D, 93, 124051, doi: 10.1103/PhysRevD.93.124051
  * Schnetter et al. (2004) Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Classical and Quantum Gravity, 21, 1465, doi: 10.1088/0264-9381/21/6/014
  * Shibata et al. (2000) Shibata, M., Baumgarte, T. W., & Shapiro, S. L. 2000, Astrophys. J., 542, 453, doi: 10.1086/309525
  * Shibata et al. (2005) Shibata, M., Taniguchi, K., & Uryū, K. 2005, Phys. Rev. D, 71, 084021, doi: 10.1103/PhysRevD.71.084021
  * Sun et al. (2022) Sun, L., Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2022, Phys. Rev. D, 105, 104028, doi: 10.1103/PhysRevD.105.104028
  * Takami et al. (2014) Takami, K., Rezzolla, L., & Baiotti, L. 2014, Phys. Rev. Lett., 113, 091104, doi: 10.1103/PhysRevLett.113.091104
  * Takami et al. (2015) —. 2015, Phys. Rev. D, 91, 064001, doi: 10.1103/PhysRevD.91.064001
  * The CoRe Collaboration (2022) The CoRe Collaboration. 2022, core-watpy, 0.1.1.  http://www.computational-relativity.org/
  * The LIGO Scientific Collaboration & The Virgo Collaboration (2017) The LIGO Scientific Collaboration, & The Virgo Collaboration. 2017, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Tootle et al. (2022) Tootle, S., Ecker, C., Topolski, K., et al. 2022, SciPost Phys., 13, 109, doi: 10.21468/SciPostPhys.13.5.109
  * Tootle et al. (2021) Tootle, S. D., Papenfort, L. J., Most, E. R., & Rezzolla, L. 2021, Astrophys. J. Lett., 922, L19, doi: 10.3847/2041-8213/ac350d
  * Ujevic et al. (2023) Ujevic, M., Gieg, H., Schianchi, F., et al. 2023, Phys. Rev. D, 107, 024025, doi: 10.1103/PhysRevD.107.024025
  * Weih et al. (2020) Weih, L. R., Hanauske, M., & Rezzolla, L. 2020, Phys. Rev. Lett., 124, 171103, doi: 10.1103/PhysRevLett.124.171103
  * Weih et al. (2020) Weih, L. R., Hanauske, M., & Rezzolla, L. 2020, Phys. Rev. Lett., 124, 171103, doi: 10.1103/PhysRevLett.124.171103
  * Wolfram Research (2020) Wolfram Research, I. 2020, Mathematica, Version 12.1.  https://www.wolfram.com/mathematica
  * Zappa et al. (2023) Zappa, F., Bernuzzi, S., Radice, D., & Perego, A. 2023, Mon. Not. R. Astron. Soc., 520, 1481, doi: 10.1093/mnras/stad107

## Supplemental Material

### Toy-model analogy

Although it is possible to characterise the new spectral feature $f_{0}$
simply in terms of the instantaneous GW frequency at the time
$t^{\psi_{4}}_{0}$ at which the Weyl scalar has its first minimum, it is more
interesting to associate with this definition also a physical interpretation.
As mentioned in the main text, $t^{\psi_{4}}_{0}$ effectively corresponds to
when the two stellar cores have reached the minimum separation and are about
to bounce-off each other. At this instant, the radial velocity of the two
stellar cores and the angular velocity of the HMNS has a minimum, so that the
corresponding amplitude of $\psi_{4}$ has a minimum and instantaneous GW
frequency a local maximum.

To illustrate this behaviour, we employ the toy model developed in Takami et
al. (2015) (and also employed in other studies, e.g., Ellis et al. (2018);
Lucca et al. (2021)) to describe the dynamics of the postmerger remnant and
which consists of an axisymmetric disk rotating rapidly at a given angular
frequency, say $\Omega_{2}$, to which two spheres reproducing the two stellar
cores are connected but are also free to oscillate via a spring that connects
them (see Fig. 17 Takami et al. (2015)). In such a system, the two spheres
will either approach each other, decreasing the moment of inertia of the
system, or move away from each other, increasing the moment of inertia.
Because the total angular momentum is essentially conserved, the system’s
angular frequency will vary between a minimum value $\Omega_{1}$
(corresponding to the time when the two spheres are at the largest separation)
and a maximum value $\Omega_{3}$ (corresponding to the time when the two
spheres are at the smallest separation). Overall, the mechanical toy model
will rotate with an angular frequency that is a function of time and bounded
by $\Omega_{1}$ and $\Omega_{3}$, where more time is spent and hence more
spectral power is accumulated. If a damping is introduced, the excursion of
the oscillations between the spheres will decrease over time and eventually
stop; when this happens, the toy model will simply rotate at the frequency
$\Omega_{2}$, as does the HMNS after the transient period and before other
deformations (e.g., the $\ell=2,m=1$ mode) affect it.

The Lagrangian of the toy model can be found in Takami et al. (2014) and from
it it is possible to derive the differential equation for the radial
displacement $r(t)$

$\ddot{r}+\frac{4k(r-r_{0})}{m}-\Bigg{[}\frac{c_{1}}{r^{2}+MR^{2}/(2m)}\Bigg{]}^{2}r+\frac{2br}{m}=0\,,$
(1)

where $M$ and $m$ are the masses of the disc and of the spheres, $R$ is the
radius of the disc, $c_{1}$ an integration constant related to the total
angular momentum, $r_{0}$ the natural displacement of mass and $b$ is
responsible for dissipative effects. The dissipation due to the emission of
GWs and which results in the two spheres getting closer to one another with
time and with a decreasing radial displacement is introduced by varying the
natural displacement $r_{0}$ via via an exponential function
$\exp({-{t}/{\tau}})$ with a suitably chosen relaxation time $\tau$. Note that
this modification carries straight from the Lagrangian to the final ODE
without additional time derivatives, as it is present in the derivative of the
Lagrangian with respect to the canonical coordinate $r(t)$ and not its
conjugate.

The system is sufficiently simple that it is possible to compute the GW
amplitude and frequency derived from the quadrupole formula so that the strain
components read

$h_{+}=\frac{2m}{d}[\\{\dot{r}^{2}+r(\ddot{r}-2r\Omega^{2})\\}\cos(2\varphi)-r(4\dot{r}\Omega+r\dot{\Omega})\sin(2\varphi)~{}]\,,$
$h_{\times}=\frac{2m}{d}[\\{\dot{r}^{2}+r(\ddot{r}-2r\Omega^{2})\\}\sin(2\varphi)+r(4\dot{r}\Omega+r\dot{\Omega})\cos(2\varphi)~{}]\,,$

with $d$ the distance from the source to the detector. The $\psi_{4}$
polarisations are obtained by differentiating twice with respect to time.
Introducing the abbreviations
$A\,\coloneqq\dot{r}^{2}+r(\ddot{r}-2r\Omega^{2})$ and $B\,\coloneqq
r(4\dot{r}\Omega+r\dot{\Omega})$, as well as recalling that
$\Omega(t)=\dot{\varphi}(t)$ we write them now explicitly

$\psi_{4,+}=\frac{2m}{d}[-\sin(2\varphi)\\{4\dot{\varphi}(\dot{A}-B\dot{\varphi})+\ddot{B}+2A\ddot{\varphi}\\}+\cos(2\varphi)\\{-4\dot{\varphi}(\dot{B}+A\dot{\varphi})+\ddot{A}-2B\ddot{\varphi}\\}~{}]\,,$
$\psi_{4,\times}=-\frac{2m}{d}[\cos(2\varphi)\\{4\dot{\varphi}(\dot{A}-B\dot{\varphi})+\ddot{B}+2A\ddot{\varphi}\\}+\sin(2\varphi)\\{-4\dot{\varphi}(\dot{B}+A\dot{\varphi})+\ddot{A}-2B\ddot{\varphi}\\}~{}]\,,$

where $\psi_{4,+}:=\partial^{2}_{t}h_{+}$ and
$\psi_{4,\times}:=\partial^{2}_{t}h_{\times}$.

The amplitude and the instantaneous frequency of the signal are calculated
using the usual definitions

$\displaystyle|\psi_{4}(t)|$ $\displaystyle\coloneqq$
$\displaystyle\sqrt{\psi_{4,\times}^{2}(t)+\psi_{4,+}^{2}(t)}\,,$ (2)
$\displaystyle f_{\rm GW}^{\psi_{4}}(t)$ $\displaystyle\coloneqq$
$\displaystyle\frac{1}{2\pi}\frac{d}{dt}\arctan\left(\frac{\psi_{4,\times}(t)}{\psi_{4,+}(t)}\right)\,.$
(3)

Figure 5: Solution of the toy model for a representative set of parameters.
The left panel reports the GW frequency (top row) from the toy model and from
an actual simulation (marked by the index “NR”), the GW amplitudes (middle
row) and the radial displacement $r(t)$ (bottom row). The right panel shows
instead the relation between the GW quantities and other dynamical quantities
of the toy model, such as the time derivatives of the angular velocity and of
the radial displacement. Note that the GW frequencies and amplitudes are anti-
correlated and the minimum in the radial displacement at $t^{\psi_{4}}_{0}$
corresponds to the maximum frequency $f^{\psi_{4}}_{0}$ and minimum amplitude.
Also, the minima of the amplitude (maxima of the frequency) correspond to
moments in time when the angular frequency $\Omega$ attains a local maximum,
and the radial displacement a local minimum.

The rather lengthy final expressions read:

$\displaystyle\frac{d}{2m}|\psi_{4}(t)|$ $\displaystyle=$
$\displaystyle\big{[}~{}4\sin(2\varphi)\dot{A}\dot{\varphi}+4\cos(2\varphi)\dot{B}\dot{\varphi}+4A\cos(2\varphi)\dot{\varphi}^{2}-4B\sin(2\varphi)\dot{\varphi}^{2}-\cos(2\varphi)\ddot{A}+\sin(2\varphi)\ddot{B}$
(4)
$\displaystyle+2B\cos(2\varphi)\ddot{\varphi}+2A\sin(2\varphi)\ddot{\varphi}~{}\big{]}^{2}+\big{[}~{}4\cos(2\varphi)\dot{A}\dot{\varphi}-4\sin(2\varphi)\dot{B}\dot{\varphi}-4B\cos(2\varphi)\dot{\varphi}^{2}-4A\sin(2\varphi)\dot{\varphi}^{2}$
$\displaystyle+\sin(2\varphi)\ddot{A}+\cos(2\varphi)\ddot{B}+2A\cos(2\varphi)\ddot{\varphi}-2B\sin(2\varphi)\ddot{\varphi}~{}\big{]}^{2},$
$\displaystyle 2\pi\frac{f_{\rm GW}^{\psi_{4}}}{F^{\psi_{4}}}$
$\displaystyle=$ $\displaystyle
48\dot{A}^{2}\dot{\varphi}^{3}+48\dot{B}^{2}\dot{\varphi}^{3}+32A^{2}\dot{\varphi}^{5}+32B^{2}\dot{\varphi}^{5}-32A\dot{\varphi}^{3}\ddot{A}+6\dot{\varphi}\ddot{A}^{2}-32B\dot{\varphi}^{3}\ddot{B}+6\dot{\varphi}\ddot{B}^{2}-24B\dot{\varphi}\ddot{A}\ddot{\varphi}$
(5)
$\displaystyle+24A\dot{\varphi}\ddot{B}\ddot{\varphi}+24A^{2}\dot{\varphi}\ddot{\varphi}^{2}+24B^{2}\dot{\varphi}\ddot{\varphi}^{2}+4B\dot{\varphi}^{2}\dddot{A}-\ddot{B}\dddot{A}-2A\ddot{\varphi}\dddot{A}-4A\dot{\varphi}^{2}\dddot{B}+\ddot{A}\dddot{B}-2B\ddot{\varphi}\dddot{B}$
$\displaystyle+2\big{(}-4A^{2}\dot{\varphi}^{2}+A\ddot{A}+B(-4B\dot{\varphi}^{2}+\ddot{B})\big{)}\dddot{\varphi}+2\dot{B}\big{(}3\ddot{B}\ddot{\varphi}+6\dot{\varphi}^{2}(-3\ddot{A}+4B\ddot{\varphi})-2\dot{\varphi}\dddot{B}+A(40\ddddot{\varphi}+6\ddot{\varphi}^{2}-4\dot{\varphi}\dddot{\varphi})\big{)}$
$\displaystyle+\dot{A}\big{(}6\ddot{A}\ddot{\varphi}+12\dot{\varphi}^{2}(3\ddot{B}+4A\ddot{\varphi})-4\dot{\varphi}\dddot{A}-4B(20\dot{\varphi}^{4}+3\ddot{\varphi}^{2}-2\dot{\varphi}\dddot{\varphi})\big{)},$
$\displaystyle\frac{-1}{F^{\psi_{4}}}$ $\displaystyle:=$ $\displaystyle
16\dot{A}^{2}\dot{\varphi}^{2}+16\dot{B}^{2}\dot{\varphi}^{2}+16A^{2}\dot{\varphi}^{4}+16B^{2}\dot{\varphi}^{4}-8A\dot{\varphi}^{2}\ddot{A}+\ddot{A}^{2}-8B\dot{\varphi}^{2}\ddot{B}+\ddot{B}^{2}+4(-B\ddot{A}+A\ddot{B})\ddot{\varphi}$
(6)
$\displaystyle+4(A^{2}+B^{2})\ddot{\varphi}^{2}+8\dot{A}\dot{\varphi}(-4B\dot{\varphi}^{2}+\ddot{B}+2A\ddot{\varphi})+8\dot{B}\dot{\varphi}(4A\dot{\varphi}^{2}-\ddot{A}+2B\ddot{\varphi})\,.$

Figure 5 reports the solution of the system when considering the following set
of representative parameters: $M=m=10$, $k=0.2$, $b=0.25$, $c_{1}=25$, $R=22$,
$r_{0}=5$, $\dot{r}_{0}=0.01$, $\tau=1000$. More specifically, the left panel
reports the GW frequency (top row) as computed from the toy model and from an
actual simulation (marked by the index “NR”), the corresponding GW amplitudes
(middle row) and the radial displacement $r(t)$ (bottom row). Note how the GW
frequencies and amplitudes are anti-correlated and the minimum in the radial
displacement at $t^{\psi_{4}}_{0}$ corresponds to the maximum frequency, i.e.,
$f^{\psi_{4}}_{0}$, and minimum amplitude, as expected in a condition of quasi
time-symmetry (in the toy model $f^{h}_{0}$ is actually a local minimum of
$f^{h}_{\rm GW}(t)$, a case often found in NR data). The right panel of Fig.
5, on the other hand, is used to show the relation between the GW quantities
of the toy model and other dynamical quantities such as the time derivatives
of the angular velocity, $\partial_{t}\Omega(t)$ and
$\partial^{2}_{t}\Omega(t)$, and of the radial displacement,
$\partial_{t}r(t)$ and $\partial^{2}_{t}r(t)$ (the top row of the right panel
contains the same information as in the top and middle row of the left panel).
It is clear that for the first few oscillations the minima of the amplitude
(maxima of the frequency) correspond to moments in time when the angular
frequency $\Omega$ attains a local maximum, and the radial displacement a
local minimum.

We find it rather remarkable that the model originally developed to provide
justification for the dominant postmerger frequency can be used to mimic the
presence of the $f_{0}$ frequency with no additional adjustments, appearing as
a result of chosen parameters and initial conditions. Finally, while we do not
consider it here, it would be useful to generalise the model to admit unequal
sphere masses, i.e., $m/M<1$, and assess the impact it has on the $f_{0}$
frequency and on the evolution of the $m=2$ and $m=1$ modes.

Figure 6: The same as in Fig. $2$ in the main text but for the GW frequencies
$f^{h}_{\rm mer}$ and $f_{2}$. Note how $\mathcal{F}_{1}$ is equally accurate
in modelling other frequencies, also when obtained from the strain as in the
case for $f^{h}_{\rm mer}$ and $f_{2}$.

### On the effectiveness of the fit

To illustrate effectiveness of the fitting functional form $\mathcal{F}_{1}$
in capturing the spectral properties of the signal in a quasi-universal
manner, we report in Fig. $2$ the same 3D representation as of the GW
frequencies shown in Fig. $2$ but for the GW frequencies $f^{h}_{\rm mer}$ and
$f_{2}$; note that $f^{h}_{\rm mer}$ is referred to as $f_{\rm max}$ in Takami
et al. (2014); Rezzolla & Takami (2016) and that the $f_{2}$ frequency
corresponds to the peak of the PSD computed using the definition in Tootle et
al. (2022), with the time integration encompassing the full GW signal.
Clearly, the fitting function $\mathcal{F}_{1}$ is equally accurate in
modelling other frequencies, quite independently of whether they have been
computed from the Weyl scalar or from the GW strain, as it is as it is the
case for $f^{h}_{\rm mer}$ and $f_{2}$. We also note that the equal-mass limit
of the fitting function $\mathcal{F}_{1}$ for the frequency $f^{h}_{\rm mer}$,
recovers quite accurately the simpler fit examined in Takami et al. (2014)
which expressed the merger frequency as a first-order expression of the tidal
deformability $f^{h}_{\rm mer}=a_{0}+b({\kappa_{2}^{{}_{T}}})^{{1}/{5}}$. When
setting $q=1$ and $n={1}/{5}$ in $\mathcal{F}_{1}$, we obtain
$b=b_{0}+b_{1}+b_{2}=-0.199$ (cf., Table 1), thus resulting in a relative
difference that is $\lesssim 2\%$ with respect to the value $b=-0.195$ found
in Takami et al. (2014). At the same time, we also note that the fit of
$f^{h}_{\rm mer}$ is overall slightly better than that for $f^{\psi_{4}}_{\rm
mer}$ (see Table 1); while we do not believe this to be statistically very
significant, it may be due to the slightly larger range in which
$f^{\psi_{4}}_{\rm mer}$ is measured.

The only exception to the remarkably good fit of the frequencies is found in
the $f^{\psi_{4}}_{0}$ frequency, where a noticeable deviation from the trend
is visible for binary systems with equal or very-unequal masses, regardless of
the EOS employed (see lower part of the bottom-right panel of Fig. $2$). We
attribute this decrease in accuracy to the already described difficulties of
measuring this frequency using the $\psi_{4}$ GW signal, which are even more
severe when using $|h|$. From a physical point of view, the moment of time
symmetry corresponds to the time when the non-axisymmetric quadrupolar
deformations of the HMNS are at a minimum, which leads to a severely
suppressed GW signal (indeed the amplitude is at a minimum). This is
particularly severe for binaries with $q\simeq 1$, where the $\ell=2,m=2$
deformation is the largest and is significantly suppressed at the moment of
time-symmetry, and for binaries with $q\ll 1$, where the $\ell=2,m=2$
deformation is further decreased by the mass asymmetry. Notwithstanding these
difficulties, relative fitting error for $f^{\psi_{4}}_{0}$ [i.e.,
$\max(\Delta f/f$)] is at most $4.2\%$ and only $1.4\%$ on average [i.e.,
$\langle|\Delta f|/f\rangle$].

Freq. | $a_{0}$ | $b_{0}$ | $b_{1}$ | $b_{2}$ | $n$ | $\langle\Delta f/f\rangle$ | $\max(\Delta f/f)$ | $\chi^{2}$ | $\chi^{2}_{\rm red}$ | $R^{2}$
---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  | $[\%]$ | $[\%]$ |  | $\times 10^{-3}$ |
$f^{\psi_{4}}_{\rm mer}$ | $4.589$ | $-0.581$ | $\phantom{-}0.543$ | $-0.236$ | $0.20\phantom{{}^{\dagger}}$ | $0.65$ | $4.52$ | $0.135$ | $1.19$ | $0.934$
$f^{h}_{\rm mer}$ | $4.201$ | $-0.330$ | $\phantom{-}0.198$ | $-0.067$ | $0.20\phantom{{}^{\dagger}}$ | $0.32$ | $2.99$ | $0.035$ | $0.30$ | $0.957$
$f^{\psi_{4}}_{\rm mer}$ | $6.067$ | $-2.142$ | $\phantom{-}0.970$ | $-0.410$ | $0.07^{\dagger}$ | $0.59$ | $4.24$ | $0.114$ | $1.00$ | $0.938$
$f^{h}_{\rm mer}$ | $4.457$ | $-0.578$ | $\phantom{-}0.262$ | $-0.087$ | $0.14^{\dagger}$ | $0.32$ | $2.94$ | $0.034$ | $0.30$ | $0.957$
$f^{\psi_{4}}_{0}$ | $6.550$ | $-2.099$ | $\phantom{-}0.518$ | $-0.304$ | $0.06^{\dagger}$ | $1.38$ | $4.22$ | $0.477$ | $4.92$ | $0.611$
$f_{2}$ | $4.617$ | $-0.170$ | $-0.264$ | $\phantom{-}0.160$ | $0.20^{\dagger}$ | $0.45$ | $1.34$ | $0.030$ | $0.48$ | $0.926$

Table 1: Best-fit values for the coefficients of the functional form
$\mathcal{F}_{1}$. Also reported are the maximal and average relative
difference, as well as the $\chi^{2}$, $\chi^{2}_{\rm{red}}$ and $R^{2}$
coefficients of the fit. Indicated with a $\dagger$ are the best-fit values of
$n$ when this coefficient is constrained by the fit.

### Growth of the SNR ratio

In the main text we have discussed how to use the amplitude ratio
$|\psi_{4}^{21}|/|\psi_{4}^{22}|$ as an effective proxy for the relevance of
the $m=1$ mode deformation. We next demonstrate that a clear correspondence
exists between finding $|\psi_{4}^{21}|/|\psi_{4}^{22}|\sim 1$ and the ratio
of the SNRs in the two modes. We start by recalling that given the power
spectral density of the $m$-th mode of the GW strain decomposition
$\tilde{h}_{m}(f)$ (such as the one presented in Fig. 1 for $m=2$) and a noise
spectral density of the detector $S_{n}(f)$, the corresponding $m$-th SNR is
defined as

${\rm
SNR}_{m=k}:=\left[\int_{0}^{\infty}4\frac{|\tilde{h}_{m=k}(f)|^{2}}{S_{n}(f)}df\right]^{\frac{1}{2}}\,.$
(7)

Clearly, the ratio of the SNRs, ${\rm SNR}_{m=1}/{\rm SNR}_{m=2}$ and the rate
at which it evolves depends on the time $t_{i}$ when the signal starts to be
considered. Because there is no $m=1$ signal during the inspiral, the SNR
ratio would be intrinsically dominated by the $m=2$ component of the signal if
$t_{i}$ was chosen to be the time the signal entered the detector. Hence, to
fairly assess the growing contribution of the $m=1$ mode, we compute a time-
windowed SNR over a running window of width $\Delta T=5\,$ms. We should remark
that this approach is logically and mathematically equivalent to what is done
when computing spectrograms (see, e.g., The LIGO Scientific Collaboration &
The Virgo Collaboration, 2017) and hence determines, at any given time, the
characteristic frequency at which the GW is emitted. In essence, for any time
$\bar{t}$, we compute the SNR as defined in Eq. (7), where the signal in the
time domain is in the interval $t\in[\bar{t}-\Delta T/2,\bar{t}+\Delta T/2]$.
This time-windowed SNR provides an “instantaneous” measure of the ability of a
detector to measure a signal of given strength and is mathematically
equivalent to what is routinely done when computing spectrograms in GW data
analysis.

As demonstrated in Fig. 7, the SNR ratio computed in this way grows to be of
$\mathcal{O}(1)$ at the same time when
$|\psi_{4}^{21}|/|\psi_{4}^{22}|\simeq\mathcal{O}(1)$, thus supporting the
effectiveness of the mode ratio in acting as a proxy for the SNR ratio.
Finally, we note that the estimates provided in Fig. 7 are similar in spirit
to the SNR estimates suggested by Lehner et al. (2016b). We stress that the
results shown should not be interpreted as pointing out to a global-in-time
dominance of the $m=1$ mode; rather, they suggest a an enhanced importance of
this mode that is local-in-time and appears only long past the merger. When
considering the full GW signal, the $m=2$ will always provide the largest
contribution, by far.

Figure 7: Instantaneous SNR ratio as computed for the selected set of
configurations in Fig.$4$ of the main text, where the time-window has been
chosen to be $\Delta T=5\,$ms. The detector’s sensitivity is based on the O3
observing run of LIGO and the angle of observation does not favour either of
the modes. Note that ${\rm SNR}_{m=1}/{\rm SNR}_{m=2}\simeq\mathcal{O}(1)$
around the same time as $|\psi_{4}^{21}|/|\psi_{4}^{22}|\simeq\mathcal{O}(1)$.


# Simultaneous Optimized Orthogonal Matching Pursuit with Application to ECG
Compression

Laura Rebollo-Neira
Department of Applied Mathematics and Data Science
Aston University
B3 7ET, Birmingham, UK

###### Abstract

A greedy pursuit strategy which finds a common basis for approximating a set
of similar signals is proposed. The strategy extends the Optimized Orthogonal
Matching Pursuit approach to selecting the subspace containing the
approximation of all the signals in the set. The method, called Simultaneous
Optimized Orthogonal Matching Pursuit, is stepwise optimal in the sense of
minimizing at each iteration the mean error norm of the joint approximation.
When applied to compression of electrocardiograms, significant gains over
other transformation based compression techniques are demonstrated on the MIT-
BIH Arrhythmia dataset.

Keywords: Simultaneous Optimized Orthogonal Matching Pursuit; Sparse
Representation; ECG compression.

##  1 Introduction

Important signals in everyday life such as natural images, audio, and
electrocardiogram records, are in general highly compressible. This implies
that the original signal, available as a large set of numerical values, can be
transformed into a set of much smaller cardinality or a set containing a large
proportion of zero values. The transformation, which should not compromise the
informational content of the data, is frequently called sparse representation.
Traditional methods for sparse representation of signals are realized by
applying an orthogonal transformation and disregarding the least relevant
points in the transformed domain. Subsequently the signal is recovered by
means of the inverse transformation. However, alternative transformations,
which are not orthogonal but adapted to a signal at hand, have been shown to
render high level of sparsity. Such transformations aim at representing a
signal as a superposition of elements, which are called ‘atoms’ and are
selected from a large set called ‘dictionary’. The superposition is said to be
sparse if it involves a number of atoms much smaller than the number of
numerical values representing the original signal.

Given a dictionary, the problem of finding the sparsest approximation of a
signal, up to some acceptable error, is an NP-hard problem [1]. In practice it
is addressed by tractable methodologies known as Pursuit Strategies. Such
methodologies can be grouped for the most part in two broad categories.
Namely, Basis Pursuit and Greedy Pursuit Strategies. The Basis Pursuit (BP)
approach endeavors to obtain a tractable sparse solution by minimization of
the 1-norm [2]. Greedy algorithms seek for a sparse solution by stepwise
selection of dictionary’s atoms. When dealing with real data the latter are in
general more convenient. From the seminal Matching Pursuit (MP) [3] and
Orthogonal Matching Pursuit (OMP) [4] methods, a number of Greedy Pursuit
Strategies have been developed to improve the process of sparsely representing
single signals [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Due to complexity
issues and memory requirements, most of these techniques are to be applied by
segmenting the signal and approximating each segment independently of the
others. Nonetheless, when the segments bear similarity with each other, for
some applications it is convenient to look for the dictionary’s atoms suitable
to represent all the segments simultaneously. The greedy Pursuit Strategy
which has been dedicated to simultaneously approximate a set of signal is
based on OMP [4] and has been termed Simultaneous Orthogonal Matching Pursuit
(SOMP) [16]. Since in this work we extend the Optimized Orthogonal Matching
Pursuit method [6] to simultaneously approximate a set of signals, we term the
new approach Simultaneous Optimized Orthogonal Matching Pursuit (SOOMP).

The difference between SOMP and the SOOMP approach introduced in this work is
equivalent to the difference between OMP and OOMP methods, both for
approximating single signals. OOMP is stepwise optimal in the sense of
minimizing at each iteration the norm of the residual error. Whilst OMP
minimizes the norm of the error only with respect to the coefficients of the
atomic superposition, OOMP minimizes the norm of the error with respect to
those coefficients and the selection of a new atom. In the case of multiple
signals SOOMP is designed to minimize the mean value error norm. An additional
advantage arises from the proposed implementation. Based on adaptive
biorthogonalization, the SOOMP method produces at each iteration the common
dual basis to the basis of selected atoms. This allows to calculate the
coefficients of the representation of each signal in the set simply by
computation of inner products. The practical relevance of the approach is
illustrated by using it for compression of electrocardiogram (ECG) records.

An ECG signal represents a sequence of heartbeats which, if properly segmented
and aligned, are suitable to be simultaneously approximated. This property is
shown to benefit compression. Reliable comparison with other compression
techniques is made possible by recuse to an adaptive quantization procedure
that facilitates to reconstruct the whole ECG record at the required quality.
The compression results are shown to significantly improve upon results
produced by different transformation based approaches.

The paper is organized as follows: Sec. 2 introduces the problem and the
mathematical notation. Sec. 3 establishes the proposed SOOMP approach for
simultaneous approximation of a set of similar signals. Sec. 4 applies the
method for compressing digital ECG records and produces reliable comparisons
with previously reported results. The conclusions are presented in Sec. 5.

##  2 Mathematical introduction of the problem

In order to pose in mathematical terms the problem to be addressed we need to
introduce the notation used throughout the paper as well as some preliminary
background.

The sets of real, integer, and natural numbers are indicated by
$\mathbb{R},\mathbb{Z}$, and $\mathbb{N}$, respectively. Boldface letters are
used to indicate Euclidean vectors or matrices whilst standard mathematical
fonts indicate components, e.g.,
$\mathbf{f}\in\mathbb{R}^{N},\,N\in\mathbb{N}$ is a vector of components
$f(i),\,i=1,\ldots,N$ and $\mathbf{C}\in\mathbb{R}^{Q\times k}$ is a matrix of
elements $C(i,j),\,i=1,\ldots,Q,\,j=1,\ldots,k$ which when not leaving room
for ambiguity will also be represented as $C(:,j),j=1,\ldots,k$ A set of $Q$
signals of equal length, say $L$, to be simultaneously approximated in a
common subspace, is represented as a set of vectors
$\\{{\mathbf{f}}{\\{q\\}}\in\mathbb{R}^{L},\,q=1,\ldots,Q\\}$. The inner
product is indicated as $\left\langle\cdot,\cdot\right\rangle$, e.g for
$\mathbf{f}\\{1\\}\in\mathbb{R}^{L}$ and $\mathbf{f}\\{2\\}\in\mathbb{R}^{L}$

$\left\langle\mathbf{f}\\{1\\},\mathbf{f}\\{2\\}\right\rangle=\sum_{i=1}^{L}f\\{1\\}(i)f\\{2\\}(i).$

The 2-norm induced by the inner product is denoted as $\|\cdot\|_{2}$, e.g.
for ${\mathbf{f}}{\\{q\\}}\in\mathbb{R}^{L}$

$\|{\mathbf{f}}{\\{q\\}}\|_{2}=\sqrt{\left\langle{\mathbf{f}}{\\{q\\}},{\mathbf{f}}{\\{q\\}}\right\rangle}=\sqrt{\sum_{i=1}^{L}\left(f\\{q\\}(i)\right)^{2}}.$

A set of $M$ vectors

$\mathcal{D}=\left\\{\mathbf{d}_{n}\in\mathbb{R}^{L}\,;\|\mathbf{d}_{n}\|_{2}=1\right\\}_{n=1}^{M},$

such that ${\mbox{\rm{span}}}\left(\mathcal{D}\right)=V_{k}$ and
dim$(V_{k})\leq M$, is called a dictionary for $V_{k}$ and its elements are
called atoms.

In our context a signal $\mathbf{f}$ is considered to be an element of some
inner product space $V$. Each given signal is assumed to be well approximated
by an element, say $\mathbf{f}_{k}$, belonging to a finite dimensional
subspace $V_{k}\subset V$. This assumption implies that, within a tolerance
$\rho$ much larger than the numerical errors in the calculations,
$\mathbf{f}_{k}\in V_{k}$ is accepted to be a good approximation of
$\mathbf{f}\in V$ if $\|\mathbf{f}-\mathbf{f}_{k}\|^{2}_{2}<\rho.$ Examples of
signals fulfilling this definition are, amongst others, images, audio signals,
and electrocardiograms. These are all signals with acceptable approximations
which, without affecting their informational content, do not necessarily
produce a highly accurate point-wise reproduction of the signals. These type
of signals are suitable for lossy compression.

Since this work concerns approximation of similar signals we need to make an
assumption on the signals that will be considered. We say that a finite set of
$Q$ signals $\\{{\mathbf{f}}{\\{q\\}}\in V\\}_{q=1}^{Q}$ are similar if they
can be well approximated in a subspace $V_{k}$ of dimension $k$, with $k$
significantly smaller than the dimention of $V$. This is equivalent to
assuming that there exists a common basis
$\\{{\boldsymbol{\alpha}}_{n}\\}_{n=1}^{k}$ for $V_{k}$ such that each signal
${\mathbf{f}}{\\{q\\}}$ is approximated as

$\tilde{\mathbf{f}}\\{q\\}=\sum_{n=1}^{K}c^{q}_{n}{\boldsymbol{\alpha}}_{n},\quad
q=1,\ldots,Q.$

The quality of the approximated set will be assessed in mean value

${\overline{{\cal{E}}}}=\sum_{q=1}^{Q}p(q)\|{\mathbf{f}}{\\{q\\}}-\tilde{\mathbf{f}}\\{q\\}\|_{2}^{2},$

where $p(q)\geq 0$ with $\sum_{q=1}^{Q}p(q)=1$.

##  3 Strategy for simultaneous approximation of a set of signals

Given a set of similar signals
$\\{{\mathbf{f}}{\\{q\\}}\in\mathbb{R}^{L}\\}_{q=1}^{Q}$ and a dictionary, the
aim is to simultaneously approximate all the signals in the set
$\\{{\mathbf{f}}{\\{q\\}}\in\mathbb{R}^{L}\\}_{q=1}^{Q}$ within a common
subspace
$V_{k}={\mbox{\rm{span}}}\left(\\{\mathbf{d}_{\ell_{n}}\\}_{n=1}^{k}\right)$.
In other words, each signal ${\mathbf{f}}{\\{q\\}}\in\mathbb{R}^{L}$ is to be
approximated as a $k$-term atomic superposition

${\mathbf{f}}{\\{q\\}}^{k}=\sum_{n=1}^{k}c\\{q\\}(n)\mathbf{d}_{\ell_{n}},\quad
q=1,\ldots,Q,$ (1)

where the atoms $\mathbf{d}_{\ell_{n}},\,n=1,\ldots,k$ in (1) are selected
from the given dictionary according to the criterion of optimality that will
be established by Proposition 1 in the next subsection. Let us suppose for the
moment that these atoms are known. Assigning a weight $p(q)\geq 0$ to the
signal ${\mathbf{f}}{\\{q\\}}$, with $\sum_{q=1}^{Q}p(q)=1$, the coefficients
$\mathbf{c}\\{q\\}\in\mathbb{R}^{k}$ in (1) are required to minimize the mean
of the square norm of the errors in the approximation of the set of signals,
i.e.

$\mathbf{c}\\{q\\},\ldots,\mathbf{c}\\{Q\\}=\operatorname*{arg\,min}_{\begin{subarray}{c}\mathbf{c}\end{subarray}^{\prime}\\{q\\},\ldots,\mathbf{c}^{\prime}\\{Q\\}}\sum_{q=1}^{Q}p(q)\|{\mathbf{f}}{\\{q\\}}-\sum_{n=1}^{k}c^{\prime}\\{q\\}(n)\mathbf{d}_{\ell_{n}}\|_{2}^{2}.$
(2)

Since $p(q)\geq 0$ the above minimization is equivalent to finding the
components $c\\{q\\}(n),\,n=1,\ldots,k$ of each vector $\mathbf{c}\\{q\\}$
such that

$c\\{q\\}(1),\ldots,c\\{q\\}(n)=\operatorname*{arg\,min}_{\begin{subarray}{c}c\end{subarray}^{\prime}\\{q\\}(1),\ldots,c^{\prime}\\{q\\}(n)}\|{\mathbf{f}}{\\{q\\}}-\sum_{n=1}^{k}c^{\prime}\\{q\\}(n)\mathbf{d}_{\ell_{n}}\|_{2}^{2}\,\quad
q=1\ldots,Q.$ (3)

Accordingly, the minimization with respect to the coefficients in (1) can be
implemented by adaptive biorthogonalization [17], as proposed within the OOMP
algorithm for a single signal [6],

$c\\{q\\}(n)=\left\langle\boldsymbol{\beta}_{n}^{k},{\mathbf{f}}{\\{q\\}}\right\rangle,\quad
q=1,\ldots,Q,$ (4)

with vectors $\boldsymbol{\beta}_{n}^{k}$ calculated as will be described in
the next section.

The selection of the atoms $\mathbf{d}_{\ell_{n}},\,n=1\ldots,k$ in the
decomposition (1) such that

$\sum_{q=1}^{Q}p(q)\|{\mathbf{f}}{\\{q\\}}-\sum_{n=1}^{k}c\\{q\\}(n)\mathbf{d}_{\ell_{n}}\|_{2}^{2}\,\,\,\text{is
minimized}$

poses an intractable problem (for a dictionary of $M$ atoms there are
$\frac{M!}{(M-k)!k!}$ possibilities to be checked). We address the selection
in a tractable manner by extending the OOMP strategy to simultaneously
approximate a set of similar signals. The extended strategy is refereed to as
SOOMP (Simultaneous OOMP).

###  3.1 SOOMP algorithm

The algorithm is initialized by setting:
${\mathbf{r}}{\\{q\\}}^{0}={\mathbf{f}}{\\{q\\}}$,
${\mathbf{f}}{\\{q\\}}^{0}=0$, $\Gamma=\emptyset$ and $k=0$. The first atom is
selected as the one corresponding to the index $\ell_{1}$ such that

$\ell_{1}=\operatorname*{arg\,max}_{n=1,\ldots,M}\sum_{q=1}^{Q}p(q)\left|\left\langle\mathbf{d}_{n},{\mathbf{r}}{\\{q\\}}^{0}\right\rangle\right|^{2}.$
(5)

This first atom is used to assign
$\mathbf{w}_{1}=\boldsymbol{\beta}_{1}=\mathbf{d}_{\ell_{1}}$, calculate
${\mathbf{r}}{\\{q\\}}^{1}={\mathbf{f}}{\\{q\\}}-\mathbf{d}_{\ell_{1}}\left\langle\mathbf{d}_{\ell_{1}},{\mathbf{f}}{\\{q\\}}\right\rangle$
and iterate as prescribed below.

  * 1)

Upgrade the set $\Gamma\leftarrow\Gamma\cup\ell_{k+1}$, increase $k\leftarrow
k+1$, and select the index of a new atom for the approximation as

$\ell_{k+1}=\operatorname*{arg\,max}_{\begin{subarray}{c}n=1,\ldots,M\\\
n\notin\Gamma\end{subarray}}\sum_{q=1}^{Q}p(q)\frac{\left|\left\langle\mathbf{d}_{n},{\mathbf{r}}{\\{q\\}}^{k}\right\rangle\right|^{2}}{1-\sum_{i=1}^{k}|\left\langle\mathbf{d}_{n},\mathbf{\widetilde{w}}_{i}\right\rangle|^{2}},\quad\text{with}\quad\mathbf{\widetilde{w}}_{i}=\frac{\mathbf{w}_{i}}{\|\mathbf{w}_{i}\|}.$
(6)

  * 2)

Compute the corresponding new vector $\mathbf{w}_{k+1}$ as

$\begin{split}\mathbf{w}_{k+1}=\mathbf{d}_{\ell_{k+1}}-\sum_{i=1}^{k}\frac{\mathbf{w}_{i}}{\|\mathbf{w}_{i}\|^{2}}\left\langle\mathbf{w}_{i},\mathbf{d}_{\ell_{k+1}}\right\rangle,\end{split}$
(7)

including for numerical accuracy the re-orthogonalization step:

$\mathbf{w}_{k+1}\leftarrow\mathbf{w}_{k+1}-\sum_{i=1}^{k}\frac{\mathbf{w}_{i}}{\|\mathbf{w}_{i}\|^{2}}\left\langle\mathbf{w}_{i},\mathbf{w}_{k+1}\right\rangle.$
(8)

  * 3)

Upgrade vectors $\boldsymbol{\beta}_{n}^{k}$ as

$\boldsymbol{\beta}_{k+1}^{k+1}=\frac{\mathbf{w}_{k+1}}{\|\mathbf{w}_{k+1}\|^{2}},\quad\boldsymbol{\beta}_{n}^{k+1}=\boldsymbol{\beta}_{n}^{k}-\boldsymbol{\beta}_{k+1}^{k+1}\left\langle\mathbf{d}_{\ell_{k+1}},\boldsymbol{\beta}_{n}^{k}\right\rangle,\quad
n=1,\ldots,k.$ (9)

  * 4)

Update ${\mathbf{r}}{\\{q\\}}^{k}$ as

$\displaystyle{\mathbf{r}}{\\{q\\}}^{k+1}$ $\displaystyle=$
$\displaystyle{\mathbf{r}}{\\{q\\}}^{k}-\left\langle\mathbf{w}_{k+1},{\mathbf{f}}{\\{q\\}}\right\rangle\frac{\mathbf{w}_{k+1}}{\|\mathbf{w}_{k+1}\|^{2}}.$
(10)

  * 5)

If for a given value $\rho$ the condition
$\sum_{q=1}^{Q}p(q)\|{\mathbf{r}}{\\{q\\}}^{k+1}\|_{2}^{2}<\rho$ has been met
stop the selection process. Otherwise repeat steps 1) - 5).

Once the iterations have finished calculate the coefficients for the
decomposition (1) as

$\mathbf{c}\\{q\\}_{n}=\left\langle\boldsymbol{\beta}_{n}^{k},{\mathbf{f}}{\\{q\\}}\right\rangle,\quad
n=1,\ldots,k,\,q=1,\ldots,Q.$

For $q=1,\ldots,Q$ calculate the final approximation of each signal
${\mathbf{f}}{\\{q\\}}$ as

${\mathbf{f}}{\\{q\\}}^{k}={\mathbf{f}}{\\{q\\}}-{\mathbf{r}}{\\{q\\}}^{k}.$

######  Remark 1.

The set of vectors $\boldsymbol{\beta}_{n}^{k},\,n=1,\ldots,k,\,q=1,\ldots,Q$
as given in (9) fulfills that

${\mathbf{f}}{\\{q\\}}^{k}=\sum_{n=1}^{k}\left\langle\beta_{n}^{k},{\mathbf{f}}{\\{q\\}}\right\rangle\mathbf{d}_{\ell_{n}}=\hat{P}_{V_{k}}{\mathbf{f}}{\\{q\\}},\quad
q=1,\ldots,Q,$

where $\hat{P}_{V_{k}}{\mathbf{f}}{\\{q\\}}$ is the orthogonal projector of
${\mathbf{f}}{\\{q\\}}$ onto
$V_{k}={\mbox{\rm{span}}}\\{\mathbf{d}_{\ell_{n}}\\}_{n=1}^{k}$. Please find
the proof in [6], or as a particular case of the more general proof in [17].

######  Proposition 1.

The recursive selection of the indices $\ell_{1},\ldots,\ell_{k}$, as proposed
in (6), is stepwise optimal. It minimizes, at each iteration, the mean of the
square distance between the set of signals
${\mathbf{f}}{\\{q\\}},\,q=1,\ldots,Q$ and their corresponding approximations
${\mathbf{f}}{\\{q\\}}^{k},\,q=1,\ldots,Q$.

###### Proof.

For $k=0$ it is clear that $\ell_{1}$ selected as in (5) minimizes the mean of
the square distance $\overline{\mathcal{E}^{1}}$ as given by

$\overline{\mathcal{E}^{1}}=\sum_{q=1}^{Q}p(q)\|{\mathbf{f}}{\\{q\\}}-{\mathbf{f}}{\\{q\\}}^{1}\|_{2}^{2}=\sum_{q=1}^{Q}p(q)(\|{\mathbf{f}}{\\{q\\}}\|_{2}^{2}-|\left\langle\mathbf{d}_{\ell_{1}},{\mathbf{f}}{\\{q\\}}\right\rangle|^{2}).$

Let us assume that the indices $\ell_{1},\ldots,\ell_{k}$ selected as proposed
in (6) minimize, in the specified stepwise sense, the mean square distance

$\overline{\mathcal{E}^{k}}=\sum_{q=1}^{Q}p(q)\|{\mathbf{f}}{\\{q\\}}-{\mathbf{f}}{\\{q\\}}^{k}\|_{2}^{2}.$

We shall prove by induction that if the atoms
$\mathbf{d}_{\ell_{1}},\ldots,\mathbf{d}_{\ell_{k}}$ are fixed, at iteration
$k+1$ the atom $\mathbf{d}_{\ell_{k+1}}$ selected as in (6) minimizes
$\overline{\mathcal{E}^{k+1}}$. The proof stems from the fact that at
iteration $k$ the approximation ${\mathbf{f}}{\\{q\\}}^{k}$ of each signal
${\mathbf{f}}{\\{q\\}}$ is the orthogonal projection of
${\mathbf{f}}{\\{q\\}}$ onto the subspace
$V_{k}={\mbox{\rm{span}}}\\{\mathbf{d}_{\ell_{n}}\\}_{n=1}^{k}$ (c.f. Remark
1).

Consider that $V_{k}$ is augmented by one element, say
$\mathbf{d}_{\ell_{k+1}}\notin V_{k}$, so that
$V_{k+1}=V_{k}\oplus\mathbf{d}_{\ell_{k+1}}$, where $\oplus$ indicates direct
sum. The orthogonal projection of each signal
${\mathbf{f}}{\\{q\\}},\,q=1,\ldots,Q$ onto $V_{k+1}$ can be expressed as

${\mathbf{f}}{\\{q\\}}^{k+1}=\hat{P}_{V_{k+1}}{\mathbf{f}}{\\{q\\}}=\hat{P}_{V_{k}}{\mathbf{f}}{\\{q\\}}+\frac{\mathbf{w}_{k+1}}{\|\mathbf{w}_{k+1}\|_{2}^{2}}\left\langle\mathbf{w}_{k+1},{\mathbf{f}}{\\{q\\}}\right\rangle\quad\text{with}\quad\mathbf{w}_{k+1}=\mathbf{d}_{\ell_{k+1}}-\hat{P}_{V_{k}}\mathbf{d}_{\ell_{k+1}}.$

Thus

$\displaystyle\|{\mathbf{f}}{\\{q\\}}-{\mathbf{f}}{\\{q\\}}^{k+1}\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\|{\mathbf{f}}{\\{q\\}}-\hat{P}_{V_{k+1}}{\mathbf{f}}{\\{q\\}}\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\|{\mathbf{f}}{\\{q\\}}-\hat{P}_{V_{k}}{\mathbf{f}}{\\{q\\}}-\frac{\mathbf{w}_{k+1}}{\|\mathbf{w}_{k+1}\|_{2}^{2}}\left\langle\mathbf{w}_{k+1},{\mathbf{f}}{\\{q\\}}\right\rangle\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\|{\mathbf{f}}{\\{q\\}}-\hat{P}_{V_{k}}{\mathbf{f}}{\\{q\\}}\|_{2}^{2}-\frac{|\left\langle\mathbf{w}_{k+1},{\mathbf{f}}{\\{q\\}}\right\rangle|_{2}^{2}}{\|\mathbf{w}_{k+1}\|_{2}^{2}}.$

Since $\|{\mathbf{f}}{\\{q\\}}-\hat{P}_{V_{k}}{\mathbf{f}}{\\{q\\}}\|_{2}^{2}$
is optimized and fixed at iteration $k$, it is true that at iteration $k+1$
the index of the atom which minimizes $\overline{\mathcal{E}^{k+1}}$ fulfils

$\displaystyle\ell_{k+1}$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,max}_{\begin{subarray}{c}n=1,\ldots,M\\\
n\notin\Gamma\end{subarray}}\sum_{q=1}^{Q}p(q)\frac{|\left\langle\mathbf{w}_{k+1},{\mathbf{f}}{\\{q\\}}\right\rangle|_{2}^{2}}{\|\mathbf{w}_{k+1}\|_{2}^{2}}$
(11) $\displaystyle=$
$\displaystyle\operatorname*{arg\,max}_{\begin{subarray}{c}n=1,\ldots,M\\\
n\notin\Gamma\end{subarray}}\sum_{q=1}^{Q}p(q)\frac{\left|\left\langle\mathbf{d}_{\ell_{n}}-\hat{P}_{V_{k}}\mathbf{d}_{\ell_{n}},{\mathbf{f}}{\\{q\\}}\right\rangle\right|^{2}}{\|\mathbf{d}_{\ell_{n}}-\hat{P}_{V_{k}}\mathbf{d}_{\ell_{n}}\|_{2}^{2}}.$

The proof is concluded using the self-adjoint properties of $\hat{P}_{V_{k}}$
to write:

$\displaystyle\left|\left\langle\mathbf{d}_{\ell_{n}}-\hat{P}_{V_{k}}\mathbf{d}_{\ell_{n}},{\mathbf{f}}{\\{q\\}}\right\rangle\right|$
$\displaystyle=$
$\displaystyle\left|\left\langle\mathbf{d}_{\ell_{n}},{\mathbf{f}}{\\{q\\}}\right\rangle-\left\langle\hat{P}_{V_{k}}\mathbf{d}_{\ell_{n}},{\mathbf{f}}{\\{q\\}}\right\rangle\right|$
$\displaystyle=$
$\displaystyle\left|\left\langle\mathbf{d}_{\ell_{n}},{\mathbf{f}}{\\{q\\}}-\hat{P}_{V_{k}}{\mathbf{f}}{\\{q\\}}\right\rangle\right|$
$\displaystyle=$
$\displaystyle\left|\left\langle\mathbf{d}_{\ell_{n}},{\mathbf{r}}{\\{q\\}}^{k}\right\rangle\right|.$

Moreover, since all atoms are normalized and the set
$\\{\mathbf{\widetilde{w}}_{i}\\}_{i=1}^{k}$ is an orthonormal basis for
$V_{k}$ we have

$\hat{P}_{V_{k}}\mathbf{d}_{\ell_{n}}=\sum_{i=1}^{k}\mathbf{\widetilde{w}}_{i}\left\langle\mathbf{\widetilde{w}}_{i},\mathbf{d}_{\ell_{n}}\right\rangle,$

so that

$\|\mathbf{d}_{\ell_{n}}-\hat{P}_{V_{k}}\mathbf{d}_{\ell_{n}}\|_{2}^{2}=1-\|\hat{P}_{V_{k}}\mathbf{d}_{\ell_{n}}\|_{2}^{2}=1-\sum_{i=1}^{k}|\left\langle\mathbf{\widetilde{w}}_{i},\mathbf{d}_{\ell_{n}}\right\rangle|^{2},$

which shows the equivalence between (11) and (6). ∎

##  4 Application to compression of ECG records

A digital ECG signal represents a sequence of heartbeats. In a typical record
each heartbeat is characterized by a combination of three graphical
deflections, known as QRS complex, and two lateral and less visually
noticeable P and T waves. A short segment of a typical ECG record is
illustrated in Fig.1.

Figure 1: A short segment of an ECG record

In order to simultaneously approximate all the beats in a record we need to
segment and align the beats to meet the requirement of being similar. The
procedure is discussed in the next subsection.

###  4.1 Segmentation and alignment of heartbeats

The QRS complex is segmented once the central R peak is detected. This can be
effectively done by the Pan Tompkins method [18]. In our numerical examples we
use the off-the-shelf MATLAB implementation of this algorithm [19]. Since the
distance between peaks in a record is not uniform, the length of the segmented
beats should be passed to the decoder. The segmented peaks are placed in
arrays ${\mathbf{f}}{\\{q\\}},\,q=1,\ldots,Q$ of equal length $L$ by padding
with zeros. Fig.4 illustrates the resulting configuration with 80 heartbeats.

Figure 2: Configuration resulting by segmentation and alignment of 80
heartbeats for illustration purposes.

The segmented and aligned heartbeats are simultaneously approximated using the
SOOMP approach by assigning the same weight to each heartbeat, i.e. setting
$p(q)=\frac{1}{Q},\,q=1,\ldots,Q$. The adopted dictionary is the Cohen-
Daubechies-Feauveau $\mathbf{CDF97}$ dictionary, of redundancy approximately
two, introduced in [20, 21]. Given a partition $x_{i},\,i=1,\ldots,N$ of the
interval $[c,d]$ the dictionary is constructed as follows [20, 21].

$\mathcal{D}={\cal V}_{0}\cup{\cal W}_{0}\cup{\cal W}_{1}\cup{\cal
W}_{2}\cup{\cal W}_{3}\cup{\cal W}_{4},$ (12)

with

${\cal
V}_{0}=\\{\phi(x_{i}-\frac{k}{2})|_{[c,d]},\,k\in\mathbb{Z},\,i=1,\ldots,N\\},$
(13)

and

${\cal
W}_{j}=\\{2^{j/2}\psi(2^{j}x_{i}-\frac{k}{2})|_{[c,d]},\,k\in\mathbb{Z},\,i=1,\ldots,N\\},$
(14)

where $\psi(2^{j}x_{i}-\frac{k}{2})|_{[c,d]}$ indicates the restriction of the
function $\psi(2^{j}x_{i}-\frac{k}{2})$ to the interval $[c,d]$. The prototype
functions $\phi(x)$ and $\psi(x)$ are plotted in the left and right graphs of
Fig.3 respectively. The MATLAB codes for producing numerically both functions
and building the dictionary (12) are described in [21]. The codes have been
made available in [23] together with of the complete MATLAB software for
reproducing the numerical examples in this work.

Figure 3: Cohen-Daubechies-Feauveau scaling and wavelet functions [22].

The left graph of Fig. 4 is the two dimensional image of the segmented and
aligned heartbeats corresponding to the record 100 in the MIT-BIH Arrhythmia
database [24].

The coefficients in the decompositions (1) are placed in a two dimensional
array $\mathbf{C}\in\mathbb{R}^{Q\times k}$.

Figure 4: The graph on the left is the image of the magnitude of the aligned
heartbeats in record 100. The graph on the right is the image of the magnitude
of the coefficients in the simultaneous approximation of the heartbeats.

The right graph in Fig.4 shows the magnitude of the $\mathbf{C}$ entries
corresponding to the approximation of the heartbeats represented in the left
graph of the same figure. We notice that the size of the array in the left
graph ($2273\times 426$) is reduced to a smaller array of size $2273\times 19$
containing the coefficients in (1) with respect to the common basis
$\mathbf{d}_{\ell_{1}},\ldots\mathbf{d}_{\ell_{k}}$. As it is clear from the
location of the brightest pixels in the right image of Fig.4 the coefficients
of largest magnitude are all concentrated in the first vertical lines. This
implies that to favor compression it is convenient to apply an orthogonal
transformation to map the coefficients in the vertical direction to smaller
values, which eventually will be quantized to zero. Thus, by applying the
discrete cosine transform on each columns of $\mathbf{C}$, we create the
transformed array $\mathbf{B}\in\mathbb{R}^{Q\times k}$ with the following
entries

$\mathbf{B}(:,n)=\widehat{\mathrm{dct}}\mathbf{C}(:,n),\quad n=1,\ldots,k,$
(15)

where $\widehat{\mathrm{dct}}\mathbf{C}(:,n)$ indicates the one dimensional
discrete cosine transform operating on the $n$-th column of array
$\mathbf{C}$. The magnitude of the transformed points is shown Fig. 5.

Figure 5: Magnitude of the array $\mathbf{B}$ (c.f. (15)).

###  4.2 Encoding

At the encoding step the $Q\times k$ array $\mathbf{B}$ is expressed as a
vector $\mathbf{b}=(b(1),\ldots,b(K))$ of $K=Q\cdot k$ components, adopting
the column-major order. The encoding of this vector follows the procedure
outline in [25]. The components of $\mathbf{b}$ are converted to integer
numbers by a mid-tread uniform quantizer as follows:

$b^{\Delta}(i)=\left\lfloor\frac{b(i)}{\Delta}+\frac{1}{2}\right\rfloor,\quad
i=1,\ldots,K,$ (16)

where $\lfloor x\rfloor$ indicates the largest integer smaller or equal to $x$
and $\Delta$ is the quantization parameter. For comparison with results in
other publications in the numerical examples the quantization parameter
$\Delta$ is set to produce the required quality of the reconstructed signal.

The absolute value of the elements (16) are placed in a smaller vector, say
$\mathbf{b^{\prime}}=(b^{\prime}(1),\ldots,b^{\prime}(K^{\prime}))$, after the
elimination of zeros. The signs are encoded separately in a vector
$\mathbf{s}=(s(1),\ldots,s(K^{\prime}))$ using a binary alphabet: 1 for $+$
and 0 for $-$.

Assuming that the nonzero values in (16) occur at the positions
$\jmath_{i},\ldots,\jmath_{K^{\prime}}$, these indices are re-ordered in
ascending order
$\jmath_{i}\rightarrow\tilde{\jmath}_{i},\,i=1,\ldots,K^{\prime}$, i.e.
$\tilde{\jmath}_{i}<\tilde{\jmath}_{i+1},\,i=1,\ldots,K^{\prime}$. This
induces new order in the coefficients,
$\mathbf{b^{\prime}}\rightarrow\mathbf{\tilde{b}^{\prime}}$ and in the
corresponding signs $\mathbf{s}\rightarrow\mathbf{\tilde{s}}$. Defining
$\delta(i)=\tilde{\jmath}_{i}-\tilde{\jmath}_{i-1},\,i=2,\ldots,K^{\prime}$
the array
$\boldsymbol{{\delta}}=(\tilde{\jmath}_{1},\delta(2),\ldots,\delta(K^{\prime}))$
stores the indices $\tilde{\jmath}_{1},\ldots,\tilde{\jmath}_{K^{\prime}}$
with unique recovery.

Finally the vectors
$\mathbf{\tilde{b}^{\prime}},\mathbf{\tilde{s}},\boldsymbol{{\delta}}$, as
well as the length of the heartbeats $\mathbf{h}$, are compressed using
adaptive Huffman coding implemented by the off-the-shelf MATLAB function
Huff06 [26]. The additional numbers which have to be passed to the decoder
are:
(i) The indices $\ell_{i},\,i=1,\ldots,k$ of the selected dictionary’s atoms
forming the common basis.
(ii) The quantization parameter $\Delta$.
(iii) The mean value of the 1D ECG record (if not previouly substracted).
(iv) The number of rows and columns of $\mathbf{C}$, i.e. $Q$ and $k$.

###  4.3 1D ECG signal recovery

At the decoding stage, after reverting Huffman coding, the locations
$\tilde{\jmath}_{1},\ldots,\tilde{\jmath}_{K}$ of the nonzero entries in the
transformed array after quantization are readily obtained. This allows the
recovery of the array $\mathbf{B}^{r}$ as follows.

  * (i)

Set $b^{r}(i)=0,\,i=1,\ldots,K$ and
$b^{r}(\tilde{\jmath}_{i})=(2\tilde{s}(i)-1)\tilde{b}^{\prime}(i)\Delta,\,i=1,\ldots,K$.

  * (ii)

Reshape the vector $\mathbf{b}^{r}$ to produce a 2D array $\mathbf{B}^{r}$ of
size $Q\times k$. The array $\mathbf{C}^{r}$ is recovered from the
$\mathbf{B}^{r}$ one by inverting the $\widehat{\mathrm{dct}}$ transformation
(c.f. (15)).

  * (iii)

Each row of the recovered array $\mathbf{C}^{r}$ gives the coefficients in the
decomposition (1) of the approximated heartbeats, i.e.
$\mathbf{f}^{\rm{r}}\\{q\\}=\sum_{i=1}^{k}\mathbf{C}^{r}(q,i)\mathbf{d}_{\ell_{i}},\quad
q=1,\ldots,Q.$

  * (iv)

Finally the reconstructed beats $\mathbf{f}^{\rm{r}}\\{q\\}$ are assembled in
a 1D record using the distance between heartbeats that was stored in the
vector $\mathbf{h}$.

The achieved compression ratio $\mathrm{CR}$, which is defined as

$\text{CR}=\frac{{\text{Size of the uncompressed file}}}{{\text{Size of the
compressed file}}},$ (17)

depends on the required quality of the recovered signal. In the numerical
examples the quality of the recovered records is assessed by the
$\mathrm{PRDN}$ metric defined as follows

$\mathrm{PRDN}=\frac{\|\mathbf{f}-\mathbf{f}^{\rm{r}}\|}{\|\mathbf{f}-\overline{\mathbf{f}}\|}\times
100\%$ (18)

where $\mathbf{f}$ is the whole ECG record, $\mathbf{f}^{\rm{r}}$ is the
reconstructed record from the compressed file and ${\overline{\mathbf{f}}}$ is
the mean of $\mathbf{f}$. It is pertinent to stress the importance of adopting
the normalized metric (18) for comparison of reconstruction quality. The
subtraction of $\overline{\mathbf{f}}$ avoids dependence on the signal
baseline.

##  5 Numerical Tests

For the numerical test we use the MIT-BIH Arrhythmia database [24]. Each of
the records is of 30 min length, consisting of $N=650000$ 11-bit samples at a
frequency of 360 Hz.

For comparison purposes we compress the subset of records reported in [27],
[28], and [29] and reproduce the values of $\mathrm{PRDN}$ in those
publications. This is achieved as follows: The SOOMP method is applied to
approximate the set of heartbeats in each record up to 80% the target
$\mathrm{PRDN}$. The quantization parameter $\Delta$ is automatically fixed,
by a bisection algorithm, in order to reproduce the target $\mathrm{PRDN}$ for
the whole record within two decimal places.

The first, second and third columns of Table 1 reproduce the results published
in [27]. The comparison is relevant because the approach [27] is also based on
approximation of heartbeats using a dictionary. The techniques are very
different though. Whilst our dictionary does not have to be stored because it
is numerically generated, the dictionary in [27] is part of the ECG record to
be compressed. Moreover, the method for finding the sparse representation is
different and so is the procedure to store the parameters that should be
passed to the decoder.

Table 1: Comparison with the results in [27]. The first collumn lists the records considered in [27]. The second column displays the values of $\mathrm{PRDN}$ and the third collum their $\mathrm{CRs}$. Our $\mathrm{CRs}$ for the same $\mathrm{PRDN}$ are shown in the forth column. The fifth column shows the $\mathrm{CRs}$ obtained with the fast approach [25]. Record | $\mathrm{PRDN}$ |  $\mathrm{CR}$[27] |  $\mathrm{CR}$ prop. |  $\mathrm{CR}$ [25]
---|---|---|---|---
100 | 18.03 | 78.20 | 143.99 | 36.51
100 | 17.22 | 75.12 | 139.47 | 35.25
101 | 14.66 | 80.24 | 102.58 | 31.26
101 | 12.91 | 76.46 | 82.31 | 30.31
102 | 18.54 | 58.54 | 58.49 | 33.89
102 | 18.16 | 48.47 | 58.13 | 33.48
103 | 12.57 | 46.32 | 90.91 | 30.84
103 | 11.57 | 44.33 | 86.27 | 29.61
109 | 13.70 | 24.86 | 145.80 | 51.23
109 | 9.97 | 23.53 | 97.73 | 36.91
111 | 26.20 | 31.05 | 121.09 | 38.29
111 | 19.51 | 29.44 | 60.49 | 32.20
112 | 16.58 | 34.06 | 91.48 | 35.05
112 | 15.99 | 35.49 | 85.59 | 34.32
113 | 14.08 | 37.42 | 90.76 | 32.49
113 | 9.82 | 32.55 | 55.30 | 27.68
115 | 9.76 | 38.26 | 62.31 | 24.52
115 | 9.18 | 36.57 | 57.32 | 23.74
117 | 14.42 | 38.94 | 120.89 | 36.94
117 | 13.38 | 37.13 | 105.97 | 35.74
119 | 32.19 | 16.26 | 153.33 | 90.40
119 | 16.36 | 15.24 | 78.81 | 48.08
121 | 17.36 | 26.67 | 111.74 | 46.45
121 | 15.63 | 25.29 | 100.72 | 41.11
Average | 17.33 | 41.9 | 107.78 | 40.65
Average | 14.14 | 39.97 | 84.00 | 34.04

Our compression results are shown in the forth column of Table 1. These
results demonstrate a significant gain in $\mathrm{CR}$ for the same recovery
quality. For further comparison we apply the fast compression algorithm [25],
which does not require peak segmentation or Huffman coding. This method has
been already shown to improve the average $\mathrm{CR}$ for the 48 records in
the MIT-BIH Arrhythmia dataset with respect to the results in [30], [31], and
[32], for a broad rage of average qualities. For comparison with [27] in Table
1 the compression is realized to reproduce the $\mathrm{PRDN}$ listed in the
second column for each record.

Table 2: Same description as in Table 1 but the comparison is with the results of Table I in [28]. Record | $\mathrm{PRDN}$ |  $\mathrm{CR}$[28] |  $\mathrm{CR}$ prop. |  $\mathrm{CR}$ [25]
---|---|---|---|---
100 | 11.46 | 39.81 | 64.47 | 23.03
101 | 14.13 | 42.04 | 95.53 | 30.91
102 | 19.94 | 41.09 | 63.69 | 35.09
103 | 6.72 | 41.24 | 39.17 | 21.05
107 | 13.27 | 41.84 | 71.40 | 38.15
109 | 7.31 | 38.25 | 59.34 | 28.46
111 | 13.94 | 41.73 | 41.59 | 25.14
115 | 8.04 | 42.71 | 47.79 | 22.04
117 | 10.00 | 46.75 | 51.56 | 26.17
118 | 15.33 | 39.60 | 66.44 | 28.26
119 | 9.67 | 41.97 | 43.98 | 30.71
213 | 13.63 | 32.58 | 62.24 | 25.72
222 | 22.44 | 40.69 | 31.95 | 24.77
232 | 20.76 | 42.28 | 59.66 | 35.52
Average | 13.33 | 40.90 | 57.09 | 28.22

The first, second and third columns of Table 2 reproduce the results published
in [28], which are achieved with an approach based on the Singular Value
Decomposition (SVD).

Our compression ratios ($\mathrm{CRs}$) are shown in the the forth column of
Table 2. The fifth column shows the $\mathrm{CRs}$ produced by the fast
compression algorithm [25].

Table 3: Same description as in Table 1 but the comparison is with the results of Table I in [29]. Record | $\mathrm{PRDN}$ |  $\mathrm{CR}$[29] |  $\mathrm{CR}$ prop. |  $\mathrm{CR}$ [25]
---|---|---|---|---
100 | 11.55 | 50.70 | 65.47 | 23.14
101 | 11.29 | 58.54 | 64.13 | 27.86
103 | 9.16 | 54.16 | 61.37 | 25.70
107 | 14.53 | 53.12 | 78.31 | 40.34
109 | 11.83 | 46.43 | 120.08 | 44.67
111 | 16.40 | 53.39 | 50.11 | 28.34
115 | 8.94 | 56.77 | 58.76 | 23.49
117 | 12.43 | 66.15 | 90.39 | 34.62
119 | 10.28 | 56.03 | 46.99 | 32.17
214 | 17.03 | 50.84 | 87.66 | 54.04
223 | 17.49 | 45.38 | 86.06 | 40.80
Avegave | 12.81 | 53.77 | 73.57 | 34.11

The first, second and third columns of Table 3 reproduce the results published
in [29], which are also obtained with a Singular Value Decomposition based
approach. Our $\mathrm{CRs}$ are shown in the forth column of this table. The
fifth column shows the $\mathrm{CRs}$ produced by the fast compression
algorithm [25].

Note: The MATLAB software for reproducing the tables is available on
http://www.nonlinear-approx.info/examples/node017.html

##  6 Conclusions

The Optimized Orthogonal Matching Pursuit approach has been extended with the
purpose of selecting a common basis for the simultaneous approximation of a
set of similar signals. The extended approach, termed Simultaneous Optimized
Orthogonal Matching Pursuit, minimizes at each iteration the mean error norm
of the joint approximation. The applicability of the method was illustrated by
the simultaneous approximation of heartbeats in ECG records taken from the
MIT-BIH Arrhythmia database. The particular records were selected for
comparison purposes as in [27], [28], and [29]. It was demonstrated that
simultaneous approximations of heartbeats can be used for compressing a whole
record. The adopted compression strategy was shown to improve upon compression
results achieved by other methods for the same reconstruction quality. The
comparison was made possible by means of an iterative quantization procedure
which delivers the required quality. We feel confident that the proposed
technique could be of assistance to other signal processing applications which
benefit from the use of a common basis for the approximation of a set of
similar signals.

## References

  * [1] B. K. Natarajan, “Sparse Approximate Solutions to Linear Systems”, SIAM Journal on Computing, 24, 227–234 (1995).
  * [2] S. S. Chen, D. L. Donoho, and M. A Saunders, “Atomic Decomposition by Basis Pursuit”, SIAM Journal on Scientific Computing, 20, 33–61 (1998).
  * [3] S. G. Mallat and Z. Zhang, “Matching Pursuits with Time-Frequency Dictionaries”, IEEE Trans. Signal Process., 41, 3397–3415 (1993).
  * [4] Y.C. Pati, R. Rezaiifar, and P.S. Krishnaprasad,  “Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition,”  Proc. of the 27th ACSSC,1, 40–44 (1993).
  * [5] R. Gribonval, “Fast Matching Pursuit with a multiscale dictionary of Gaussian Chirps”, IEEE Trans. Signal Process., 49, 994–1001 (2001).
  * [6] L. Rebollo-Neira and D. Lowe, “Optimized orthogonal matching pursuit approach”, IEEE Signal Process. Letters, 9, 137–140 (2002).
  * [7] M. Andrle, L. Rebollo-Neira, E. Sagianos, “Backward-optimized orthogonal matching pursuit approach”, IEEE Signal Proc. Let.,11, 705–708 (2004).
  * [8] M. Andrle and L. Rebollo-Neira, “A swapping-based refinement of orthogonal matching pursuit strategies”, Signal Processing, 86, 480–495 (2006).
  * [9] D. L. Donoho , Y. Tsaig , I. Drori , and J. Starck, “Stagewise Orthogonal Matching Pursuit”, IEEE Transactions on Information Theory, 58, 1094–1121 (2006).
  * [10] D. Needell and J.A. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples”,  Applied and Computational Harmonic Analysis, 26, 301–321 (2009).
  * [11] D. Needell and R. Vershynin, “Signal Recovery From Incomplete and Inaccurate Measurements via Regularized Orthogonal Matching Pursuit”, IEEE Journal of Selected Topics in Signal Processing, 4, 310–316 (2010).
  * [12] Y. Eldar, P. Kuppinger, and H. Biölcskei, “Block-Sparse Signals: Uncertainty Relations and Efficient Recovery”, IEEE Trans. Signal Process., 58, 3042–3054 (2010).
  * [13] L. Rebollo-Neira and J.Bowley, Sparse representation of astronomical images, Journal of The Optical Society of America A, 30, 758–768 (2013).
  * [14] L. Rebollo-Neira, R. Matiol, and S. Bibi, “Hierarchized block wise image approximation by greedy pursuit strategies,” IEEE Signal Process. Letters, 20, 1175–1178 (2013).
  * [15] L. Rebollo-Neira, “Cooperative Greedy Pursuit Strategies for Sparse Signal Representation by Partitioning”, Signal Processing, 125, 365–375 (2016).
  * [16] J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Algorithms for simultaneous sparse approximation. Part I: Greedy pursuit”, Signal Processing, 86, 572–588 (2006).
  * [17] L. Rebollo-Neira, “Constructive updating/downdating of oblique projectors: a generalisation of the Gram-Schmidt process”, J. Phys. A: Math. Theor. 40, 6381–6394 (2007).
  * [18] J. Pan and W.J. Tompkins,“A Real-Time QRS Detection Algorithm”, IEEE Trans. on Biomedical Engineering, BME-32, 230–236 (1985).
  * [19] H. Sedghamiz, Complete Pan Tompkins Implementation ECG QRS detector (https://www.mathworks.com/matlabcentral/fileexchange/45840-complete-pan-tompkins-implementation-ecg-qrs-detector), MATLAB Central File Exchange. Retrieved June 4, 2024.
  * [20] L. Rebollo-Neira and D. C̆erná, “Wavelet based dictionaries for dimensionality reduction of ECG signals”, Biomedical Signal Processing and Control, 54, 101593 (2019).
  * [21] D. C̆erná and L. Rebollo-Neira, “Construction of wavelet dictionaries for ECG modelling”, MethodsX, 8, 101314 (2021).
  * [22] A. Cohen, I. Daubechies, and J.C. Feauveau, “Biorthogonal bases of compactly supported wavelets”, Comm. Pure and Appl. Math., 45, 485–560 (1992).
  * [23] http://www.nonlinear-approx.info/examples/node017.html
  * [24] https://physionet.org/physiobank/database/mitdb/
  * [25] L. Rebollo-Neira, “Effective high compression of ECG signals at low level distortion”, Scientific Reports, 9, No 4564 (2019).
  * [26] Karl Skretting. Huffman Coding and Arithmetic Coding (https://www.mathworks.com/matlabcentral/fileexchange/2818-huffman-coding-and-arithmetic-coding), MATLAB Central File Exchange. Retrieved June 4, 2024.
  * [27] A. Adamo, G. Grossi, R. Lanzarotti, J. Lin, “ECG compression retaining the best natural basis k-coefﬁcients via sparse decomposition”, Biomedical Signal Processing and Control, 15, 11–17 (2015).
  * [28] T. Yu Liu, K. Jen Lin, and H. Chun Wu, “ECG Data Encryption Then Compression Using Singular Value Decomposition”, IEEE Journal of Biomedical and Health Informatics, 22, 707 – 713 (2018).
  * [29] L. Zheng, Z. Wang, J. Liang, S. Luo, and S. Tian, “Effective compression and classiﬁcation of ECG arrhythmia by singular value decomposition”,Biomedical Engineering Advances, 2, 100013 (2021).
  * [30] S.J. Lee, J. Kim, and M. Lee,“A Real-Time ECG Data Compression and Transmission Algorithm for an e-Health Device”, IEEE Trans. on Biomedical Engineering, 58, 2448 – 2455 (2011).
  * [31] J.L. Ma, T.T. Zhang, and M. C. Dong, “A Novel ECG Data Compression Method Using Adaptive Fourier Decomposition With Security Guarantee in e-Health Applications”, IEEE Journal of Biomedical and Health Informatics, 19, 986 – 994 (2015).
  * [32] C. Tan, L. Zhang and H. Wu,“A Novel Blaschke Unwinding Adaptive Fourier Decomposition based Signal Compression Algorithm with Application on ECG Signals”, IEEE Journal of Biomedical and Health Informatics 23, pp 672 – 682 (2018).