# Filtering through a topological lens: homology for point processes on the time-frequency plane

Juan Manuel Miramont, Kin Aun Tan, Soumendu Sundar Mukherjee, Rémi Bardenet, Subhroshekhar Ghosh

## ► To cite this version:

Juan Manuel Miramont, Kin Aun Tan, Soumendu Sundar Mukherjee, Rémi Bardenet, Subhroshekhar Ghosh. Filtering through a topological lens: homology for point processes on the time-frequency plane. 2025. hal-05030402

HAL Id: hal-05030402

https://hal.science/hal-05030402v1

Preprint submitted on 11 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Filtering through a topological lens: homology for point processes on the time-frequency plane

Juan Manuel Miramont\*†, Kin Aun Tan\*‡, Soumendu Sundar Mukherjee\*§, Rémi Bardenet\*†, Subhroshekhar Ghosh\*‡

## Abstract

We introduce a very general approach to the analysis of signals from their noisy measurements from the perspective of Topological Data Analysis (TDA). While TDA has emerged as a powerful analytical tool for data with pronounced topological structures, here we demonstrate its applicability for general problems of signal processing, without any a-priori geometric feature. Our methods are well-suited to a wide array of time-dependent signals in different scientific domains, with acoustic signals being a particularly important application. We invoke time-frequency representations of such signals, focusing on their zeros which are gaining salience as a signal processing tool in view of their stability properties. Leveraging state-of-the-art topological concepts, such as stable and minimal volumes, we develop a complete suite of TDA-based methods to explore the delicate stochastic geometry of these zeros, capturing signals based on the disruption they cause to this rigid, hyperuniform spatial structure. Unlike classical spatial data tools, TDA is able to capture the full spectrum of the stochastic geometry of the zeros, thereby leading to powerful inferential outcomes that are underpinned by a principled statistical foundation. This is reflected in the power and versatility of our applications, which include competitive performance in processing a wide variety of audio signals (esp. in low SNR regimes), effective detection and reconstruction of gravitational wave signals (a reputed signal processing challenge with non-Gaussian noise), and medical time series data from EEGs, indicating a wide horizon for the approach and methods introduced in this paper.

The interplay of randomness and structure is a defining characteristic of modern day data science and machine learning. The size and dimensionality of data we can access and store is vast and ever-growing, potentially overwhelming existing computational capacities. This has motivated the quest for simple, low-dimensional structure in the randomness of massive datasets, which is an effective ansatz to capture how the real world manifests itself in big data.

**Topological Data Analysis.** An effective tool that has emerged in recent years to identify and analyze such structure is *Topological Data Analysis* (abbrv. TDA). Leveraging ideas and techniques from algebraic topology, this suite of methods aims to understand latent geometric structures that are often present in large-scale datasets. By developing a principled multi-scale approach to understand the topological structure of data point clouds, TDA isolates geometric features (such as holes) that *persist* across a large number of scales (a technique referred to as *persistent homology*). These persistent features may be identified to be inherent or intrinsic to the dataset that encapsulates its latent geometric structure, which in turn reveals the fundamental real-world information hidden in the ambient noise and many confounding attributes.

Originating in the late 20th century, TDA has quickly emerged to demonstrate practical utility in a wide range of applications, starting from materials science [1, 2] and genomics [3] to fingerprinting proteins [4] and audio signals [5], biomedicine [6], neuroscience [7] and spatial networks [8, 9], to provide a very limited and partial list of successful use cases. However, while TDA provides an attractive combination of powerful and sophisticated mathematical tools brought to the service of real-world data, its many successful applications have mostly been focused on situations where the data is fundamentally *geometric* in its nature, or at least has an underlying notion of *shape* in a

<footnote>\*Corresponding authors: J. M. Miramont, K. A. Tan, S. S. Mukherjee</footnote>

<footnote>E-mails: jmiramontt[at]univ-lille[dot]fr; kinaun[dot]tan[at]u[dot]nus[dot]edu; ssmukherjee[at]isical[dot]ac[dot]in</footnote>

<footnote>†Université de Lille, CNRS, Centrale Lille, UMR9189-CRISTAL, Av. Paul Langevin, 59650 Villeneuve-d'Ascq, France.</footnote>

<footnote>‡Department of Mathematics, National University of Singapore, 10 Lower Kent Ridge Road, Singapore 119076.</footnote>

<footnote>§Statistics and Mathematics Unit, Indian Statistical Institute, 203 Barrackpore Trunk Road, Kolkata, India 700108.</footnote>
reasonably straightforward sense [10, 11]. However, most real-world data generating processes do not
embody a clear geometric (or shape) structure in any direct manner. This naturally leads to the
tantalizing question of bringing the formidable arsenal of TDA techniques to bear on wider classes of
data analytical problems without an a priori geometric structure.

In this work, we focus on a piece of this puzzle, proposing a principled approach towards the large-scale application of TDA methods to a very general class of signal processing problems. In particular, our approach is readily applicable to a broad bandwidth of time-indexed signals, including as a special case all kinds of noisy signals arising in acoustics. We show the range of applicability extends much farther, encompassing gravitational wave signals (that are known for their notorious difficulty due to extremely low SNR and non-Gaussianity of noise), and general classes of time-series data (such as medical data from EEG measurements). Our methods are underpinned with a principled statistical foundation, including in particular specialized ingredients from multiple hypothesis testing to inform our signal reconstruction procedure. Further connections to related work [5, 12], are addressed later in the Discussion (Sec. 8).

**Time-frequency signal processing.** A more in-depth discussion calls forth a brief digression into time-frequency signal processing, a suite of techniques that distinguishes itself by studying square-integrable functions of time (known as *finite-energy signals*), tracking frequency as it varies over time in the signal.

One can think of someone whistling a melody they are reading on a musical score: time-frequency
analysis aims at reconstructing the musical score from an audio recording of the whistler. This recon-
structed musical score, known as a time-frequency representation of the signal, can then serve in many
downstream tasks, such as detecting whether someone was indeed whistling, and if yes, identifying the
song from which the melody is taken. Time-frequency analysis is not limited to audio processing: a
modern success story of the field is the detection of gravitational waves, which have a well-understood
time-frequency signature in the interferometric signals collected by the LIGO and VIRGO detectors
[13]; other applications abound in a multitude of fields such as medical image processing [14], time
series [15], analysis of financial data [16], biomedical signals [17, 15] and biomechanics of postures [15].

One of the dominant time-frequency representations is the short-time Fourier transform (STFT). It is a complex-valued function of two variables, time and frequency, which collects Fourier transforms of the signal multiplied by a sliding window. Following the seminal work [18], a string of recent works [19, 20, 21, 22, 23, 24, 25, 26, 27] has focused on extracting information solely from the zeros of this transform, informally time-frequency pairs that are points of *pure silence*; see also [28] for an introductory survey. This approach is motivated in part by results that indicate that these zeros almost characterize the whole transform [29, 18], and in part by experiments hinting that detection using zeros comes with better statistical power than maxima-based methods at low signal-to-noise ratio (SNR) [30, 31]. In short, the spectrogram zeros form a highly parsimonious, albeit non-linear, summary of the landscape of the STFT that is nonetheless informative enough to perform meaningful signal processing.

**Zeros, hyperuniformity and statistical mechanics.** Spectrogram zeros bring into play con-
nections to the theory of point processes and strongly correlated particle systems from statistical
mechanics, which we expound on below. To wit, when the signal is modeled by a random process,
the zero set of its STFT becomes a point process, i.e. a random configuration of points [20, 21, 24].
When the signal is pure white Gaussian noise, this point process turns out to correspond to well-
known models of strongly correlated particle systems in statistical mechanics, namely the zeros of the
so-called planar (real) Gaussian analytic function (and its relatives). In particular, such systems are
known to exhibit a very high degree of homogeneity in their spatial distribution, a property known as
*superhomogeneity* or *hyperuniformity* [32, 33, 25, 34, 35, 36]. An illustrative sample can be found in
Fig. 1a.

When a deterministic signal is added, the modulus of the STFT becomes large in a time-frequency
region informally called the *support* of the signal, effectively creating signal-shaped holes in the pat-
tern of zeros; see e.g. Fig. 1b. Holes in a hyperuniform, grid-like pattern [33, 35] like the zeros of the
planar Gaussian analytic function are intuitively easier to detect than similar holes in a completely
random point pattern such as a Poisson point process. This seminal observation of [18] has motivated
the use of techniques from spatial statistics to detect and extract signals by detecting and locating
unusual holes in the pattern of zeros. Yet, such methods have so far relied on relatively rudimentary
test statistics (such as pairwise distances [20, 30]) that often fall far short of capturing the intricate
stochastic geometry of the zeros, often resulting in a collection of arbitrary algorithmic choices, leading to an ill-posed question of defining a significant hole from a collection of generic spatial statistics. This, therefore, calls for a principled, automated approach to detecting and reconstructing holes of arbitrary shape in the configuration of zeros of a spectrogram, which is a problem that we will focus on.

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Spectrogram of a realization of complex white Gaussian noise (white dots indicate the spectrogram zeros). (b) Same as (a) but with an additive mixture of signal and noise.</figcaption></figure>

**TDA meets signal processing.** The inherent capability of TDA to access the finer topological properties of a point cloud, and the intrinsic requirement for this capability in signal analysis based on spectrogram zeros, make the intersection of these two domains a rich and natural overlap for symbiotic development – a direction that we pursue in this work. To sum up, the crucial problem in zero-based signal processing is to algorithmically detect and reconstruct statistically significant holes in an efficient manner; on the other hand, the algebraic topology of the two-dimensional spectrogram zero point cloud is fundamentally determined by its holes (in technical terms, its homology), which smoothly dovetails the algorithmic developments in persistent homology. Reciprocally, this overlap also raises fundamental and interesting statistical problems in TDA (such as the determination of the number and identities of statistically significant holes from persistence diagrams), some of which we take up in this paper.

Furthermore, as discussed at the beginning, this brings the powerful tools and techniques of TDA to bear on a very general class of time-indexed signal processing problems where a pronounced geometric shape structure is not present in the *a priori* description of the problem. We hope that this work will inaugurate a systematic development of the mutual interconnections between TDA and signal processing, motivating the growth of new tools and techniques in both fields.

# 1 Spectrograms and their zeros

The short-time Fourier transform (STFT) is a fundamental concept in time-frequency analysis [37, 38]. For a signal $f \in L^2(\mathbb{R})$ and a window $g \in L^2(\mathbb{R})$, the STFT of $f$ with window $g$ is defined as

$$V_g(f)(u, v) = \int_{-\infty}^{+\infty} f(t) \overline{g(t-u)} e^{-2i\pi tv} dt, \quad (1)$$

where the variables $u, v \in \mathbb{R}$ are usually referred as *time* and *frequency*. Taking $\|g\| = 1$, $V_g$ is an isometry, and the *spectrogram* $(u, v) \mapsto |V_g(f)(u, v)|^2$ is considered a measure of how the *energy* $\|f\|^2 = \|V_g(f)\|^2$ is spread across the time-frequency plane $\mathbb{R}^2$ [39].

In the fundamentally important case when the window is taken to be the Gaussian $g: t \mapsto 2^{1/4} \exp(-\pi t^2)$, $V_g$ maps $S'(\mathbb{R})$ to analytic functions of $z = u + iv$ [40, 20]. In particular, since analytic functions have isolated zeros, one can talk about the random configuration of isolated points (in probabilistic terms, the *point process*) formed by the zeros of $z = u + iv \mapsto V_g(f + \xi)(u, v)$, where $\xi$ is a Gaussian white noise. Finally, we define the signal-to-noise ratio (SNR) of a signal-plus-noise mixture $h = f + \xi$ as

$$\text{SNR}(h) = 10 \log_{10} \left( \frac{\|f\|^2}{\mathbb{E}[\|\xi\|^2]} \right) \quad (\text{dB}). \quad (2)$$
<figure><img src="image_3.png" /><figcaption>**Figure 2:** (a) The spectrogram of the gravitational wave event GW150914 (top), recorded by the Laser Interferometer Gravitational-wave Observatory (LIGO) in Hanford, WA, USA. and the corresponding signal (bottom). (b) Filtration of $\alpha$-complexes showing the birth of the most persistent cycle in red line. (c) Persistence diagram. The most persistent component is marked in red and the corresponding cycle is superimposed to the spectrogram for reference. $d_{(1)}$ is the longest distance to the identity diagonal.</figcaption></figure>

Figure 1a and Fig. 1b show the log spectrogram of resp. a realization of Gaussian white noise and a deterministic signal plus Gaussian white noise, with white dots indicating the zeros of the spectrogram. As first noted by [18], the zeros are absent in the signal’s time-frequency support, so that identifying the regions with low density of zeros is tantamount to identifying areas of the time-frequency plane where the signal is strong.

# 2 Tools from topological data analysis

Persistent homology is arguably the main workhorse of topological data analysis, providing a principled understanding of the multiscale geometry of a data point cloud, and significantly, enabling algebraic manipulations that are amenable to computer implementation. In the limited scope of the main text, we provide a succinct and conceptual overview focusing mostly on the 2D setting relevant for us, directing the interested reader to Section A of the Appendix for a more detailed and in-depth discussion. For even more extensive coverage, we refer to any of the excellent texts [41, 42, 43, 44], to provide a partial list of references.

## 2.1 Persistence and Persistence Diagrams

A very canonical approach to TDA proceeds via developing a *geometric graph* (more precisely, a *simplicial complex*) based on the given data point cloud, by connecting data points that are closer to each other than a threshold $r \ge 0$. For contextualization, we refer the reader to the illustrations in Fig. 2b. The threshold $r$ will be varied continuously, giving rise to a multiscale sequence of nested geometric graphs (*a filtration*), which at any particular threshold value gives a skeletal representation of the dataset at that scale.

Starting with the isolated data points when the threshold $r$ is equal to 0, the graph grows and connectivity increases with increasing $r$, eventually culminating in the fully connected graph (c.f. Fig. 2b). At some point, cycles start forming due to edges being connected; this stage may be thought of as the *birth* of a cycle (to be thought of as a *hole*). At this point, the $r/2$ disks centered at the endpoints of each edge of the cycle overlap pairwise, but all the disks pertaining to the cycle do not overlap jointly. At a suitably higher threshold $r$, all these disks do indeed overlap, in which case we declare the *hole* to be *filled* and the cycle to have died. It is reasonable to surmise that there are two kinds of cycles aka *holes* – those that get *filled* i.e. die soon after they are created, and those that *persist* for a long time between their birth and death. Persistent homology theory identifies the latter variety of cycles
or holes to be fundamental characteristics of the topological structure of the dataset, relegating the former type to be an artifact of random noise. The so-called *Persistence Diagram* (*abbrv.* PD) plots the times of birth and death of the different cycles against each other in a planar diagram; as such, the points in the PD that are far from the diagonal are understood to carry fundamental structural information on the dataset. Figure 2c displays the PD for the gravitational wave signal from Fig. 2a, where is possible to appreciate that the most persistent hole, i.e. the point further from the identity in the PD, corresponds with the hole created by the signal in the pattern of zeros.

## 2.2 Persistent Homology and Learning Holes

The homological ingredient in Persistent Homology (*abbrv.* PH) brings into play a crucial linear algebraic toolbox by considering the abstract vector spaces spanned by the cycles, edges and vertices in this (sequence of) geometric graphs, and by looking at the *boundary maps* that connect each of these objects to their extremities. Under the standard *homotopy equivalence* of cycles in algebraic topology, one can obtain natural (finitely generated) Abelian group structures on these vector spaces with the boundary maps inducing homomorphisms connecting them, with the groups corresponding to different thresholds being further connected via canonical maps arising out of natural embeddings of the simplicial complexes. An in-depth description of this can be found in Sec. A.1 of the Appendix.

The upshot of this is that efficient linear algebraic methods can be invoked to study these homology groups, wherein a non-trivial homology group (of first order) corresponds to the presence of holes that equal in number to the rank of this group (counted modulo homotopy equivalence). In particular, there are effective tools to compute the *generators* of these homology groups, which themselves correspond to (homotopy equivalence classes of) holes. This motivates our broad approach to the identification of non-trivial holes in a point cloud: identify points on the PD that are far from the diagonal (in a statistically significant sense). Such an *outlier* in the PD would pertain to a cycle in the relevant homology group, which would lead us to a hole in the underlying point cloud.

## 2.3 Optimal Cycle Representations

A key issue in the above approach to learning holes is the inherent ambiguity in the homological cycles due to them being merely homotopy equivalence classes of cycles in the original point cloud (as opposed to being the cycles themselves). This leads to the question of identifying suitable representative cycles in these homological equivalence classes, culminating in the notion of *Minimum* and *Stable* Volumes [45, 46], that we will use in our algorithmic considerations.

### 2.3.1 Minimum Volumes

It is intuitively appealing to search for a representative cycle for the equivalence class that encloses the smallest possible area of $\mathbb{R}^2$ in the spectrogram, that portion would likely correspond to the signal support. One way to formalize this minimizing class representer is the notion of *minimum-volume cycle*, introduced in [47] and generalized in [45]. While the original definition takes the form of the solution to a constrained optimization problem, the solution of which has to be numerically approximated, our particular setup allows us to easily obtain minimum-volume cycles through the so-called *persistence tree* (PT) of the point cloud, a notion introduced by [47]. Notice that in our application, we are more interested in the area enclosed by the minimum-volume cycles, henceforth termed minimum volume (*abbrv.* MV). For an algorithm to generate the PT, we refer the reader to the Algorithm 1 in the Appendix and [45]. A detailed discussion on MVs is taken up in Sec. A.2 of the Appendix.

### 2.3.2 Stable Volumes

It is a well-known issue in computational topology that a small change in the input signal, such as random noise, can produce large variations of the cycles –and their volumes–, even for the cycles associated to long-living components [48, 49, 46]. Approaches to *stabilize* the components of interest have been proposed [49, 46]. Intuitively, one way to robustify the construction of minimum volumes is to take the intersection of minimum volumes corresponding to slightly altered inputs, where the alteration is small enough that we can identify a volume from one signal to the other. Formally, we focus here on the so-called *stable volumes* (*abbrv.* SV) recently proposed by Obayashi [46], which can
be defined directly from the persistence tree described in Sec. A.2 of the Appendix. Moreover, a detailed discussion on stable volumes and its computation is undertaken in Sec. B of the Appendix as well.

# 3 Ingredients from Statistical Inference

We give in this short section a few elements from statistical inference that we later use in signal detection and reconstruction. We heavily lean on multiple hypothesis testing to make informed decisions in our signal processing algorithms. In simple hypothesis testing, one has one null hypothesis $H_0$ and an associated $p$-value for a given observation. We say we reject $H_0$ when the $p$-value is lower than a prescribed significance level $\alpha$. Intuitively, $\alpha$ is the probability of incorrectly rejecting $H_0$ (i.e. type I error) that, as a researcher, one is willing to accept.

Multiple hypothesis testing is a more general scenario where several tests are performed for a number of null hypotheses. Say we perform $K$ tests, this is akin to having $K$ null hypotheses $H_{1,0},..., H_{K,0}$ and determining if any of the $K$ $p$-values is below some per-test significance level $\alpha_k$.

A relevant concept in this case is the family-wise error rate (FWER), which is the worst-case probability of rejecting a true hypothesis among $H_{1,0},..., H_{K,0}$, i.e. the probability of making any type I error while conducting all $K$ tests [50]. Controlling the FWER means ensuring that the chance of falsely rejecting any null hypothesis is less than $\alpha$ [51, 52]. If $\alpha_k \equiv \alpha$, then from basic probability we know that the FWER is higher than $\alpha$ itself. Thus, some sort of correction is needed. A fairly used approach is the Bonferroni correction, where $\alpha_k \equiv \frac{\alpha}{K}$, which we shall use in the following sections.

For a more elaborate discussion of the concepts involved in multiple testing, we refer the reader to Sec. C of the Appendix.

# 4 Signal detection via Persistence Diagrams

We formulate several detection tests based on different choices of test statistics. We first describe the simple and multiple testing procedures for a generic test statistic, to be instantiated later with test statistics based on TDA.

## 4.1 Simple and multiple hypothesis testing

Signal detection can be stated as a binary test [37, 20, 24]:

$$ \begin{cases} H_0 : h = \xi, & \text{i.e. pure noise} \\ H_1 : h = f + \xi, & \text{i.e. signal plus noise} \end{cases} $$

where the null model proposed by $H_0$ simply states that the observed data $h$ is purely noise, whereas the alternative states that the data is actually a mixture of signal and noise.

### 4.1.1 Simple hypothesis test

Let $X$ be a generic test statistic that takes higher values for signals than for pure noise. We take $B$ independent Monte Carlo simulations of white noise and compute the test statistics $X_1,..., X_B$. A $p$-value for testing $H_0$ is then computed as $p = \frac{1}{B} \sum_{j=1}^B \mathbb{I}(X_{(j)} > X_{\text{obs}})$, where $X_{(j)}$ are the top order statistics of $X$ and $X_{\text{obs}}$ is the test statistic obtained from the sample to be tested.

### 4.1.2 Multiple hypothesis testing - simultaneous tests

We devise an alternative test for detection by considering the $K$ most persistent components in the PD rather than just the most significant one as before. Start with a suitable choice of $K$. Again, from $B$-many simulated persistence diagrams, with statistics $X_{(k),j}, 1 \le j \le B$, we compute a series of $p$-values $p_1,p_2,...,p_K$ as

$$ p_k := \frac{1}{B} \sum_{j=1}^{B} \mathbb{I}(X_{(k),j} > X_{(k),\text{obs}}). \qquad (3) $$
We conduct these tests at varying levels of significance $\alpha_k$. A signal is detected if and only if $p_k < \alpha_k$ for some $1 \le k \le K$. We choose $\alpha_k = \frac{\alpha}{K}$ according to the Bonferroni correction. More details of this testing strategy can be seen in Sec. C.2.2 of the Appendix.

The following subsections define the test statistics we derive from TDA: distance to the diagonal in the PD and energy from the most persistent components.

## 4.2 Distance-based test statistic

Formally, for a cycle c, let $d(c) = (y(c)-x(c))/\sqrt{2}$ denote the distance of the birth-death pair $(x(c), y(c))$ to the diagonal line in the PD. We compute the order statistics $d_{(1)} \ge d_{(2)} \ge \dots$ of these distances, where $d_{(1)}$ is the longest distance to the diagonal in the PD (as shown in Fig. 2c), $d_{(2)}$ is the second longest and so on. Possible test statistics in this case are $d_{(1)}$ (for simple hypothesis testing) or the K most significant distances (for multiple hypothesis testing).

## 4.3 Energy-based test statistic

Aside from $d_{(k)}$, we have derived an alternative test statistic based on the energy $\mathcal{E}_{(k)}$ contained in the K most persistent components according to the distances to the diagonal line in the PD.

In practice, spectrograms are approximated by Riemann sums at a finite grid of times and frequencies, yielding a matrix $S = (s_{mn}) \in \mathbb{R}_+^{M \times N}$, where $m \in \{0, 1, \dots, M-1\}$ is the frequency bin and $n \in \{0, 1, \dots, N-1\}$ is the time index. Thus $S = s_{mn} \approx |V_g(f)(\frac{n}{N}T_s, \frac{m}{2MT_s})|^2$ for a sufficiently small sampling period $T_s$.

We compute $\mathcal{E}_{(k)}^{\text{MV}}$ (resp. $\mathcal{E}_{(k)}^{\text{SV}}$), the energy contained in the minimum (resp. stable) volume of the k-th most significant component, as

$$ \mathcal{E}_{(k)}^{\text{MV}} = \frac{\sum_{m,n} (1_{mn}^{\text{MV}(k)} s_{mn})}{\sum_{m,n} s_{mn}}, \quad \mathcal{E}_{(k)}^{\text{SV}} = \frac{\sum_{m,n} (1_{mn}^{\text{SV}(k)} s_{mn})}{\sum_{m,n} s_{mn}}, $$

where $(1_{mn}^{\text{MV}(k)})$ and $(1_{mn}^{\text{SV}(k)})$ are the binary masks indicating the time-frequency region in the spectrogram occupied by the minimum volume and the stable volume of the k-th most significant component, respectively.

## 4.4 Energy vs distance based tests

Depending on the nature of the dataset, hypothesis testing based on energies might outperform hypothesis testing based on distances to the diagonal. For instance, signals containing a noticeable amount of energy but in a short duration and narrow frequency range occupy a small region in the time-frequency plane that is less detectable using distance statistics. The comparison between the two choices of test statistics is further elaborated and demonstrated in the experiments sections (Secs. 6 and 7).

## 4.5 Accumulated persistence function

In [8], Biscio et al. proposed a Monte Carlo envelope test based on a function of the PD called the *accumulated persistence function* (APF). APF condenses the information from the persistence diagram obtained from the observation, or a simulation, into a single real-valued function. We use a Monte Carlo *global* test [53, 54] based on contrasting the statistics of the APF of simulations and data in order to detect a signal. Details of this test can be seen in the Sec. D of the Appendix.

# 5 Signal reconstruction via Persistence Diagrams

Finding persistent holes in the zeros of the spectrogram allows the identification of regions where the signal is stronger than noise, the so-called *signal domain D* in the TF plane [18, 20]. In this section we explain how to reconstruct a signal from its noisy samples based on an estimation of D.
<figure><img src="image_4.png" /><figcaption>Figure 3: Example of component identification on a bat echolocation signal. The number superimposed on each component indicates the relevance of the component in the persistence diagram shown on the right.</figcaption></figure>

## 5.1 Signal reconstruction

In order to properly estimate $f$, we first need to determine the number of components, i.e. persistent holes in the pattern of spectrogram zeros, which we use to approximate $D$. As an example, Figure 3 shows the minimum-volume cycles corresponding the three most persistent components of a bat echolocation signal, and their corresponding birth-death pairs in the PD. Ideally, $D$ corresponds to the union of the associated minimum (or stable) volumes. Once $D$ is estimated, the signal itself can be estimated as [18]:

$$ \hat{f}(t) = \frac{1}{g(0)} \int_{-\infty}^{+\infty} V_g(h)(t, \nu) \mathbb{I}_D(t, \nu) e^{2i\pi t\nu} d\nu. \quad (4) $$

This leads to the question of how to estimate the number of holes, which we address in the following section.

## 5.2 Estimation of the number of holes

Counting the number of holes is tantamount to the problem of detecting several signals. Hence, we use once again multiple hypothesis testing but now with the goal of estimating an appropriate number of relevant components, rather than just claim the existence of a signal as we did in Sec. 4.

### 5.2.1 Multiple hypothesis testing - sequential tests

As before, from the $B$-many noise realizations we obtain independent test statistics $X_{(k),1}, \dots, X_{(k),B}$ and compute $K$ $p$-values $p_k$ using (3). Now, in contrast to what we do for simultaneous tests, we conduct these tests sequentially, i.e. one-by-one, at varying levels of significance $\alpha_k$. If $p_1$ is significant, there is evidence for the presence of signal (as in the simple hypothesis case). For $k \ge 1$, if $p_k < \alpha_k$ we go on to conduct the test for $k + 1$. We stop at $k_*$, the smallest $k$ such that $p_k \ge \alpha_k$, declaring that our estimate of $n_{\text{holes}}$, the total number of holes or persistent components to extract, is $k_* - 1$. Thus $\hat{n}_{\text{holes}} = \min\{k \ge 0 : p_{k+1} \ge \alpha_{k+1}\}$. Notice that the $p$-values are computed for up to $K$ components. The tests can be carried out by taking the test statistics to be the distances $d_{(k)}$ to the diagonal line in the PD, or the energies $\mathcal{E}_{(k)}^{MV}$ (or $\mathcal{E}_{(k)}^{SV}$) with $p$-values computed in a similar way as before.

### 5.2.2 Choice of levels of significance

We consider three different choices of $\alpha_k$ for sequential testing: (i) Polynomial decay: $\alpha_k = \frac{\alpha}{k^m}$, $m \ge 0$; (ii) Geometric decay: $\alpha_k = \alpha\beta^k$, $\beta \in (0, 1)$ and (iii) Bonferroni correction. Both polynomial and geometric decay can be shown to control the FWER under the assumption that the $p$-values are independent. Further details of such corrections can be seen in Sec. C.2.3 of the Appendix.
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Statistical power of the detection tests described in Sec. 4 as a function of the SNR. Each color and marker indicate a different test, the name of which is indicated on top of the figure. Error bars indicate the 95% Clopper-Pearson confidence intervals, Bonferroni-corrected across the 32 comparisons (4 SNR levels × 8 approaches). The spectrogram of the probed signal with SNR = 10 dB is superimposed.</figcaption></figure>

# 6 Applications: Detection of synthetic and real-world signals

This section describes the application of the proposed approaches for signal detection based on TDA and the zeros of the spectrogram. First, the results of the tests for signal detection considering a synthetic signal are reported. Then we illustrate the performance of the detection tests on gravitational wave signals. Alpha complexes and persistence diagrams were computed using GUDHI's Python front-end [55]. The code to reproduce the experiments in this paper is provided<fnref n="1" />.

## 6.1 Synthetic chirp detection

### 6.1.1 Test signal specification

Following similar experiments from [20, 30, 28], we first demonstrate the performance of the detection tests described in Sec. 4 on a discrete signal of length $N = 1024$ consisting of a linear chirp of duration $N/2$ centered in the middle of the observation window. Figure 4 shows the statistical power of different tests for several SNRs. The spectrogram of the probed signal is also shown in the figure.

### 6.1.2 Detection tests

We compared a total of eight tests: six variations of TDA-based tests (Sec. 4), a test based on the APF (Sec. 4.5), and a Monte Carlo test based on ranked envelopes (taken from [31]). We consider the latter two as a baseline for signal detection tests based on spectrogram zeros.

Two TDA-based tests rely on the distance to the diagonal in the PD (Sec. 4.2). One uses the distance to the diagonal of the furthest birth-death pair in the PD, i.e. the first or most relevant component (denoted “First component-dist” in Fig. 4). The other searches for a statistically relevant component among the $K = 30$ components with the largest distance to the diagonal, using a Bonferroni correction to control the type I error (“Bonferroni-dist” in Fig. 4).

The remaining four tests based on TDA use the energy associated with the components from the PD as test statistic (Sec. 4.2).

The energy was computed from the MVs or the SVs (defined in Sec. 2.3) corresponding to the first $K$ persistent components.

<footnote n="1">¹https://github.com/rbardenet/persistence</footnote>
### 6.1.3 Detection results

Figure 4 shows that the proposed tests based on the energy of the first $K = 30$ significant components have the highest detection power from all the evaluated tests ("Bonferroni-energy-SV" and "Bonferroni-energy-MV"). From these, the tests based on stable volumes ("Bonferroni-energy-SV") are more powerful than the one based on minimum volumes, specially for low SNRs. These tests have a better performance than the baseline ("MC Rank test") for all the considered signal lengths. The two other tests based on the energy ("First component-energy-MV" and "First component-energy-SV") have a very similar performance.

In comparison, the tests based solely on the distance to the diagonal ("First component-dist" and "Bonferroni-dist") and the test based on the APF show lower power than the baseline for the detection of a single component, even for high SNRs.

## 6.2 Gravitational waves

### 6.2.1 The Gravitational Waves problem

Gravitational waves (GWs) are spacetime perturbations produced by astrophysical events such as merging black holes. Their detection enables tests of general relativity and studies of extreme cosmic phenomena. However, detecting GWs is a formidable challenge because of their exceptionally low SNR.

For illustration purposes, we have selected GW data from the *Kaggle Gravitational Waves Detection competition* [56]. The dataset consists of realistic simulations of measurements from interferometers. Each data sample contains a signal spanning 2 seconds and sampled at 2048 samples-per-second.

### 6.2.2 Stratifying GW signals via adaptive thresholding

Since SNR values are not available, we compute a hard threshold [57, 58, 59] on the spectrogram of each signal in order to segment the dataset.

We define the normalized spectrogram maximum (NSM), $\hat{\zeta} = \hat{\gamma}/\hat{\sigma}$, where $\hat{\gamma} = \max_{u,v} |V_g(f)(u,v)|$ is the maximum modulus of the STFT across the observation window on the time-frequency plane, and $\hat{\sigma}$ is the standard deviation of the noise, estimated from the signal using a median absolute deviation estimator [57, 58].

The NSM serves as a practical proxy for SNR; signals with higher NSM values are expected to have higher SNR values. We then consider the class of GW signals, which is the union of the subsets $Z_i$ for $i \in \{2, 3, 4, 5, 6, 7\}$, where each $Z_i$ contains signals with $\hat{\zeta} \ge i$. Note that $Z_i \subset Z_j$ if $i \ge j$.

We collected 50,000 signal samples for $Z_2$, from which $|Z_3| = 49,995$, $|Z_4| = 14,542$, $|Z_5| = 1,528$, $|Z_6| = 670$, and $|Z_7| = 255$. Even with a threshold of $\hat{\zeta} = 5$, we retained a large enough number of signals for robust analysis. Examples of spectrograms of the signals in the different subsets can be found in Sec. E of the Appendix. Most signals in $Z_7$ are visually identifiable in their spectrograms, while $Z_5$ and even $Z_6$ already contain numerous more subtle signals.

We observe that our NSM thresholding is only used to create subsets of the data for comparison purposes; it plays no role in the detection procedure whatsoever. We also note that our approach is a general signal analysis technique that is being applied to the challenging problem of GW data as a particular application, as opposed to most state-of-the-art tools that are highly specialized and tailored to the GW problem. We envisage that combining our methods with specialized GW data analysis tools have the potential for producing powerful signal processing outcomes for GW data.

### 6.2.3 Gravitational waves detection

Figure 5 shows the Receiver Operating Characteristic (ROC) curves for the tests based on $\mathcal{E}_{(\ell)}^{\text{SV}}$, $\ell = 1, \dots, K$, the energies contained in stable volumes associated with the $K$ most significant components (see Sec. 4.3).

For this application, we chose $K = 5$ and applied simultaneous hypothesis tests with Bonferroni correction (see Sec. 4.1.2).

ROC curves were computed using class-balanced datasets for each $Z_i$. As one would expect, signals from datasets corresponding to higher thresholds are likelier to have higher SNRs, making their features easier to detect.
<figure><img src="image_6.png" /><figcaption>**Figure 5:** ROC curves for different thresholds $\zeta \in \{2, 3, 4, 5, 6, 7\}$ using energy as test statistics. The values shown in the parentheses are the area under the ROC curves for each threshold.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Example of gravitational wave detection and estimation. (a) Detail of the spectrogram of the signal, the boundary of the stable volume corresponding to the detected component in red (the complete spectrogram is superimposed for reference). (b) Persistence diagram. The red square indicates the detected component. (c) Time domain comparison of the original and the detected, filtered signal. Notice the difference in the scale of the vertical axis.</figcaption></figure>

We achieved an area under the ROC curve of 0.97 for the highly curated dataset $Z_7$ (see Fig. 5). For subset $Z_5$, which includes a large number of more challenging examples, the area under the ROC curve is still high, to the tune of 0.84. Including the examples with the lowest SNRs reduces the performance, as our generic approach is not tailored for GW detection—unlike other competing methods that rely on matched filtering. Future work will focus on enhancing the detection power of our tests (see the discussion in Sec. 8).

## 6.3 Detection summary

In Sec. 6.1 we compared detection tests based on two different choices of test statistics. The first uses the distance to the diagonal in the PD as a test statistic, whereas the second uses the energy of the most persistent components as a test statistic. The latter turned out to be more effective for signal detection (see Fig. 4). When comparing the power of the distance-based tests ("dist" in Fig. 4) with
<figure><img src="image_8.png" /><figcaption>Figure 7: (a)-(c) Synthetic signal reconstruction performance measured by SNR gain, i.e. SNR-QRF (dB) for, using minimum and stable volumes. Bars indicate the average over 200 realizations. Error bars indicate the 95% confidence interval, Bonferroni-corrected across the 15 comparisons (3 SNR levels × 5 approaches). The proposed approach is indicated as TDA. DT stands for the Delaunay triangulation approach.</figcaption></figure>

those based on energy (“energy” in Fig. 4), it seems that large values of the spectrogram play a key role in improving detection tests based solely on the zeros of the spectrogram.

# 7 Applications: Signal Reconstruction and Audio Improvement

In this section we demonstrate how we can reconstruct signals by identifying the volumes in the TF plane corresponding to the most persistent components. We first show an example of the reconstruction of a GW from the dataset used in Sec. 6.2. Then we quantitatively characterize the performance of our approach with synthetic and audio signals.

## 7.1 Gravitational Wave Reconstruction

Figure 6a shows the spectrogram of a gravitational wave taken from the dataset. Only one component was detected, as expected, and Figure 6a shows that the boundary of the corresponding stable volume (in red), which correctly encloses the signal domain of the gravitational wave with its characteristic shape of an exponential chirp [60, 13]. Figure 6b shows the associated PD, where the birth-death pair corresponding to the detected component is highlighted in red. The original and estimated signals are shown in Fig. 6c.

In contrast to [18], our approach based on persistent homology does not depend on a parameter that fixes the minimum scale of the holes (such as a minimum length of the edges of the triangles). Instead, one needs to determine the number of components to be estimated from the signal. This is automatically and adaptively done by the detection tests proposed in Sec. 5.2, so that the method is fully based on the data.

## 7.2 Synthetic Signal Reconstruction

We show quantitative results for signal reconstruction using three synthetic signals with different time-frequency structures. Their spectrograms can be seen in Fig. 7.

### 7.2.1 Performance metric

The reconstruction performance of synthetic signals was measured using the quality reconstruction factor

$$QRF = 10 \log_{10} \left( \frac{\|f\|_2^2}{\|f - \hat{f}\|_2^2} \right) \text{ (dB)},$$
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Output audio quality measured by PESQ [61] after applying variations of the proposed approach (indicated as **TDA**) and the Delaunay triangulation (DT) based method for a short cello recording. Error bars indicate the 95% confidence intervals, Bonferroni-corrected across the 32 comparisons (4 SNR levels × 8 approaches).</figcaption></figure>

where *f* is the original noiseless signal and $\hat{f}$ is a denoised approximation [19]. QRF can be understood as an estimation of the SNR of the output signal after reconstruction, which should be higher than the SNR of the noisy input signal to indicate an improvement. Therefore the difference QRF – SNR (dB) is equivalent to the gain in terms of SNR after processing. Notice that one can only compute the QRF if the noiseless reference is available, hence QRF is only reliably computable from synthetic signals.

### 7.2.2 Results

Figure 7 compares the gain in SNR i.e. QRF – SNR (dB), using MVs and SVs for several input SNRs and three different signals. Notice that we also vary the significance of the detection tests used to count the number of components ($\alpha$ = 0.05 and $\alpha$ = 0.15). Using SVs and a higher $\alpha$ results in a better reconstruction in the case of low SNRs (for instance, 0 and 10 dB). In most cases the results obtained with TDA are better or comparable to those obtained with DT [18].

### 7.2.3 Discussion of results

The situation where one *signal* component corresponds to more than one *homology* component can occur in practice, particularly for low SNRs. Detection tests that automatically determine the number of relevant homology components to extract alleviate this issue. However, such tests are, in practice, limited by the SNR. If detection tests are not effective, the user can always determine the number of homology components to extract by hand. Even depending on *a priori* detection of components based on detection tests, signal estimation performance with the approach based on TDA is comparable or better than previous approaches based on spectrogram zeros (compare with the DT method in Fig. 7). Increasing $\alpha$ results in a larger number of homology components detected, and it could be an acceptable strategy when a signal is known to be present beforehand, with reducing the noise the primary goal rather than detecting the signal (as done in [18]).

## 7.3 Audio improvement

We now explore the reconstruction performance of an audio signal corresponding to a cello excerpt, sampled at 8 kHz and contaminated with white Gaussian noise. Figure 8 shows the spectrogram of the signal with four harmonics. In contrast to the synthetic case, there is no ground-truth, noiseless version of the real-world signal available.

### 7.3.1 Performance metric

We compare the audio quality of the estimated signal by means of the *perceptual evaluation of speech quality* (PESQ) metric [62, 61]. Although originally conceived to assess speech quality over telecommu-
nication networks, PESQ has been used as a general audio quality metric, even for non-speech signals [63]. A large value of PESQ is then interpreted as indicative of better audio quality.

### 7.3.2 Results

Figure 8 shows the PESQ values for several SNRs (from 0 dB to 30 dB). As a reference, the value of PESQ for the signal not contaminated with noise, is also shown. We explored variations of our TDA-based approach method, modifying the procedure by which the number of components to extract is determined (see step 4 in Algorithm 2 described in the Appendix).

The three $\alpha$ corrections from Sec. 5.2, i.e. Bonferroni, geometric decay and polynomial decay, were pairwise combined with the two test statistics described in Sec. 4, i.e. $d_{(k)}$ and $\mathcal{E}_{(k)}^{SV}$. The significance level was fixed to $\alpha = 0.15$ in this case.

For SNR > 10 dB, Figure 8 shows that when using “TDA-distance” the results tend to be better than when the energy is used as a test statistic (“TDA-energy”). The choice of $\alpha$ correction seems to play a less relevant role. In the low SNR regime, the use of $\mathcal{E}_{(k)}$ in combination with stable volumes increases the performance, similar to the results obtained for synthetic signals.

### 7.3.3 Discussion of results

The results obtained for the proposed approach are comparable to those obtained with the DT approach (“DT” in Figure 8) in the case of low SNR, while for higher SNRs the results are better than with DT. Figure 10a shows the minimum-volume cycles identified for the cello signal. Real-world signals might present an intricate structure in the spectrogram as well non-white noise (contrary to what it was assumed in the Monte Carlo simulations used to determine the number of components). Such factors can produce instabilities in the MVs. However, stable volumes can alleviate this issue, as depicted in Fig. 10b.

An important difference with preexisting approaches is that using TDA one can simultaneously identify individual components and their domain (see Fig. 3 or Fig. 10 in the main document) without a post-processing heuristic that groups triangles into components (as done in [19]).

<figure><img src="image_10.png" /><figcaption>Figure 9: (a) Spectrograms of EEG with no seizure present. (b) EEG with a seizure. The minimum-volume cycles corresponding to the three more persistent components are shown in each case. The number superimposed on each component indicates the relevance of the component in the persistence diagram shown on the right of each subfigure.</figcaption></figure>

# 8 Further applications, discussion and future work

## 8.1 Further applications

A key feature of our approach is the general scope of its applicability. While its utilization on acoustic signals and gravitational wave signals has been discussed in earlier sections, we hint at the scope and versatility of our methods here with a possible usage in neuroscience, as a sample of other future
<figure><img src="image_11.png" /><figcaption>**Figure 10:** Extracted volumes for the cello signal. (a) Minimum-volume cycles. (b) Boundary of stable volumes. The number indicates the relevance of the component in the persistent diagram. (c) Persistence diagram.</figcaption></figure>

applications to non-acoustic domains. This application consists in the characterization of TF *signa-tures* using TDA and the tools described before to extract *features* that can be fed to downstream classification (or regression) models.

To demonstrate the potential of such a framework for the analysis of time-series, we describe an example based on non-acoustic signals from a intracranial electroencephalograms (EEGs) dataset [64, see datasets C, D and E]. These EEGs can be classified as having a seizure present, i.e. pathological signals, or reflecting seizure-free brain activity. Seizures have been associated with a more synchronized activity in EEGs, whereas typical brain activity could be succinctly described as noise-like [64]. Because of these differences, when a seizure is present, EEGs tend to be more organized and *rhythmic*, which translates into a different signature in the TF plane consisting of more persistent holes in the pattern of spectrogram zeros than seizure-free EEGs. Figure 9 shows the difference between two signals obtained from the dataset and their corresponding PDs. Figure 9a shows the spectrogram corresponding to a seizure-free EEG, whereas Fig. 9b corresponds to an EEG with a seizure. Comparing the persistence diagrams from Figs. 9a and 9b, the latter has more persistent holes, linked to a pattern of harmonic-like components in the spectrogram.

## 8.2 Discussion and future work

To our knowledge, this work is arguably the first to systematically apply TDA as a broad-based methodological toolbox for spectrogram zero-based signal processing. The limited amount of existing literature that appear to be related to our work often involves specialised applications and particular kinds of signals. For instance, [5] investigates the problem of audio obfuscation, where topological features based on Betti numbers are invoked. In a different direction, [65] proposes a groupoid $C^*$ algebra-based method to combine dimensionality reduction and generic persistent homology for recovering topological features of the input signal, while [66] works on filtered level sets of the spectrogram for feature extraction. As alluded before, persistence diagrams have also been converted to a functional statistic called the *accumulated persistence function* in spatial statistics [8].

While the APF can be used in conjunction with classical spatial statistical tests, our approach is more tailored to time-frequency detection and reconstruction, which we have numerically demonstrated (cf. the detection power of APF in 4).

The signal detection methods described in this paper, as well as the detection tests introduced in [20, 30, 24], fall under a different paradigm than traditional signal detection methods where a model of the signal is available, like matching filter techniques and variations thereof [60, 67]. In this sense, the signal detection test defined in Sec. 4 constitutes a signal-agnostic alternative, relying on a model of the noise and Monte Carlo simulations rather than a parametric model of the signal. Such tests might find future applications to the detection of so-called *continuous* gravitational waves, for which matching possible signal templates is computationally ineffective [68]. Furthermore, extending the proposal to the continuous wavelet transform and its zeros [21] could further improve detection performance in
certain applications. Connections with recent work on hypothesis testing based on the PD [12] are to
be explored in future work, where applicability of certain *universal* null models could be considered in
order to increase the power of our tests.

9 Conclusion

In this work, we propose a very general approach to the detection and reconstruction of signals from noisy measurements, leveraging a robust interplay between time-frequency analysis based on spectrogram zeros and state-of-the-art concepts and techniques from topological data analysis (TDA). Our approach is undergirded by a solid statistical foundation based on the theory and methods of multiple hypothesis testing, resulting in the possibility of providing scientifically principled predictions augmented with confidence guarantees. In terms of ideas, our method interfaces with a range of disparate scientific domains – time-frequency analysis (which is a staple tool in the processing of a wide variety of time-dependent signals), hyperuniformity of the spectrogram zero patterns (which is of interest in statistical physics), the theory of stable and minimum volumes and that of persistence trees (of recent interest in topological data analysis), and conceptual ingredients like FWER and Bonferroni correction (which are of fundamental importance in statistical inference).

In terms of applications, we are able to bring our methods to bear on a highly versatile class of
time-varying signals, including acoustic data from musical instruments (of key interest in the field of
audio and speech recognition), Gravitational Wave (GW) data from the celebrated LIGO experiment
(a burgeoning area of interest in physics and astronomy), and futuristic applications towards EEG time
series data (of interest in medical science). Our topologically-based methods appear to demonstrate
good performance with acoustic data especially in low SNR regimes, and allow for the possibility of
signal reconstruction in highly challenging problems such as gravitational waves (where state-of-the-art
methods focus largely on the question of signal detection).

We demonstrate the possibility of invoking the powerful machinery of TDA for very general classes
of signal processing problems, without the requirement of the a-priori presence of prominent geometrical
structures in the signals (which is a feature of most classical TDA applications). We believe that this
would provide engineers with a new and powerful suite of signal processing tools on one hand, and
topological data analysts with a novel and interesting class of problems on the other. Our work opens
up the possibility of integrating topology-inspired methods with dedicated state-of-the-art techniques
for specific signal processing tasks (such as neural network-based approaches for GW) for improved
outcomes, and suggests natural extensions to other time-varying signal classes in statistics, medicine
and data science.

Acknowledgment

The authors wish to thank Curtis Condon, Ken White, and Al Feng of the Beckman Institute of
the University of Illinois for the bat data used in Fig. 3 and for permission to use it in this paper,
and Alvin Chua of the National University of Singapore for helpful discussions regarding gravitational
wave data. JMM and RB were supported by ERC grant BLACKJACK ERC-2019-STG-851866 and
ANR grant BACCARAT ANR-20-CHIA-0002. KAT was supported by the NUS Research Scholarship.
SSM was partially supported by the INSPIRE research grant DST/INSPIRE/04/2018/002193 from
the Dept. of Science and Technology, Govt. of India, and a Start-Up Grant from Indian Statistical
Institute. SG was supported in part by the Singapore MOE grants R-146-000-250-133, R146-000-312-
114, A-8002014-00-00 and MOE-T2EP20121-0013.
# A Tools from topological data analysis

This section gives an in-depth description of the main tools we take from topological data analysis and that we use throughout the main document, namely persistent homology, persistent diagram, minimum volumes and stable volumes.

## A.1 Persistent homology

We extract a few definitions from [43], to which we refer for more details. Let $S \subset \mathbb{R}^2$ be a finite set of points, and $r > 0$. Let $B_u(r)$ be the closed ball of center $u \in \mathbb{R}^2$ and radius $r$. We consider

$$ \mathcal{A}(r) = \left\{ \sigma \subseteq S : \bigcap_{u \in \sigma} (B_u(r) \cap V_u) \neq \emptyset \right\}, $$

where $V_u$ is the (closed) cell associated to $u$ in the Voronoi diagram of $S$. In other words, $V_u$ is the set of points in $\mathbb{R}^2$ that are closer to $u$ than to any other point in $S$.

For $r > 0$, $\mathcal{A}(r)$ is called the $\alpha$-complex of $S$ of scale $r$ [43, Section III.4]. The $\alpha$-complex is an example of a simplicial complex; Other common choices include the Delaunay and Vietoris-Rips complexes [43, Chapter III]. Elements of any $\mathcal{A}(r)$ are subsets of $S$, and are called simplices. A simplex formed by $p+1$ affinely independent points is called a $p$-simplex, and is identified with its convex hull: a 0-simplex is thought of as a point, a 1-simplex as a segment with positive length, a 2-simplex as a non-flat triangle. 3-simplices and above do not exist in dimension 2. For simplicity, we assume henceforth that no three points of $S$ are collinear, so that all triplets define a non-degenerate triangle.

Figure 11a illustrates $\mathcal{A}(r)$ for different values of $r$, with $S$ the set of zeros of a spectrogram. Momentarily considering $r$ as the age of the construction, the basic intuition is that as $r$ grows, holes appear and die in $\mathcal{A}(r)$, and that holes that stay alive for a long time are likely to be indicative of some underlying structure, in our case a signal. Persistent homology is a set of tools to define and keep track of these holes.

Define $\mathcal{A}(0) = S$. When $r \ll 1$, $\mathcal{A}(r) = \mathcal{A}(0)$. As $r$ grows, $\mathcal{A}(r)$ takes on a finite number of values, progressively including all subsets of $S$ at some finite time $r_{max} > 0$. In short, for any $0 < r < s$,

$$ S = \mathcal{A}(0) \subset \mathcal{A}(r) \subsetneq \mathcal{A}(s) \cdots \subset \mathcal{A}(r_{\max}) = \mathcal{P}(S). \qquad (\text{A-1}) $$

The increasing sequence of sets $\mathcal{A}(r)$, $r \ge 0$, is called a filtration. Since each complex differs from the previous one by one simplex, we can equivalently describe the filtration A-1 by an ordered list of simplices $((r_1, \sigma_1), \dots, (r_J, \sigma_J))$ for some $J \ge 1$, where $r_i$ is the value of $r$ for which $\sigma_i$ is added to the filtration.

Denote a formal weighted sum $c = \sum a_i \sigma_i$ of $p$-simplices $\sigma_1, \dots, \sigma_k$, with $a_1, \dots, a_k \in \mathbb{Z}/2\mathbb{Z}$, as a $p$-chain. Let $r > 0$. The $p$-chains of $\mathcal{A}(r)$, with the term-by-term sum modulo 2, form a group$^2$ $C_p(r)$. A natural morphism is then the boundary operator $\partial_{p+1} : C_{p+1}(r) \to C_p(r)$, defined as follows. For a $(p+1)$-simplex $\sigma$, define its faces as all its subsets of cardinality $p$. The boundary $\partial_{p+1}\sigma$ is then the formal sum of all the faces of $\sigma$, with constant weight 1. We then extend $\partial_{p+1}$ to $C_p(r)$ by linearity, as

$$ \partial_{p+1} \sum a_i \sigma_i := \sum a_i \partial_{p+1} \sigma_i. $$

The kernel of $\partial_p$ is the subgroup $Z_p(r) < C_p(r)$ of so-called $p$-cycles, and the image of $\partial_{p+1}$ is the subgroup $B_p(r) < C_p(r)$ of $p$-boundaries. A fundamental result is that $B_p(r) < Z_p(r)$, i.e., that a $p$-boundary is a $p$-cycle.

While introducing algebraic notions might seem an unnecessary overhead at first sight, it crucially allows the translation of topological questions into algebraic manipulations that are amenable to computer implementation. In particular, a key notion is the $p$-th homology group $H_p(r) = Z_p(r)/B_p(r)$. In words, an element of $H_p(r)$ (i.e., an equivalence "class" in the quotient) represents a set of $p$-chains obtained by taking a $p$-cycle of $H_p(r)$, and adding to it any $p$-boundary in $H_p(r)$. For illustration, Figure 12 gives three cycles in $Z_1(r)$ for $r$ large enough. The two cycles in red, which go around the same hole, are in the same equivalence class: one can indeed write the "larger" red cycle as the sum of the "smaller" one and the boundary "in between" the two cycles. On the contrary, one cannot

<footnote>$^2$We reserve sans serif characters for groups.</footnote>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11:** (a) Simplicial complexes in a filtration for three values of *r*. Red points, black segments, and shaded triangles are, correspondingly, 0-simpleces, 1-simpleces and 2-simpleces. For each point **u**, the intersection $B_u(r)$ and $V_u$ is shown in green, where $B_u(r)$ is a closed ball centered at **u** with radius *r*, and $V_u$ is the Voronoi cell associated to **u**; see Sec. A.1. (b) Spectrogram of a mixture of a two-component signal and real white Gaussian noise. The zeros are shown as white dots. (c) Minimum-volume cycles (continuous lines) and a non-minimum cycle (dashed line) superimposed to the log-spectrogram of a signal. (d) and (e) show persistence diagrams $D_0$ (red dots) and $D_1$ (blue dots). Each dot correspond to a birth-death pair. (d) Persistence diagrams for a mixture of signal and noise, the spectrogram of which can be seen in the right panel. For $D_1$, the two squared birth-death pairs are significantly further away from the diagonal, corresponding to long-lived *holes* caused by the two-component signal. (e) PD for white Gaussian noise.</figcaption></figure>
obtain the magenta cycle by summing any of the two red cycles and a boundary: the magenta cycle corresponds to a different class in $H_1(r)$. Intuitively, equivalence classes in $H_p(r)$ represent the holes we are trying to detect: for $p = 1$, a class contains all loops that encircle a given hole.

Now, as $r$ grows, $H_p(r)$ gains some classes and loses some, and the persistence diagram is a way of keeping track of these events. Formally, let $r_0 > 0$. For $c \in H_p(r_0)$, denote by $x(c) \in [0, r_0]$ the smallest $r > 0$ for which $c \in H_p(r)$, i.e., $x(c)$ is the birth time of the class. Similarly, let $y(c)$ be the smallest $r > 0$ for which $c \notin H_p(r)$, i.e., the death time of the class. The set of all birth-death time pairs

$$ \mathcal{D}_p = \{(x(c), y(c)) : c \in H_p(r) \text{ for some } r > 0\} \subset \mathbb{R}^2 $$

is called the $p$-th persistence diagram of $S$. Plotting the diagrams $\mathcal{D}_0$ and $\mathcal{D}_1$, as in Figure 11, gives topological information on the pattern $S \subset \mathbb{R}^2$. Figure 11 shows the 0- and 1-persistence diagrams corresponding to a real white Gaussian noise realization, whereas Figure 11d shows the diagrams corresponding to the signal with two components used in Figure 11a. On the one hand, $\mathcal{D}_0$ (red dots in Figure 11) shows the birth and death times of connected components, all born at $r = 0$, and progressively dying until only one remains indefinitely. On the other hand, $\mathcal{D}_1$ (blue dots in Figure 11) gives a visual diagnostic of the number and size of holes in the pattern $S$. Points far away from the diagonal indicate long-living holes, which we interpret as being indicative of a signal: in Figure 11e two birth-death pairs are significantly further away from the diagonal than the remaining pairs. Such pairs are typically representative of signal components. In contrast, Figure 11d shows that the pairs from $\mathcal{D}_1$ stay close to the diagonal (notice the difference in the scale between the axes in Figure 11e and Figure 11d.

In Section 4 of the main text, we design a detection test based on the persistence diagram, while in Section 5 we investigate the possibility to identify the support of the signal once the test has rejected the hypothesis of pure noise. For the latter task, we need one more notion from persistence homology.

## A.2 Minimum volumes

Once an anomalous point $(x(c), y(c))$ has been detected in the persistence diagram $\mathcal{D}_1$, the signal processor would like to find a particular cycle corresponding to $c$, and use it to identify a candidate support for the signal that hopefully led to detection. Since one can add any boundary without changing $c$, it is intuitively appealing to search for a member of $c$ that encloses the smallest possible portion of $\mathbb{R}^2$, that portion likely corresponding to the signal support. One way to formalize this minimizing class representer is the notion of minimum-volume cycle, introduced in [47] and generalized in [45].

While the original definition takes the form of the solution to a constrained optimization problem that has to be numerically approximated, our special interest in the alpha-complex and the case of $S \subset \mathbb{R}^2$ allows us to easily obtain minimum-volume cycles through the persistence tree (PT) $G = (\mathcal{V}, E)$ of $S$ defined by Algorithm 1; the PT construction was introduced by [47].

The input to Algorithm 1 is the list of simplices ordered by appearance in the filtration $((r_1, \sigma_1), \dots, (r_J, \sigma_J))$. The set of nodes $\mathcal{V}$ is initialized with $\sigma_\infty = (\mathbb{R}^2 \cup \{\infty\}) \setminus \mathcal{P}(S)$, which represent the complement of the simplicial complex in the one-point compactification of $\mathbb{R}^2$ [45]. Each node/edge also contains extra information in $\mathbb{R}$ by saving the value of $r_i$ corresponding to $\sigma_i$, with $r_\infty = \infty$. Algorithm 1 simply moves through the filtration in inverse order, starting from $\sigma_J$. When $\sigma_j$ is a 2-simplex, it is automatically added to $\mathcal{V}$.

If $\sigma_j$ is a 1-simplex, then is added to the list of edges $E$ as long as it joins two unconnected nodes from $\mathcal{V}$. The Root($\sigma_s, \mathcal{V}, E$) procedure traverses the tree $(\mathcal{V}, E)$ upstream, starting from $\sigma_s$, searching for a node without parents, i.e. the root corresponding to $\sigma_s$. Note that the notation $\sigma \xrightarrow{r_j} \sigma_k$ indicates the way the persistence tree should be read, meaning here that $\sigma_k$ is a descendant of $\sigma$, and the edge $(r_j, \sigma_j)$ joining both nodes is a 1-simplex $\sigma_j$ appearing at $r = r_j$ in the filtration. [45] shows that $(b, d) \in \mathcal{D}_1(r)$ if and only if there exists $\sigma \xrightarrow{r_k} \sigma_d$, i.e. an edge with $r = r_b$ pointing at a node $\sigma_d$. The persistence diagram can thus be read from $G$.

Figure 12a shows a simple example of the simplicial complexes obtained for some values of $r$. The corresponding persistence tree is depicted in Fig. 12b. In the tree, the 1-simplex $\sigma_b = [0, 4]$, corresponding to the birth of a non-trivial cycle at $r = 0.097$, points to the 2-simplex $\sigma_d = [0, 1, 4]$, the appearance of which marks the death of this homological class at $r = 0.152$. Formally, a minimum-volume cycle generating the pair $(b, d) \in \mathcal{D}_1$ is defined as $\partial_2 MV(\sigma_d)$, where $MV(\sigma_d)$ is the volume
<caption>**Algorithm 1** Building the persistence tree [47, 45]</caption>

**Input:** A filtration as an ordered list $((r_1, \sigma_1), \dots, (r_J, \sigma_J))$.

**Output:** The persistence tree $G = (\mathcal{V}, E)$.

1: $\mathcal{V} = \{(r_\infty, \sigma_\infty)\}, E = \emptyset$

2: **for** $j = J, \dots, 1$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\sigma_j$ is a 2-simplex **then**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Add node $(r_j, \sigma_j)$ to $\mathcal{V}$

5: &nbsp;&nbsp;&nbsp;&nbsp;**else if** $\sigma_j$ is a 1-simplex **then**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Get $\sigma_s$ and $\sigma_t$, the two 2-simplices sharing $\sigma_k$ as an edge.

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\sigma_{s'} = \text{Root}(\sigma_s, \mathcal{V}, E)$$

8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\sigma_{t'} = \text{Root}(\sigma_t, \mathcal{V}, E)$$

9: **end if**

10: **if** $s' == t'$ **then**

11: &nbsp;&nbsp;&nbsp;&nbsp;Continue

12: **else if** $s' > t'$ **then**

13: &nbsp;&nbsp;&nbsp;&nbsp;Add edge $(\sigma_{s'} \xrightarrow{r_j} \sigma_{t'})$ to $E$.

14: **else if** $s' < t'$ **then**

15: &nbsp;&nbsp;&nbsp;&nbsp;Add edge $(\sigma_{t'} \xrightarrow{r_j} \sigma_{s'})$ to $E$.

16: **end if**

17: **end for**

obtained as

$$MV(\sigma_d) = \{\sigma_d\} \cup \text{des}(\sigma_d, E), \qquad (A-2)$$

$\text{des}(\sigma_d, E)$ being the descendants of the node $\sigma_d$ in the persistence tree $G$ (dashed blue line in the PT of Fig. 12b). As shown in [45], $MV(\sigma_d)$ is the volume with the smallest size, measured in number of 2-simplices, for the component associated with a given $(b, d) \in \mathcal{D}_1$. Hence, we will refer to (A-2) as a minimum volume (MV).

Fig. 12b, the MV corresponding to $\sigma_d = [0, 1, 4]$ is displayed in blue in the bottom section of the right panel. This notion corresponds to the intuition conveyed by Figure 12a: at the iteration where the class died, all the volume is then “filled”. We can recover the MV by tracking the 2-simplices in it in the reverse order of appearance in the filtration, starting from $\sigma_d$ and adding one 2-simplex $\sigma_i$ at the time for $r_i > r_b$. This is exactly what we do by taking the descendants of $\sigma_d$ from the persistence tree in (A-2).

Returning now to an example based on the spectrogram, Fig. 11c illustrates these concepts by superimposing different cycles on top of a spectrogram. The continuous red and magenta lines show the minimum-volume cycles corresponding to the significant pairs $(x(c), y(c))$ highlighted in Fig. 11d (the color indicates the corresponding pair in the persistence diagram). Both red lines, continuous and dashed, show two cycles in the same class, i.e. represented by the same pair in the persistence diagram. However, the dashed cycle is not a minimum-volume cycle.

Finally, note that if one takes two points in the persistence diagram, nothing prevents the two corresponding minimum volumes to intersect. In our particular application of identifying signal components in the time-frequency plane, we prefer restricting ourselves to non-overlapping volumes. To achieve that, we simply restrict to minimum volumes (A-2) such that $\sigma_d$ is a child of the root $\sigma_\infty$. By construction, none of these minimum volumes intersect.

# B Stable volumes

We noticed that a small change in the input signal, like the difference between two noise realizations, can produce large variations of the cycles –and their volumes–, even for the cycles associated to long-living components [48, 49, 46]. Such instability in the presence of noise is a known issue in computational topology, and some approaches to stabilize the components of interest have been proposed [49, 46].

Intuitively, one way to robustify the construction of minimum volumes is to take the intersection of minimum volumes corresponding to slightly altered inputs, where the alteration is small enough that we can identify a volume from one signal to the other. Formally, we focus here on the so-called stable
<figure><img src="image_14.png" /><figcaption>**Figure 12:** (a) Alpha-complexes in a filtration, for different values of $r$. The minimum-volume cycle born at $r = 0.097$ is shown in blue, dashed line. (b) Persistence tree corresponding to the filtration (top), where $[-1, -1, -1]$ represents $\sigma_\infty$, and the minimum volume (bottom, in blue) corresponding to the class born at $r = 0.097$. The blue, dashed line in the persistence tree links the 2-simplices comprised in the minimum volume. (c) and (d) show minimum volumes and stable volumes, respectively, for comparison. (e) Shows the volume size of the stable volumes for different values of $\epsilon$.</figcaption></figure>
volumes recently proposed by Obayashi [46], which can be defined directly from the persistence tree in Sec. A.2.

More precisely, given a filtration ((r₁, σ₁), ..., (rJ, σJ)), a pair (b, d) ∈ D₁ and its corresponding persistence tree, the stable volume with noise bandwidth ϵ ≥ 0 is defined as

$$SV_{\epsilon}(\sigma_{b}, \sigma_{d}) = \{\sigma_{d}\} \cup \left( \bigcup_{\sigma \in C_{\epsilon}(\sigma_{b}, \sigma_{d})} \mathrm{des}(\sigma, G) \right), \quad (A-3)$$

where C_ϵ(σ_b, σ_d) = {σ ∈ X⁽²⁾ | σ_d r_k σ and r_k ≥ (1 + ϵ) r_b}, X⁽²⁾ is the set of all 2-simplices. Obayashi [46] shows that (A-3) is indeed the intersection of minimum volumes corresponding to a certain sophisticated type of perturbation of the filtration, which is a trade-off between actually representing noise in the input and allowing the identification of minimum volumes across perturbations.

By comparing (A-3) to (A-2), it can be interpreted that in order to get a stable volume, one needs first to prune the branches stemming from the node σ_d in G for r < (1 + ϵ)r_b in the filtration. Then the stable volume is computed by aggregating the remaining descendants of σ_d, in a similar fashion as the minimum volumes defined in Sec. A.2.

Figures 12c and 12d illustrate the difference between minimum and stable volumes, obtained for the same signal. By comparing both cases, one can see that stable volumes (Fig. 12d) produce a more accurate estimation of the region of interest, in this case the domain corresponding to each of the four impulses in the signal.

The bandwidth parameter ϵ plays a key role in stabilizing a minimum-volume cycle, since ϵ = 0 makes (A-3) equivalent to (A-2). Using ϵ > 1.0 produces a stable volume the size of which –in number of 2-simplices– is smaller than the minimum-volume cycle obtained for the same point in the PD.

Figure 12e shows the monotonically decreasing relationship between the size of the volume and ϵ for a signal with SNR = 0 dB, i.e. noise and signal with the same energy.

Theoretically, ϵ depends on the SNR of the signal [46]. However, an explicit relationship between the bandwidth parameter and the level of noise in the signal is not available in this case. Hence, in practice, ϵ should be large enough to reduce the instabilities, but not as large that the resulting volume is too small to correctly represent the signal domain of interest. The next section describes a practical approach to tune ϵ.

## B.1 Tuning the bandwidth parameter

As suggested by Obayashi in [46], a possible heuristic to fix the bandwidth parameter consists in counting the number of 2-simplices in the volume obtained for a range of values of ϵ in order to build a curve like the one shown in Fig. 12e. A plateau in the curve indicates a stabilization of the volume, hence a value of ϵ in a plateau of this curve should be selected. For each component extracted from the PD, we look for a plateau in an interval [0, ϵ_max], and the minimum value of ϵ in the plateau is chosen. An ϵ_max = 0.15 was empirically chosen.

It might occur that there are not enough points in the interval [0, ϵ_max] to form a plateau. When there are less than ten points in the search interval [0, ϵ_max], the latter is expanded until a plateau can be computed. This approach is used to automatically fix ϵ when stable volumes are computed, as in Fig. 12d.

## B.2 Signal Reconstruction

Algorithm 2 summarizes the steps for signal estimation based on minimum (or stable) volumes.

The performance of the reconstruction for the synthetic signals was measured using the quality reconstruction factor (QRF, see [19]) given by QRF = 10 log₁₀ (||f||₂²/||f - f̃||₂²) (dB), where f is the original noiseless signal and f̃ is a denoised approximation. The QRF can be understood as an estimation of the SNR of the output signal after a reconstruction approach such as Algorithm 2. In practice, it should be higher than the SNR of the noisy input signal to indicate an improvement. Notice that one can only compute the QRF if the noiseless reference is available, hence the importance of working with synthetic signals at an early stage.
**Algorithm 2** Signal reconstruction

**Input:** Noisy signal $h$, significance level $\alpha$.

1. Compute the STFT $V_g(h)$ using (1).

2. Find the zeros of $|V_g(h)|^2$.

3. Compute the persistence diagram based on the zeros.

4. Estimate the number of significant components to extract from the PD using the tests described in Sec. 5.2 with significance level $\alpha$.

5. Compute the volume of the minimum-volumes cycle (or stable volume) associated with each relevant component.

6. Approximate the signal domain $D$ as the union of the computed volumes.

7. Compute a signal estimation $\tilde{f}$ using (4) of the main document.

8. return $\tilde{f}$.

## C Hypothesis testing

This section complements the contents in the main document, providing more technical details regarding the statistical tools we use for signal detection and denoising.

### C.1 Simple hypothesis

Given a noise model $\xi$ and a candidate signal $f$, signal detection can be stated as a binary hypothesis test [37, 20]:

$$ \begin{cases} H_0 : h = \xi, & \text{i.e. the observation is pure noise} \\ H_1 : h = f + \xi, & \text{i.e. the observation is a mixture of signal plus noise} \end{cases} \quad (A-4) $$

Formally, both $H_0$ and $H_1$ here can be identified with singletons containing one probability measure each. Alternately, if no candidate signal $f$ is known, one can define $H_1$ as the complement $H_0^C$ of $H_1$. A test aiming to detect a signal can be characterized by its

1. Specificity: the probability $1 - \alpha$ of correctly keeping $H_0$, where $\alpha$ is the so-called significance of the test, i.e. the probability of wrongfully rejecting $H_0$.

2. Sensitivity: the probability of correctly rejecting $H_0$, also known as the statistical power of the test.

The main idea behind our detection tests is computing the top order statistics $\Xi_{(\ell)}$, $\ell = 1, 2, \dots$, for a generic test statistic $\Xi$. In our case, $\Xi$ could be the distance to the diagonal in the PD, or the energy associated to the minimum (or stable) volume of that component in the spectrogram. One should be able to detect the presence of a signal by comparing the quantiles of the top order statistics against the same quantiles under the null model (e.g., white Gaussian noise). For practical purposes, quantiles are obtained by Monte Carlo simulations. This procedure is explained below.

From the underlying noise model, we generate $B$ independent samples. Let $\Xi_{(1),1}, \dots, \Xi_{(1),B}$ denote the corresponding top order statistics. Then a $p$-value for testing $H_0$ can simply be obtained as

$$ \frac{1}{B} \sum_{j=1}^{B} \mathbb{I}(\Xi_{(1),j} > \Xi_{(1),\text{obs}}), $$

where $\Xi_{(1),\text{obs}}$ is the top order statistics obtained from the sample to be tested. A signal is then detected if the $p$-value is below a predefined value of $\alpha \in (0, 1)$. By a symmetry argument, if $B\alpha$ is an integer, the significance of the test is guaranteed to be at least $\alpha$.

### C.2 Multiple hypothesis testing

Our applications on signal detection and reconstruction require inferring the existence of one or more signal components. This is addressed by using multiple hypothesis testing, conducting several tests simultaneously (for signal detection) or sequentially (for estimating the number of signal components).
### C.2.1 Family-wise error rate

A particularly important aspect of multiple testing is to construct procedures capable of controlling a chosen type I error-related criterion. A traditional choice is the Family-wise Error Rate (FWER), which is the worst-case probability of incorrectly rejecting any null-hypothesis. Formally, this can be defined as follows; see e.g. [52]. Let $(\mathcal{X}, \mathcal{X})$ be a mesurable space, and $\mathcal{P}$ the set of probability measures on $\mathcal{X}$. A hypothesis $H$ can be seen as a subset of $\mathcal{P}$, and $H$ is true under $P$ if $P \in H$, and false if $P \notin H$. Denote the set of true hypotheses under $P$ as $\mathcal{T}(P) = \{H \in \mathcal{H} : P \in H\}$. Then, writing $\mathcal{R}$ for the (observation-dependent) subset of $\mathcal{P}$ that is rejected, the FWER is defined as

$$ \text{FWER}(\mathcal{R}) = \sup_{P \in \mathcal{P}} P(\mathcal{R} \cap \mathcal{T}(P) \neq \emptyset). $$

### C.2.2 Simultaneous tests

With the notations of Section C.1, we may want to base our decision on several top order statistics $\Xi_{(k)}$ instead of a single one. This leads to multiple hypothesis testing, which we now describe more formally. Start with a suitable choice of $K$ most persistent components to evaluate, i.e. the components associated with the $K$ longest distances to the diagonal in the PD. Denote by $\mu_k^{(0)}$, $1 \le k \le K$, the marginal distribution of $\Xi_{(k)}$ under our null hypothesis $H_0$ that the original signal is pure white Gaussian noise. Then, we conduct $K$ tests simultaneously for the $K$ hypotheses

$$ H_{0,k} : \Xi_{(k)} \sim \mu_k^{(0)}, \quad k = 1, \dots, K. $$

Formally, $H_{0,k}$ is the set of distributions over signals under which the marginal of $\Xi_{(k)}$ is $\mu_k^{(0)}$.

Again, from the $B$-many simulated persistence diagrams, with statistics $\Xi_{(k),j}$, $1 \le j \le B$, we compute a series of $p$-values $p_1, p_2, \dots, p_K$ as

$$ p_k := \frac{1}{B} \sum_{j=1}^{B} \mathbb{I}(\Xi_{(k),j} > \Xi_{(k),\text{obs}}). $$

For each $H_{0,k}$, we perform a level-$\alpha_k$ test, for some $\alpha_k \in (0, 1)$. We define the hypothesis of pure noise to be rejected if and only if $H_{0,k}$ is rejected for at least one $1 \le k \le K$. For any distribution $P$ under which data and the $B$ simulations are i.i.d.,

$$ P(H_{0,k} \text{ is rejected for at least one } k) = P\left(\bigcup_{k=1}^{K} \{p_k < \alpha_k\}\right) < \sum_{k=1}^{K} \alpha_k. \quad (A-5) $$

Choosing

$$ \alpha_k = \alpha/K \quad (A-6) $$

for some predefined $\alpha \in (0, 1)$, (A-5) and thus the FWER are bounded by $\alpha$. The choice of (A-6) is called Bonferroni's correction, and is often the default choice in practice among so-called $\alpha$-corrections.

As an aside, note that multiple hypothesis testing can be used to estimate the number of holes, e.g. using as estimator the cardinality of $\{1 \le k \le K : p_k < \alpha_k\}$. However, we found this estimator not to be robust to the presence of spurious holes in noise, and rather recommend sequential testing.

### C.2.3 Sequential tests

In Section C.2.2, we claimed a rejection when at least one—not necessarily the first—of the hypotheses $H_{0,1}, \dots, H_{0,K}$ was violated. Instead, we could check whether $p_1 < \alpha_1$, and if not then check whether $p_2 < \alpha_2$, and continue so until either $H_{0,k}$ is rejected or we reach the final $K$th test and do not reject. This leads to another rejection region, i.e. another test, which we call a sequential test. This also naturally suggests another estimator for the number of holes. Indeed, intuitively, rejecting $H_{0,k}$ for $k = 1, \dots, k_*$ is an indication of the presence of at least $k_*$ holes in the configuration of zeros. Formally, we set

$$ \hat{n}_{\text{holes}} = \min\{0 \le k \le K-1 : p_{k+1} \ge \alpha_{k+1}\}, $$

with the convention that $\min\emptyset = 0$. The power of the sequential test, as well as the distribution of the estimator $\hat{n}_{\text{holes}}$, depend on the choice of $(\alpha_k)$.
### C.2.4 Choice of levels of significance in sequential tests

For any distribution $P$ under which data and the $B$ simulations are i.i.d., the test wrongly detects if and only if $\hat{n}_{holes} > 0$, so that

$$P(\text{false detection}) = P(\hat{n}_{holes} > 0) = P(p_1 \ge \alpha_1) = \alpha_1.$$

So, as long as $\alpha_1 \le \alpha$, any $\alpha$-correction controls the FWER of the sequential test. We consider three different choices of $\alpha$-corrections $(\alpha_k)_{1 \le k \le K}$ in the main paper: the Bonferroni correction we used in simultaneous testing (A-6), as well as

1. Polynomial decay: $\alpha_k = \frac{\alpha}{k^m}$, $m \ge 0$.

2. Geometric decay: $\alpha_k = \alpha\beta^k$, $\beta \in (0, 1)$.

Note that for the sequential test, the Bonferroni correction is overly conservative.

Finally, on top of controlling the FWER, we note that fixing an $\alpha$-correction allows to further elaborate on the distribution of the estimator $\hat{n}_{holes}$. For instance, let $k^* > 0$ and assume that $P \in H_{0,1}^C \cap H_{0,2}^C \cap \dots \cap H_{0,k^*}^C \cap H_{0,k^*+1}^C$. Intuitively, this translates to the fact that $k^*$ is the underlying number of holes in the signal. Let us compute the probability under such a $P$ to misestimate $k^*$.

For $j < k^*$,

$$P(\hat{n}_{holes} = j) = P(p_1 < \alpha_1, p_2 < \alpha_2, \dots, p_{j+1} \ge \alpha_{j+1}) = \epsilon_{j+1} \prod_{l=1}^{j} (1 - \epsilon_l),$$

where $\epsilon_l = P(p_l \ge \alpha_l | p_j < \alpha_j \text{ for } j < l)$. In particular,

$$P(\hat{n}_{holes} < k^*) = \sum_{j=0}^{k^*-1} \left( \epsilon_{j+1} \prod_{l=1}^{j} (1 - \epsilon_l) \right). \qquad (A-7)$$

On the other side, for $j > k^*$,

$$P(\hat{n}_{holes} = j) = P(p_1 < \alpha_1, p_2 < \alpha_2, \dots, p_{j+1} \ge \alpha_{j+1}) = \left( \prod_{l=1}^{k^*} (1 - \epsilon_l) \right) \cdot (1 - \tilde{\alpha}_{j+1}) \prod_{l=k^*+1}^{j} \tilde{\alpha}_l,$$

where $\tilde{\alpha}_l = P(p_l < \alpha_l | p_j < \alpha_j \text{ for } j < l)$. In particular,

$$P(\hat{n}_{holes} > k^*) = \left( \prod_{l=1}^{k^*} (1 - \epsilon_l) \right) \cdot \sum_{j>k^*} \left( (1 - \tilde{\alpha}_{j+1}) \prod_{l=k^*+1}^{j} \tilde{\alpha}_l \right) \qquad (A-8)$$

Under the polynomial decay $\alpha$-correction, $\alpha_k = \frac{\alpha}{k^m}$, $m \ge 0$. We assume that for each $l$, the probability of event $\{p_k < \alpha_k\}$ is offset by some small error $\delta_k$ after conditioning on $\{p_k < \alpha_k \text{ for } j < k\}$, i.e. $\tilde{\alpha}_k = \frac{\alpha}{k^m} + \delta_k$. Then

$$\begin{aligned}
\sum_{j>k^*} \left( (1 - \tilde{\alpha}_{j+1}) \prod_{l=k^*+1}^{j} \tilde{\alpha}_l \right) &= \sum_{j>k^*} \left( (1 - \alpha_{j+1}) \prod_{l=k^*+1}^{j} \alpha_l \right) + o(\tilde{\delta}) \\
&= \sum_{j>k^*} \left( \alpha^{j-k^*} \cdot \left( \frac{k^*!}{j!} \right)^m \cdot \left( 1 - \frac{\alpha}{(j+1)^m} \right) \right) + o(\tilde{\delta}) \\
&= \sum_{j>k^*} \left( \alpha^{j-k^*} \left( \frac{k^*!}{j!} \right)^m - \alpha^{j+1-k^*} \left( \frac{k^*!}{(j+1)!} \right)^m \right) + o(\tilde{\delta}) \\
&= \alpha \cdot \left( \frac{k^*!}{(k^*+1)!} \right)^m + o(\tilde{\delta}) \\
&= \frac{\alpha}{(k^*+1)^m} + o(\tilde{\delta}).
\end{aligned} \qquad (A-9)$$
Here $o(\tilde{\delta})$ consists of terms that approach 0 as $\delta_k \to 0$ for all $k$. Therefore, the probability of incorrectly estimating the number of components is

$$P(\hat{n}_{\text{holes}} \neq k^*) = \sum_{j=0}^{k^*-1} \left( \epsilon_{j+1} \prod_{l=1}^{j} (1 - \epsilon_l) \right) + \frac{\alpha}{(k^*+1)^m} \left( \prod_{l=1}^{k^*} (1 - \epsilon_l) \right) + o(\tilde{\delta}). \quad (\text{from (A-7), (A-8), (A-9)})$$

Similarly, under the geometric decay $\alpha$-correction $\alpha_k = \alpha\beta^k$, $\beta \in (0, 1)$, assuming that $\tilde{\alpha}_k = \alpha\beta^k + \delta_k$, we have

$$P(\hat{n}_{\text{holes}} \neq k^*) = \sum_{j=0}^{k^*-1} \left( \epsilon_{j+1} \prod_{l=1}^{j} (1 - \epsilon_l) \right) + \alpha\beta^{k^*+1} \left( \prod_{l=1}^{k^*} (1 - \epsilon_l) \right) + o(\tilde{\delta}).$$

In either $\alpha$-correction method, the probability of inaccurately estimating the number of holes depends on the $\epsilon_j$'s and $\delta_j$'s, which in turn depend on the choice of test statistics $\Xi_{(j)}$ and their distributions $\mu_j^{(0)}$. A good choice of test statistics $\Xi_{(j)}$ would be able to distinguish between noise and true signal, ensuring $\epsilon_j$ remains small for all $j$. Ideally, if the dependence among the order statistics $\Xi_{(j)}$ (for $1 \le k \le K$) is negligible, then $\delta_j$ will also be negligible. Coupled with a suitable choice of $\alpha, m$ and $\beta$, both polynomial decay and geometric decay $\alpha$-corrections will estimate the exact number of holes correctly with high probability.

## D Detection based on the accumulated persistence function

A common variant of the persistence diagram is the *rotated-rescaled persistence diagram* (RRPD). A birth-death pair $(b_i, d_i)$ appears in the RRPD diagram as coordinates $(m_i, l_i)$, where $m_i = (b_i + d_i)/2$ and $l_i = d_i - b_i$, i.e. the half-life and total life of a component. The Accumulated Persistence Function (APF) [8] is defined as a function of the half-life

$$\text{APF}: m \mapsto \sum_{i} l_i \mathbb{I}(m_i \le m), \qquad (10)$$

where we assume that the RRPD is computed from $D_0$ or $D_1$, $m \in [0, m_{\text{max}}]$. Note that APF is a monotonically increasing function. Following [53, 54], a test can be devised by first choosing a significance level $\alpha$, simulating $m$ signals coming from the null distribution, computing

$$t_{\ell} = \max_{m \in [0, m_{\text{max}}]} (\text{APF}_{\ell}(m) - \overline{\text{APF}}(m)), \quad \ell = 0, \dots, L, \qquad (11)$$

where $\text{APF}_0$ is the APF obtained from the observed data, $\text{APF}_{\ell}$ the APF obtained from the $\ell$th simulated data ($1 \le \ell \le L$), and $\overline{\text{APF}}$ is the average of all $\text{APF}_{\ell}$, $\ell = 0, \dots, L$. Finally, the latter are then sorted in ascending order $t_{(\ell)}$, $\ell = 1, \dots, L$, and the null hypothesis is rejected if $t_0 \ge t_{(\ell^*)}$, where $\ell^* = \alpha(L+1)$. Notice that $\alpha$ and $L$ should be chosen so that $\ell^*$ is an integer. By a symmetry argument, the resulting test has significance $\alpha$ [53].

## E Detection of gravitational waves

Using persistent homology to detect gravitational waves is particularly intriguing, as it differs significantly from applications involving synthetic signals with background white noise. In this section, we provide more details regarding the GW application described in the main document.

### E.1 The Kaggle Gravitational Waves Detection dataset

For illustration purposes, we use the gravitational wave data from the Kaggle Gravitational Waves Detection competition [56]. The dataset consists of realistic simulations of gravitational wave measurements from a network of three gravitational wave interferometers (LIGO Hanford, LIGO Livingston, and VIRGO). Each data sample, or *event*, contains three time series, each 2 seconds long and sampled at 2048 samples per second. The data is labeled as “0” if there is only noise and “1” if there is a noisy gravitational wave signal. The detection task is challenging, particularly for data with very low SNR.
### E.1.1 Data subsets

In order to characterize in depth the performance of the signal detection tests for gravitational waves, we divided the original dataset in different groups according to their SNR. Since an estimation of the SNR of each signal is not available, we used an adaptive hard threshold [57, 58] on the spectrogram to segment the dataset first, following the notion that when the signature of a gravitational wave in the spectrogram is clear its values will surpass a higher threshold than when the signal is deeply buried in noise. Given that the real and imaginary parts of the STFT of white Gaussian noise $\xi$ with variance $\sigma^2$ are normally distributed, and provided that the analysis window in (1) is normalized to $\|g\|_2^2 = 1$, the spectrogram ordinates follow an exponential distribution $\text{Exp}(\lambda) : t \mapsto \frac{1}{\lambda}e^{\frac{t}{\lambda}}$ with parameter $\lambda = \sigma^2$ [59, 39].

We define the normalized spectrogram maximum (NSM), $\hat{\zeta} = \hat{\gamma}/\hat{\sigma}$, where $\hat{\gamma} = \max_{u,v} |V_g(f)(u, v)|$ is the maximum modulus of the STFT across the time-frequency plane, and $\hat{\sigma}$ is the standard deviation of the noise, estimated from the signal using a median absolute deviation estimator [57, 58]

$$ \hat{\sigma} = \frac{\text{median}(|\Re\{V_g(f)\}|)}{0.6745}. $$

The NSM serves as a practical proxy for the SNR; signals with higher NSM values are expected to have higher SNR values. We then consider the set of GW signals as a union of the subsets $Z_i$ for $i \in \{2, 3, 4, 5, 6, 7\}$, where each $Z_i$ is defined to contain the signals with $\hat{\zeta} \ge i$. Note that $Z_i \subset Z_j$ if $i \ge j$.

We collected 50,000 signal samples for $Z_2$, and computed $\hat{\zeta}$ for each one, using the signal corresponding to LIGO Hanford. From this group, $|Z_3| = 49,995$, $|Z_4| = 14,542$, $|Z_5| = 1,528$, $|Z_6| = 670$, and $|Z_7| = 255$. Examples of spectrograms of the signals in the different subsets can be found in Sec. E of the SI document. As can be noted in Figure 13a, most signals in $Z_7$ are visually identifiable in their spectrograms, while $Z_5$ and $Z_6$ contain numerous more subtle signals. Note that our NSM thresholding is only used to create subsets of the data for comparison purposes; it plays no role in the detection procedure whatsoever.

To compute our ROC curves, for each value of $\hat{\zeta}$, we randomly selected an equal number of noise-only and signal-plus-noise samples in each subset of the data. This makes sure we have class-balanced datasets for each value of $\hat{\zeta}$.

Finally, note that the signal components in GW signals typically span a narrow range of frequencies for more than half of the signal duration, causing them to persist for shorter periods in the persistence tree. Hence, energy-based test statistics are expected to enhance detection performance, compared to distance-based tests. Detection was carried out by applying each test to the three interferometers and claiming a detection when the null was rejected for at least one of the three.
<figure><img src="image_15.png" /><figcaption>(a) Random subsets of spectrograms from the GW dataset with thresholds $\hat{\zeta} \in \{5, 6, 7\}$. The signals correspond to the LIGO Hanford observatory.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13:** Gravitational wave (GW) detection. (a) shows spectrograms of random subsets of signals coming from the dataset with thresholds $\hat{\zeta} \in \{5, 6, 7\}$. All the signals have a GW present, even though is not evident in most of the signals from the sets with $\hat{\zeta} = 5$ and $\hat{\zeta} = 6$ (b) Gravitational wave detection ROC curves. Test statistic: energy (Bonferroni).</figcaption></figure>
References

[1] R. Nagaoka, K. Masuzawa, M. Taniwaki, A. L. Foggiatto, T. Yamazaki, I. Obayashi, Y. Hiraoka, C. Mitsumata, and M. Kotsugi. Quantification of the coercivity factor in soft magnetic materials at different frequencies using topological data analysis. *IEEE Transactions on Magnetics*, 2024.

[2] A. Cipriani, C. Hirsch, and M. Vittorietti. Topology-based goodness-of-fit tests for sliced spatial data. *Computational Statistics & Data Analysis*, 179:107655, 2023.

[3] M. Nicolau, A. J. Levine, and G. Carlsson. Topology-based data analysis identifies a subgroup of breast cancers with a unique mutational profile and excellent survival. *Proceedings of the National Academy of Sciences*, 108(17):7265–7270, 2011.

[4] K. Xia and G.-W. Wei. Persistent homology analysis of protein structure, flexibility, and folding. *International journal for numerical methods in biomedical engineering*, 30(8):814–844, 2014.

[5] W. Reise, X. Fernández, M. Dominguez, H. A. Harrington, and M. Beguerisse-Diaz. Topological fingerprints for audio identification. *SIAM Journal on Mathematics of Data Science*, 6(3):815–841, 2024.

[6] Sohail Iqbal, H Fareed Ahmed, Talha Qaiser, Muhammad Imran Qureshi, and Nasir Rajpoot. Classification of COVID-19 via homology of CT-SCAN. *arXiv preprint arXiv:2102.10593*, 2021.

[7] Gloria Colombo, Ryan John A Cubero, Lida Kanari, Alessandro Venturino, Rouven Schulz, Martina Scolamiero, Jens Agerberg, Hansruedi Mathys, Li-Huei Tsai, Wojciech Chachólski, et al. A tool for mapping microglial morphology, morphOMICs, reveals brain-region and sex-dependent phenotypes. *Nature neuroscience*, 25(10):1379–1393, 2022.

[8] Christophe AN Biscio and Jesper Møller. The accumulated persistence function, a new useful functional summary statistic for topological data analysis, with a view to brain artery trees and spatial point process applications. *Journal of Computational and Graphical Statistics*, 28(3):671–681, 2019.

[9] Helen M Byrne, Heather A Harrington, Ruth Muschel, Gesine Reinert, Bernadette J Stolz, and Ulrike Tillmann. Topological methods for characterising spatial networks: A case study in tumour vasculature, 2019.

[10] Michael Kerber. Persistent homology: state of the art and challenges. *International Mathematische Nachrichten*, 231(15-33):1, 2016.

[11] Gunnar Carlsson. Topological methods for data modelling. *Nature Reviews Physics*, 2(12):697–708, 2020.

[12] Omer Bobrowski and Primoz Skraba. A universal null-distribution for topological data analysis. *Scientific reports*, 13(1):12274, 2023.

[13] LIGO Scientific Collaboration, Virgo Collaboration, et al. A guide to LIGO–Virgo detector noise and extraction of transient gravitational-wave signals. *Preprint*, 2020.

[14] Zhen Zeng, Rachneet Kaur, Suchetha Siddagangappa, Tucker Balch, and Manuela Veloso. From pixels to predictions: Spectrogram and vision transformer for better time series forecasting. In *Proceedings of the Fourth ACM International Conference on AI in Finance*, ICAIF '23, page 82–90, New York, NY, USA, 2023. Association for Computing Machinery.

[15] Luyuan Xie. TRLS: A time series representation learning framework via spectrogram for medical signal processing. In *2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, IEEE Signal Processing Society SigPort, 2024.

[16] Bairui Du, Delmiro Fernandez-Reyes, and Paolo Barucca. Image processing tools for financial time series classification, 2020.

[17] Hau-Tieng Wu. Current state of nonlinear-type time-frequency analysis and applications to high-frequency biomedical signals. *Current Opinion in Systems Biology*, 23:8–21, 2020.
[18] Patrick Flandrin. Time-frequency filtering based on spectrogram zeros. *IEEE Signal Processing Letters*, 22(11):2137-2141, 2015.

[19] Sylvain Meignen, Thomas Oberlin, Philippe Depalle, Patrick Flandrin, and Stephen McLaughlin. Adaptive multimode signal reconstruction from time-frequency representations. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 374(2065):20150205, 2016.

[20] Rémi Bardenet, Julien Flamant, and Pierre Chainais. On the zeros of the spectrogram of white noise. *Applied and Computational Harmonic Analysis*, 48(2):682-705, 2020.

[21] Rémi Bardenet and Adrien Hardy. Time-frequency transforms of white noises and Gaussian analytic functions. *Applied and computational harmonic analysis*, 50:73-104, 2021.

[22] Günther Koliander, Luis Daniel Abreu, Antti Haimi, and José Luis Romero. Filtering the continuous wavelet transform using hyperbolic triangulations. In *2019 13th International conference on Sampling Theory and Applications (SampTA)*, pages 1-4. IEEE, 2019.

[23] Luis Alberto Escudero, Naomi Feldheim, Günther Koliander, and José Luis Romero. Efficient computation of the zeros of the Bargmann transform under additive white noise. *Foundations of Computational Mathematics*, 24(1):279-312, 2024.

[24] Subhroshekhar Ghosh, Meixia Lin, and Dongfang Sun. Signal analysis via the stochastic geometry of spectrogram level sets. *IEEE Transactions on Signal Processing*, 70:1104-1117, 2022.

[25] Antti Haimi, Günther Koliander, and José Luis Romero. Zeros of Gaussian Weyl-Heisenberg functions and hyperuniformity of charge. *Journal of Statistical Physics*, 187(3):22, 2022.

[26] Juan M Miramont, François Auger, Marcelo A Colominas, Nils Laurent, and Sylvain Meignen. Unsupervised classification of the spectrogram zeros with an application to signal detection and denoising. *Signal Processing*, page 109250, 2023.

[27] Ali Moukadem, Barbara Pascal, Jean-Baptiste Courbot, and Nicolas Juillet. The analytic Stockwell transform and its zeros. *arXiv preprint arXiv:2407.17076*, 2024.

[28] Barbara Pascal and Rémi Bardenet. Point processes and spatial statistics in time-frequency analysis. *arXiv preprint arXiv:2402.19172*, 2024.

[29] Timothy J Gardner and Marcelo O Magnasco. Sparse time-frequency representations. *Proceedings of the National Academy of Sciences*, 103(16):6094-6099, 2006.

[30] B. Pascal and R. Bardenet. A covariant, discrete time-frequency representation tailored for zero-based signal detection. *IEEE Transactions on Signal Processing*, 2022.

[31] Juan M. Miramont, Rémi Bardenet, Pierre Chainais, and François Auger. Benchmarking multi-component signal processing methods in the time-frequency plane, 2024.

[32] Salvatore Torquato. Hyperuniformity and its generalizations. *Phys. Rev. E*, 94:022122, Aug 2016.

[33] Luis Daniel Abreu, Joao M Pereira, José Luis Romero, and Salvatore Torquato. The Weyl-Heisenberg ensemble: hyperuniformity and higher Landau levels. *Journal of Statistical Mechanics: Theory and Experiment*, 2017(4):043103, 2017.

[34] Marco Salvalaglio, Dominic J. Skinner, Jörn Dunkel, and Axel Voigt. Persistent homology and topological statistics of hyperuniform point clouds. *Phys. Rev. Res.*, 6:023107, May 2024.

[35] Charles Emmett Maher, Yang Jiao, and Salvatore Torquato. Hyperuniformity of maximally random jammed packings of hyperspheres across spatial dimensions. *Phys. Rev. E*, 108:064602, Dec 2023.

[36] Abel H G Milor and Marco Salvalaglio. Inferring traits of hyperuniformity from local structures via persistent homology. *Journal of Physics: Condensed Matter*, 37(14):145401, feb 2025.
[37] Patrick Flandrin. *Time-frequency/time-scale analysis*, volume 10. Academic press, 1998.

[38] Karlheinz Gröchenig. *Foundations of time-frequency analysis*. Springer Science & Business Media, 2001.

[39] Patrick Flandrin. *Explorations in time-frequency analysis*. Cambridge University Press, 2018.

[40] G. Ascensi and J. Bruna. Model space results for the Gabor and Wavelet transforms. *IEEE Transactions on Information Theory*, 55(5):2250–2259, 2009.

[41] W. Chachólski, T. Dyckerhoff, J. Greenlees, G. Stevenson, D. Herbera, W. Pitsch, and S. Zarzuela. *Building Bridges Between Algebra and Topology*. Advanced Courses in Mathematics - CRM Barcelona. Springer International Publishing, 2018.

[42] F. Chazal and B. Michel. An introduction to topological data analysis: fundamental and practical aspects for data scientists. *Frontiers in artificial intelligence*, 4:667963, 2021.

[43] H. Edelsbrunner and J. L. Harer. *Computational topology: an introduction*. American Mathematical Society, 2022.

[44] T. K. Dey and Y. Wang. *Computational topology for data analysis*. Cambridge University Press, 2022.

[45] I. Obayashi. Volume-optimal cycle: Tightest representative cycle of a generator in persistent homology. *SIAM Journal on Applied Algebra and Geometry*, 2(4):508–534, 2018.

[46] Ippei Obayashi. Stable volumes for persistent homology. *Journal of Applied and Computational Topology*, pages 1–36, 2023.

[47] B. Schweinhart. *Statistical Topology of Embedded Graphs*. PhD thesis, Princeton University, 2015.

[48] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer. Stability of persistence diagrams. In *Proceedings of the twenty-first annual symposium on Computational geometry*, pages 263–271, 2005.

[49] Paul Bendich, Peter Bubenik, and Alexander Wagner. Stabilizing the unstable output of persistent homology computations. *Journal of Applied and Computational Topology*, 4(2):309–338, 2020.

[50] Larry Wasserman. *All of statistics: a concise course in statistical inference*. Springer Science & Business Media, 2013.

[51] Sandrine Dudoit, Mark J Van Der Laan, and Mark J van der Laan. *Multiple testing procedures with applications to genomics*, volume 10. Springer, 2008.

[52] Magalie Fromont, Matthieu Lerasle, and Patricia Reynaud-Bouret. Family-wise separation rates for multiple testing. *The Annals of Statistics*, 44(6):2533–2563, 2016.

[53] Adrian Baddeley, Peter J Diggle, Andrew Hardegen, Thomas Lawrence, Robin K Milne, and Gopalan Nair. On tests of spatial pattern based on simulation envelopes. *Ecological Monographs*, 84(3):477–489, 2014.

[54] Mari Myllymäki, Tomáš Mrkvička, Pavel Grabarnik, Henri Seijo, and Ute Hahn. Global envelope tests for spatial processes. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 79(2):381–404, 2017.

[55] The GUDHI Project. *GUDHI User and Reference Manual*. GUDHI Editorial Board, 2015.

[56] Chris Messenger, Christopher Zerafa, Elena Cuoco, Michael J. Williams, and Walter Reade. G2Net gravitational wave detection, https://kaggle.com/competitions/g2net-gravitational-wave-detection, 2021.

[57] David L Donoho and Jain M Johnstone. Ideal spatial adaptation by wavelet shrinkage. *Biometrika*, 81(3):425–455, 1994.
[58] Stephane Mallat. *A Wavelet Tour of Signal Processing: The Sparse Way*. Academic Press, 2008.

[59] Duong-Hung Pham and Sylvain Meignen. A novel thresholding technique for the denoising of multicomponent signals. In *2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, pages 4004–4008. IEEE, 2018.

[60] Eric Chassande-Mottin and Patrick Flandrin. On the time-frequency detection of chirps. *Applied and Computational Harmonic Analysis*, 6(2):252–281, 1999.

[61] ITU-T Recommendation. Perceptual evaluation of speech quality (PESQ): An objective method for end-to-end speech quality assessment of narrow-band telephone networks and speech codecs. *Rec. ITU-T P. 862*, 2001.

[62] Rafael G. Dantas Miao Wang, Christoph Boeddeker and Ananda Seelan. Pesq (perceptual evaluation of speech quality) wrapper for python users, May 2022.

[63] Matteo Torcoli, Thorsten Kastner, and Jürgen Herre. Objective measures of perceptual audio quality reviewed: An evaluation of their application domain dependence. *IEEE/ACM Transactions on Audio, Speech, and Language Processing*, 29:1530–1541, 2021.

[64] Ralph G Andrzejak, Klaus Lehnertz, Florian Mormann, Christoph Rieke, Peter David, and Christian E Elger. Indications of nonlinear deterministic and finite-dimensional structures in time series of brain electrical activity: Dependence on recording region and brain state. *Physical Review E*, 64(6):061907, 2001.

[65] Mijail Guillemard and Armin Iske. On groupoid C*-algebras, persistent homology and time-frequency analysis. *Preprint*, 105, 2011.

[66] Holger Boche, Mijail Guillemard, Gitta Kutyniok, and Friedrich Philipp. Signal analysis with frame theory and persistent homology. In *The Conference of Sampling Theory and Applications, SampTA'13*, volume 1, page 1, 2013.

[67] Anthony D Whalen. *Detection of signals in noise*. Academic press, 2013.

[68] Anjana Ashok, Banafsheh Beheshtipour, Maria Alessandra Papa, Paulo CC Freire, Benjamin Steltner, Bernd Machenschalk, Oliver Behnke, Bruce Allen, and Reinhard Prix. New searches for continuous gravitational waves from seven fast pulsars. *The Astrophysical Journal*, 923(1):85, 2021.