A novel multi-layer modular approach for real-time fuzzy-identification of gravitational-wave signals

The first level exploits an approach derived from speech processing techniques [43–45]. The very first step is considering a continuous signal in smaller, overlapped, short data frames of predefined length (miniframes, hereafter). After preliminary tests, we settled on a optimal miniframe size of $\omega_{f} = 100$ ms, with an overlap of $\sigma_{f} = 50$ ms. Each miniframe is directly processed by the first level through a specific model, which we call ground model. The ground model is a miniframe fuzzy classifier, which is trained to distinguish between noise-only miniframes (null hypothesis, negative output of the classifier) and miniframes containing a BBH signal (alternative hypothesis, positive output of the classifier), exploiting numerical features calculated from the data of a given miniframe. The ground model is trained to return an output closer to 1 if a GW is present inside the miniframe, and closer to 0 if that is a noise-only pattern. A more detailed discussion of the ground model, including all details of its derivation, is provided in section 3.

The second level of the framework does not have direct access to the original data segment coming from the continuous stream. Instead, the signal is basically translated into a sequence of numerical values provided by the first layer. The stream of output values passed by the first level to the second level requires a dedicated interpretation. In this paper, we propose two possible level 2 layers (modules, in the modular philosophy), based on two possible interpretation models of level 1 output. The separation of the ground model derivation from its interpretation is therefore a multi-modeling procedure (see section 4.1)). A strong edge of such an approach is the absence of an a priory hypothesis about the GW transient position within the data segment.

Level 2 analyzes the sequence returned by the first layer in fixed chunks of n_f frames (local analysis). Its purpose is the search for a GW transient specific pattern in the space of values of the first layer. We call this step moving window analysis, because, at each iteration, the n_f miniframes selection moves forward by a defined number of miniframes, usually $\sigma_{mw} = 1$. This layer is crucial to the classification because it looks at the signal features in a wider time span, combining information on adjacent frames. The output of level 2 is a stream of fuzzy-like or strong integer predictions: 1 stands for a GW transient, 0 for noise. Ultimately, level 3 issues an interferometer GW alert by analyzing, in a similar manner than that used by level 2, the output of the second layer. A possible implementation of level 3 could, for example, identify clusters of ones after thresholding the fuzzy-like values of level 2. It is reasonable that clusters of ones will represent a GW in the input signal, as well clusters of zeros will mark a noise-only sequence in the original timeseries.

The framework described above is an early implementation, but promising in the continuous analysis of data, as it will be shown in section 4. Once trained, each framework level performs low complexity operations, thanks to the simplicity of the underlying models. The ground model used in the proposed implementation leads to a level-2 accuracy slightly lower than existing state-of-the-art CNN approach, such as [31, 33], but it is outstandingly computationally simple. This makes the entire framework particularly suitable to be implemented in real-time, low-latency, searches even on a single machine, thus avoiding parallel and distributed calculations and related complications.

As previously mentioned, the database used to train and test the model contains simulated interferometer timeseries, either noise-prevailing GW signals or noise-only data samples. The duration of the simulated segments was fixed to 1 s, as previously done, for example, in [30, 31, 33]. This is sufficiently large to contain a typical BBH GW signal (often called strain), for which the relevant part that enters the interferometer frequency band lasts usually a few tens of second; for their short duration, BBH signals are usually called transients. Figure 2 shows a schematic example of time-series used to derive the model. The two key steps required to create GW samples include the generation of a template BBH waveform and its sum to an interferometer noise to match with a target signal-to-noise ratio (SNR). In this way, two types of data segments are produced: noise only and signal + noise samples.

Zoom In Zoom Out Reset image size

BBH template waveforms, i.e. BBH GW signals according to existing theoretical models, are generated through the LALsuite library [46], a package for gravitational and relativistic astrophysics, which is used to perform a numerical simulation of BBH coalescence events (see appendix A).

Because the simulated source is located at a fixed distance, for a given choice of its astrophysical parameters and a given assumption on the detector noise, signal + noise segments would have exclusively a certain SNR ¹⁰ . However, to derive and test a model capable to detect GWs signals, it is crucial to probe its capabilities at different SNRs. For example, testing the accuracy of a model with different SNRs enables to probe the generalization capabilities of the framework and allows to evaluate the rise in performance with more favorable SNRs. To add detector noise to the time series, data segments containing noise only are generated according to a typical aLIGO design sensitivity curve. Noise was chosen to be exclusively Gaussian. Stationary noise is a simplified hypothesis to facilitate the comparison with previous approaches already published. A number of generated noise segments were instead used to populate the dataset of noise only samples. Finally, before summing the generated GW signals to the stationary noise, the waveform amplitude was suitably scaled to match with the target optimal SNR. The latter is defined according to the following equation [31]:

$$\begin{equation} \rho^2_{opt} = \langle u,u \rangle \end{equation} \tag{ 1 }$$

where u is the generated template and the inner product $\langle a,b \rangle$ is the noise-weighted cross-correlation defined in [23]:

$$\begin{align*}\langle a,b\rangle = 4 Re\left[\int_{f_{min}}^\infty \frac{\tilde{a}^\star(f\,)\tilde{b}(f\,)}{S_n(f\,)} \right]\end{align*}$$

being $\tilde{a}(f\,)$ ($\tilde{b}(f\,)$) the frequency-domain representation of the signal a (b) and S_n the single-sided detector noise PSD; f_min is the frequency of the gravitational-wave template at the beginning of the time-series.

Before splitting the generated segments into the datasets used to derive the models, the time series are suitably whitened in order to flatten the frequency distribution dominated by the detector noise. This is a common practice in the analysis of GW data. The stationary noise and the smoothness of the simulated aLIGO noise spectrum allowed to use a short, 0.5 s, padding. To this end, the segments had an additional 0.5 s of data before and after each generated 1 s time series. Signal content, after whitening, was truncated using a Tukey window ($\alpha = 1/8$) to its central 1 s of data. In addition, similarly to previous CNNs studies and to allow a meaningful comparison of the local results, the peak amplitude of each waveform within the final time series, was randomly placed within the range [0.75 s, 0.95 s].

For each model derived in this paper (namely the ground model, which serves the proposed level 1, and two interpretation models used in two possible versions of level 2), we build a training, a cross-validation and a testing dataset. The cross-validation dataset was used to make decisions on the complexity of the model, while the final accuracy is always estimated using testing data. The default training and cross-validation datasets consisted of 1000 noise-only data segments and 1000 signal + noise data segments at SNR 8 (i.e. $\rho_{opt} = 8$ according to equation (1)). Finally, for the testing datasets, SNR was varied from 6 to 14 in integer steps.

The frequency band of interest for GW searches with second generation interferometers ranges usually from about 10 Hz to several hundreds Hz. This is in good superposition with the typical frequency band explored in audio signal analyses. For example, the evolution of a GW transient has often been compared to a bird chirp. The general idea of our framework is thus inspired by one-dimensional signal processing techniques usually applied for the task of speech recognition, i.e. the analysis of audio signals for the purpose of recognizing voice from a persistent background noise. The underlying techniques are well-developed in the literature and are particularly powerful to classify voiced signals even in the presence of a particularly loud background. The mechanism behind the framework is that of processing data samples through smaller miniframes. In the field of audio analysis and speech-recognition, it is common use to consider the sound on a short-term basis, in order to extract information from the whole audio sequence. Therefore, an audio signal is usually divided into small segments and then the short-term processing stage is carried out. The short-term processing consists in the extraction of physical and perceptual features from each segment. Physical features are measurements computed directly from the sound wave, such as spectrum, energy, entropy, cepstral coefficients, zero crossing rate, and so on [47]. Perceptual features are values related to the perception of sounds by the human hearing, such as rhythm, pitch, timbre and noisiness. The feature extraction provides the numerical values that describe properties of an audio segment, so the real information is derived from the analysis of such values, depending on the desired task.

The literature written upon speech processing and audio analysis algorithms is well established [48], and many open-source libraries are already available for use. In this paper, we use a number of existing libraries for speech recognition and audio analysis to compute the short term features used to classify GW interferometer data segments. A full mathematical description of each feature used to derive the ground model is beyond the scope of this paper. Instead, a short resume of the most important libraries and features used for the ground model is reported in table 1. The total number of features extracted from the libraries is 184. In addition to the feature extraction process on the original, whitened, time series, we also added a pre-processing stage for the input waveform, which uses a wavelet decomposition as a low-pass filter. To this end, the wavelet decomposition of a signal is used to partially filter the original waveform. The reconstructed waveform depends on the specific wavelet used at first instance, and the spectrum of the reconstructed signal, consequently, is different from the original one. It usually appears smoothed for higher frequencies, while the most energetic tracks in the low-frequency band are highlighted from the surrounding noise. Because we used 9 wavelets, 9 different pre-processed versions of the timeseries are considered. After filtering the timeseries with the wavelet filter, the basic set of features is again extracted for each pre-processed waveform. The total number of features for each miniframe finally was $184 + 184 \cdot 9 = 1840$.

In terms of performance, the full feature set extraction on 2k s of signal took 240 s on average (the ratio is greater than 8:1), using a commercial 64-core Threadripper CPU and a simple parallel processing workaround meant to take advantage of the large number of cores. It is important to stress that this result was obtained without a dedicated (low-level) optimization of the feature extraction codes. This raw comparison suggests that the audio feature extraction, upon which the whole framework relies, can ideally keep up to the data buffer required in a real-time situation. In addition, after the model is derived, it is no longer needed to compute all the features described in this section, because the ground model relies on a strongly reduced subset of features, as it will be discussed in detail in section 3.4. In fact, one of the crucial aspects of the method used to derive the ground model in the present implementation is the so-called feature selection. Feature selection is the capability of a ML algorithm to suitably select a subset of features, among those used for model training, whose informative content is not significantly enhanced when additional features are included. Consequently, one can derive a model that exploits only a reduced number of features, even without significantly affecting the overall accuracy. Usually, a suitable trade-off between complexity (i.e. number of features or their computational complexity) and accuracy is selected. This gives the unique opportunity to explore a large dimensionality feature space, which includes a large body of features derived from speech processing literature, without enhancing the complexity of the resulting model. The latter is a requirement for the ground model of a real time detection algorithm for GWs. To this end, the ground model proposed in this paper is derived using the BP, a state-of-the-art tool for the formal modeling of data based on a hybridization of genetic programming and artificial neural networks [34, 35]. This software is particularly suitable for the task of deriving a model aimed to detect GWs in real-time because it has outstanding capabilities of feature selection, as proven for example in [37, 38].

The BP is a novel tool for the formal modeling of data based on a hybridization of genetic programming and artificial neural networks [34, 35]. Genetic programming falls in the field of evolutionary computing (EC). In computer science, EC is a branch of artificial intelligence that often deals with optimization problems through algorithms inspired by the Darwinian evolutionary theory [49]. Within this field, BP deals with the so-called symbolic regression, i.e. the process during which some data are fitted by means of a suitable analytical expression derived by the algorithm itself. Many novel algorithms for symbolic regression tasks have recently shown that genetic programming is a particularly powerful technique to solve even extremely complex problems [50]. Furthermore, some genetic programming implementations are also capable of feature selection, i.e. to suitably exclude part of the input variables, thus reducing the dimensionality of the feature space.

In the framework of evolutionary computation, a solution of a given optimization problem is initially encoded according to a predefined scheme. Each solution is usually called individual, while a set of individuals forms a population. A numerical value, called fitness, is then associated to each individual in the population. The fitness function quantifies how much a solution is a good solution for the target problem. For instance, in the problem of symbolic regression, i.e. in the task to derive a model for the description of some data, the fitness function might account for the prediction error of the model and/or for its computational complexity. The average fitness in the population is usually maximized (or either minimized, depending on the definition of the fitness function) applying some suitable evolutionary criteria, inspired by the natural selection.

As previously mentioned, in the implementation of BP, genetic programming and NNs cooperate with the aim of deriving a proper analytical model $g({\mathbf{x}})$ to fit a set of previously labeled input patterns $\{({\mathbf{x}}_i, {\mathbf{y}}_i)\ : 1\unicode{x2A7D} i\unicode{x2A7D} N_p\}$. The evolutionary computing approach followed in BP is part of a line of genetic programming research that foresees the evolution of tree structures representing formal mathematical expressions. Roughly speaking, to maximize the fitness function, the following steps are executed:

• (i)  
  an initial population of randomly constructed trees, representing solutions of the modeling problem, is generated;
• (ii)  
  each tree is evaluated through the fitness function and a fitness value is assigned to all individuals in the population;
• (iii)  
  a new population of trees is created, starting from the population at the previous step, through some evolutionary operators, e.g. copy, crossover, mutation, selection and heuristic operators;
• (iv)  
  return to step 2, until the termination criterion is met.

In the end, the best individual, i.e. the tree with the maximum fitness, is returned as the output result.

The evolutionary computing core of BP operates a global search for the optimal solution of the modeling problem. Such a global search is extremely powerful to identify a solution close to a global maximum of the fitness function, but is usually rather slow to converge to the maximum itself. NNs are therefore used to serve as a hill-climbing operator, i.e. an operator capable to perform a fast local search for the maximum of the fitness function by suitably varying the constants in the mathematical formula. This is a quite time consuming task, thus it is executed only if an individual is particularly promising (i.e. possibly close to a maximum of the fitness function). In the neural part of BP, each expression is initially transformed into a multi-layer feed-forward NN by replacing operators in the original expression with specialized neurons, while all constants are treated as the weighs of the NNs. The mathematical expression of the error gradient in the weight space is then formally calculated and an enhanced gradient descending technique is used to update the weights. Learning weights are usually different for each parameter to optimize and are dynamically varied during the training.

The fitness function used in BP takes into account the prediction error and the number of features exploited by the model (see [34] for additional details).

The software is written in highly optimized C and runs on Linux 64bit environments. It has also a parallel and distributed implementation. Each learning task is executed by more clients controlled by a process called brain_server and running on a server machine. Moreover, BP has also the capability, through a dedicated evolutionary computing algorithm, to self-adapt almost all the internal parameters (hyperparameters) required for a learning task. This makes BP a double evolutionary tool [35], capable of automatically self-tuning itself for the particular problem to solve, thus optimizing the overall performance through, each time, a dedicated set of hyperparameters. BP has been successfully used so far in a number of research fields, see e.g. [36–39], but this is the first time it is used in the field of GW physics.

Let us summarize the main points discussed so far. The detection of GW signals is based on the analysis of noisy timeseries, in which the dominant frequency range is typically in the audible band. Supported by this fact, the proposed framework is inspired by audio analysis techniques, and specifically by speech processing. One of the key aspects of the newly proposed framework is its layered implementation, which allows to easily specialize the framework for the desired real-time analysis pipeline, even combining different existing algorithms. The lowest level of the framework relies on the initial predictions operated by the ground model. The latter are performed for each individual miniframe, in which the data segment is subdivided, and the successive local analysis (level 2) is performed by analyzing the ground model output for a continuous stream of miniframes. In this paper, we propose a level 1 module whose ground model is derived through BP. BP implements a few key advantages for its use in the derivation of real-time GW search algorithms:

• BP is capable of a particularly advanced feature selection. For example, the ground model discussed in this paper was initially derived exploiting a full set of 1840 audio analysis features, but, at the end of the training process, it foresees exclusively 4 features. The strongly reduced number of features allows to minimize the computational complexity of the model, thus making it ideal for real-time applications, even if the model is derived using a large body of information on a high-dimensional feature space. In addition, feature selection can also serve as a guideline towards the correct choice of features suitable for the study of GW signals, a crucial information even for the derivation of future algorithms.
• The resulting model has a rather simple analytical expression based on numerical features, thus enabling a fast and easy implementation of the model even in low-cost systems based on commercial CPUs.

In the mechanism described in section 2, every data sample of 1 s is processed in 19 miniframes of 100 ms length each, as described in section 3.1. Given that a typical learning dataset contained 2000 data samples ¹¹ , BP had in turn $N_p = 19 \cdot 2000 = 38k$ patterns available for the training and cross-validation tasks. Each pattern consisted of 1840 inputs (the features $\{x_i\}$) and a categorical target output value, either 0, for miniframes containing noise only, or 1, if the miniframe is in a region where a GW signal is present. For a similar problem, in a high-dimensionality feature space in the presence of large noise, BP usually achieved a convergence within a time of a few hours using a single Threadripper CPU. With a similar procedure as described for example in [36, 38], after a few preliminary runs, we settled to an optimal trade-off between complexity and accuracy. An analogous procedure was also applied to set the optimal miniframe length to 100 ms; in this respect, we performed preliminary learning runs and tested the deriving preliminary models on an additional subset of the cross-validation data. The performance of the model, for each selected trade-off, was monitored by analyzing the prediction error on the cross-validation dataset, which contained the same number of patterns as in the training dataset.

The analytical implementation of the ultimately derived model is reported in table 2, together with a brief description of the corresponding numerical features. Only 4 features are exploited by the model, and it is therefore not needed to compute all the features described in table 1. The output value, $y_i = f(x^i_1, \ldots, x^i_{1840})$ is usually not defined in the range [0, 1]. In order to fulfill the condition $f({\mathbf{x}})\in[0,1]$ for all possible input vectors x, it is convenient to use the corresponding constrained function $g({\mathbf{x}}) = \max(\min(f({\mathbf{x}}), 1),0)$.

Table 2. The analytical formulation ultimately derived by BP for the miniframe output value y. Despite the complexity of the problem and the large dimensionality of the feature space used to train the model, the derived formula foresees exclusively four features and has a particularly low computational complexity. This model is used to perform a ground prediction on each miniframe, individually.

$y = - 0.071\,33 + 0.862\,74 \cdot pyAA_{31} \cdot \chi + 70.473\,28 \cdot pyAA_{33}$
where $\quad c = {\max\left ({0.224\,73}+{{\ln\left(db17_{ess3}\right)}\cdot {\textrm{erf}\left(db32_{bb4}\right)}}-{{2.6911}\cdot {pyAA_{31}}}, \; db17_{ess3} \right)}$
$\chi = \exp \left[ \min \left( 0.036\,685, db17_{ess3} \right) - {1.0623 \cdot {db17_{ess3}}} + c \; \right]$
the following 4 features have been selected:
pyAA31	The second last element of the chroma vector a given by the pyAA library.
pyAA33	The standard deviation of the 12 elements of the chroma vector
	given by the pyAA library.
$db17_{ess3}$	The energy of the waveform in band ]150 800] Hz given by Essentia library
	on the db17 pre-processed waveform.
$db32_{bb4}$	The energy of the waveform in band ]200 300] Hz given by Essentia library
	on the db32 pre-processed waveform.

^a The chroma vector is a 12-element representation of the spectral energy, where the bins represent the 12 equal-tempered pitch classes of music (in the western-type semitone spacing).

One of the key aspects of the proposed method is that the derived model has an analytical formulation, entailing only a few nodes. This allows to perform a prediction with a reduced computational complexity compared to CNN approaches. In addition, the layered framework requires suitable interpretation at each level of decision, making the framework particularly versatile for many real-time applications.

From a detailed inspection of the ground model equation of table 2, one can point out some interesting points. The features exploited by the model are exclusively: pyAA₃₁, pyAA₃₃, $db17_{ess3}$, $db32_{bb4}$. This indicates that their informative content, for the task of identifying GW signals in noisy and short time intervals, is equivalent to the informative content of the 1840 features initially selected for the learning task and reported in table 1. In other words, forcing the model to exploit all the features of table 1 would increase the complexity of the model without a statistically significant reduction of the prediction error. Selected features are briefly described in table 2. The first two features, extracted from the library pyAA, are related, respectively, to the energy of the strain in the frequency band from about 98 Hz to about 415 Hz, and to the overall variation of energy (expressed as the standard deviation) in the frequency domain from 30 Hz to 1000 Hz. The latter might usually exhibit lower values if the miniframe is strongly dominated by noise, as the noise power is equal at all frequencies after whitening the data segment, but will tend towards larger values if the miniframe contains a GW signal, whose power lies, for a given miniframe, in a narrow frequency region. The second two features are calculated from the signal initially filtered through the Wavelet filter. They represent, respectively, the energy of the pre-processed strain in the band $]150,800]$ Hz and $]200,300]$ Hz. Interestingly, selected features are related exclusively to the energetic of the strain, suggesting that energy is the most relevant information for the detection of GWs in a short-term feature approach. Furthermore, it is not surprising that all selected features are linked to a frequency range that contains the peak amplitude of BBH waveforms, in the mass range explored in the present paper.

Unlike previous classifiers based on CNNs, the ground model itself does not provide a final prediction on a fixed time length segment, because it is strictly trained on local properties of the signal. For this reason, the generation of an interferometer alert requires additional dedicated layers for the refined interpretation of the ground model information. Usually, the algorithm can be suitably adapted to the desired detection task, eventually by reaching a given trade-off between sensitivity and specificity of the classification. For example, lower levels should be extremely conservative in ruling out noise segments, in order to obtain the adequate level of sensitivity required for the successive offline analysis and parameter estimation.

In the proposed framework, regardless of the target task, either local or continuous analysis, a generic data segment from a GW interferometer is considered through n overlapped miniframes of fixed time length ($\omega_f = 100$ ms). For each miniframe, the reduced set of speech-processing features is extracted and used to compute a continuous output y (table 2), suitably constrained in the range $[0,1]$ as indicated in section 3.4. At this stage, the generic segment is matched to a stream of n consecutive outputs y_i . In this paper, we mostly focus on the local analysis, to have a more meaningful comparison with CNN approaches [30, 31, 33]. A typical stream of ground model output for a locally analyzed 1 s frame is shown in figure 3, for a signal+noise segment, and in figure 4 for a noise only segment. Both data segments used to produce the figures are extracted from the testing dataset, and are therefore good representative of the typical level 1 output. By inspecting the figures, it is clear that the observed trends for the model output, as a function of the miniframe, differ significantly for signal+noise and noise only segments. The first, figure 3, exhibits an increasing trend, with a group of consecutive miniframes having a much larger output. The latter has instead a more uniform output, reflecting the roughly uniform behavior of the detector noise in a short data segment.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1.1. L2T: a possible level-2 implementation based on a thresholded model

A first possible method for the interpretation of the ground model output, for the task of local analysis (performed by level 2), can be obtained by suitably thresholding its output. The principle is schematically explained in figures 3 and 4, where a threshold of τ = 0.68 is represented, as an example. After thresholding the output, one obtains a series of 19 categorical output values for each data segment to analyze (with values either 0 or 1). Following the different trend of model output for the two cases, one expects the appearance of clusters of 1 in the presence of a GW signal, while for noise only segments only some, spurious, isolated 1 is usually present, as seen in figures 3 and 4. One can exploit this difference in the observed patterns of the categorical output to derive a model for the thresholded interpretation of the ground model. For the thresholded interpretation model, we used maximum length of a cluster of 1, the maximum length of a cluster of 0 and the total number of 1 observed in the segment. Boundary conditions on these three variables can be derived through any ML algorithm, exploiting a dataset of previously labeled segments. We name this approach multi-modeling. To this end, we exploited 1 s data segments, previously simulated, and subdivided into training, cross-validation and testing datasets, dedicated to the derivation of the interpretation models.

Using the derived boundary conditions on the above mentioned quantities, the performance of the thresholded version of level 2 (L2T) can be evaluated using the segments of the dedicated testing dataset. Figure 5 shows the receiving operator characteristic (ROC) curves obtained using L1 and the thresholded version of L2 (L1+L2T) for optimal SNR ranging from $\rho_{opt} = 6$ to $\rho_{opt} = 14$, individually. In the graph, the Y-axis is the specificity of the classifier, i.e. the fraction of signal+noise segments that are correctly classified as containing a GW signal. This is a crucial quantity in real-time GW detection applications. X-axis is related to the capability of the model to rule out noisy segments. It corresponds to the fraction of noise-only samples incorrectly classified as containing a GW signal. For the latter, 0 indicates a complete, ideal, noise rejection, while 1 indicates that all noise-only segments in the dataset produce a zero-level trigger. Each point in the ROC curve is obtained with a certain value of the threshold τ. A threshold value τ = 1 would result in a point located at the bottom-right corner of the graph, i.e. all segments are identified as noise, while τ = 0 will identify all segments as containing a GW, top-right corner. An optimal LT2 can be obtained selecting the proper threshold, according to the required performance. For example, for an optimal SNR of $\rho_{\textrm{opt}} = 12$, one can easily rule out $90 \%$ of the noise-only sample, while identifying nearly-$100 \%$ of the GW signals. This value drops to about $70 \%$ for $\rho_{\textrm{opt}} = 8$. In general, L1+L2T exhibits a good regularity for increasing SNRs, with a significant rise in performance towards higher SNR values. Furthermore, the best trade-off between specificity and sensitivity is obtained at a similar threshold value at all SNR, consistently.

Zoom In Zoom Out Reset image size

4.1.2. L2U: an alternative level-2 implementation with an unthresholded model

An alternative way to interpret the ground model output can be achieved without thresholding the values of level 1. In fact, thresholding the output typically results in a loss of information, being the continuous output converted into a categorical one with only two possible categories. To derive the unthresholded model, we considered the ground model output values for 19 consecutive miniframes in which the data segment is subdivided. Said y_i the output of the ith miniframe, we constructed the curve identified by the set of points $\{(i,y_i) : i\in[1,19]\}$. To derive the present version of the L2U model, we decided to build a set of 260 features to represent the information on the presence of a GW signal in the noisy segments. Such features included simple statistical values (mean, median, standard deviation, kurtosis) and linear fit coefficients. The same set of features is also extracted on multiple smoothed version of the continuous output curve, obtained with different smoothing, as described in B.2. Exploiting this set of features, a second model was derived through BP, using the data segments contained in the same subset of the testing dataset previously used to derive the thresholded model. The deriving model can be used to build an alternative unthresholded level 2 (L2U). The analytical implementation of the unthresholded model is shown in table 3, together with a description of the selected features. A more detailed description of the features used by the unthresholded model is instead reported in the appendix. Similarly to the case of the ground model, also the unthresholded model exploits an extremely reduced number of features. Features are reduced from the initial 266, in the training set, to only 6. The performance of the deriving L2U in the task of identifying GW data segments is probed using the same dataset previously used to test the performance of L2T. This time, a given trade-off between sensitivity and false positive rate can be obtained by using a threshold η directly on the output of the unthresholded model. The corresponding ROC curves, in which each point corresponds to a given η value, are shown in figure 6(solid lines), compared to those of L1+L2T (dashed lines). Curves obtained at analogous SNR are drawn with the same color.

Zoom In Zoom Out Reset image size

Table 3. The analytical formulation of the unthresholded model derived via BP. A brief description of the selected features is also provided. The unthresholded model foresees only 6 features. See the appendix for a detailed description of the features selected by BP.

$Y = \tanh \Big\{ 1.83 \cdot \max \left[ \; c_1 , \; \max( c_2, b_4) \; \right] \Big\} ^{19.5073\cdot \arctan(b_6)}$
where $c_1 = 1.829 \cdot b_4 \cdot \max\Big\{1.829 \cdot b_4 \cdot \tan(b_4) , \;\; \xi + 0.947\,062 \cdot b_2 \Big\} + b_2$,
$\xi = 1.732\,176 \cdot b_4 \cdot \sin [ \max( \tan(b_3) , b_4 ) ]$
and $c_2 = \{1.829 \cdot \sinh [ \textrm{erf} ( b_3 \cdot b_1 ) ] \cdot \textrm{erf} [ \tan(b_4) ] \cdot b_5 \} + b_2$
The features:	b1 = $russo8_{3}$	b2 = $russo9_{5}$
	b3 = $sm2ampl$	b4 = $sm3weav$
	b5 = $sm5barone7_{2}$	b6 = $sm5barone9_{0}$

An inspection of the ROC curves for the two modules suggests that the performance of L1+L2U is significantly better than that obtained with L1+L2T (and therefore by thresholding the ground model output). As an example, for an optimal SNR $\rho_{opt} = 8$ and a false positive rate of $10 \%$, the sensitivity of L1+L2U raises to about $90 \%$, while L1+L2T has a sensitivity of only $70 \%$. Interestingly, for optimal SNR $\rho_{opt} = 12$, the sensitivity of L1+L2U is almost saturated to $100 \%$ at nearly all false alarm probability values, while, using the thresholded module leads to a sensitivity below $90 \%$ when the false alarm probability approaches $5 \%$. This indicates that the new set of features used for the L2U model has a much larger informative content than the simple quantities used in the L2T model. It is clear that other alternative interpretations of the ground model output (and therefore other possible level 2 implementations), other than the ones proposed here, could lead to even more accurate predictions.

The local analysis described in the previous section, i.e. the analysis of data segments of fixed length, allows a direct comparison of L1+L2 derived in the present work with other state-of-the-art approaches available in the literature. We compare the results achieved with L1+L2U to those of matched-filtering and state-of-the-art artificial intelligence approaches based on CNNs previously exploited in GW detection. We consider the CNN model described in [31]. This model is derived and tested exclusively on data segments of 1 s length, where the GW peak lies in the fractional range [0.75, 0.95] of the time series. For this reason, it is ideal to be used for the local analysis proposed in this section. Matched-filtering is a key algorithm, usually adopted for the detection of binary coalescence signals [51, 52], which is based on the coherence of the detected data segment with a given template waveform. The ranking statistics used for matched-filtering is the matched-filtering SNR numerically maximized over arrival time, phase and distance (see [31] for additional details). Given the ranking statistic of each analysis algorithm, it is possible to construct the ROC curves by evaluating, for each threshold, the number of true positives and true negatives achieved by the classifier.

In figure 7, we show the ROC curves obtained with the present framework using the most accurate module derived in this work, together with the corresponding ROC curves obtained using standard matched-filtering and the model described in [31]. We focused the comparisons on data containing signals with optimal SNR ranging from $\rho_{opt} = 6$ to $\rho_{opt} = 10$, while, for matched-filtering and the CNN model, we also show the results for $\rho_{opt} = 4$, which is well below the currently considered confident detection limit ¹² . The curves obtained at $\rho_{opt} = 6$ with matched-filtering (blue dashed line) and the CNN model (light blue dashed line) show similar values of sensitivity and specificity. The sensitivity of the CNN slightly exceeds that of matched-filtering at large false alarm probability values, while the trend is opposite at small values of false alarm probability, i.e. at large sensitivity values. The ROC curves for L1+L2U are shown with the red, yellow and green curves for $\rho_{opt} = 6$, 8, and 10, respectively. Compared to the CNN model at the same SNR, CNNs have better performance than our model for all values of the false alarm probability. For example, at a false alarm probability of 10⁻², and for $\rho_{opt} = 6$, L1+L2U is able to identify about $45 \%$ of the GW data segments in the testing dataset, while CNNs identify about $65 \%$ of the GW signals. At optimal SNR $\rho_{opt} = 4$ (dash-dotted lines), the performance of CNNs drops significantly compared to our results at $\rho_{opt} = 6$. If we focus on conventionally considered confident SNR values, i.e. $\rho_{opt}\unicode{x2A7E}8$, the performance of our simple classifier is higher than that of CNNs and matched-filtering at $\rho_{opt} = 6$. At $\rho_{opt} = 10$, we obtain an almost saturated sensitivity at all values of false alarm probability ${\unicode{x2A7E}}10^{-3}$. However, one must consider that the proposed models are outstandingly simple. Finally, it is important to stress that the situation of local analysis, which is a particular ideal benchmark for the CNN model of [31], is a particularly unfavorable condition for our framework, which is instead designed for a short-term analysis of features, thus not accounting for the duration of the strain and the position of the signal peak within the segment. The assumptions used to derive our models are ideally more advantageous for the analysis of a continuous timeseries in a realistic analysis pipeline. In addition, as shown previously in the paper, the derived models (for both level 1 and level 2) have the further advantage of an outstandingly-low computational complexity, which makes them ideal for the integration in real-time, low-latency, analysis pipelines.

Zoom In Zoom Out Reset image size

For each 1 s sample there are 19 continuous miniframe outputs y_i ($i = 1 {\ldots} 19$), usually in the range [0,1]. Let us call Γ the curve given by the set of points $(i, y_i)$:

$$\begin{align*} \Gamma = \left\{\left(i, y_i\right) : i = 1 {\ldots} 19 \right\} .\end{align*}$$

The original curve given by the miniframe outputs can be smoothed considering the mean value on the Y-axis of p points, p > 1. By definition, a mean-smoothed curve $\Gamma_p$ is the set of points $(j,Y_j)$, for $j = 1 {\ldots} (19+1-p)$, in which every point Y_j is given by the mean Y-axis value of p adjacent points of Γ:

$$\begin{align*} Y_j & = \frac{1}{p} \sum_{k = 0}^{p-1} y_{j+k} ,\end{align*}$$

$$\begin{align*} \Gamma_p & = \left\{\left(j, Y_j\right) : j = 1,{\ldots},19+1-p \right\} .\end{align*}$$

Because of this definition, the not-smoothed curve Γ is a smoothed curve Γ₁ with p = 1. For a given curve $\Gamma_x$, let us consider α its amplitude, $\alpha_x = \max\Gamma_x - \min\Gamma_x$, and its minimum value $\mu_x = \min\Gamma_x$

Each point of the curve is processed by threshold, so that every point above the threshold value τ is considered as a 'one', otherwise as a 'zero'. Within each curve, except for the case in which the threshold value is grater than the minimum value of the curve, there will be a 'first one' and a 'last one' in the sequence. Let ξ be the distance between the 'first one' and the 'last one' of the curve on the X-axis.

Considering Γ (not-smoothed curve):

$$\begin{align*} b_1 = \frac{1}{\xi}\sum_{ones} y_i \end{align*}$$

where the ones are assigned using a threshold value $\tau = \mu + 0.8\cdot\alpha$. 

$$\begin{align*} b_2 = \sum_{i = 0}^{19} \left( y_i - \left(0.9 \cdot \alpha + \mu\right)\right) \end{align*}$$

Considering a smoothed curve with p = 2:

$$\begin{align*} b_3 = \alpha_2 \end{align*}$$

Considering a smoothed curve with p = 3:

$$\begin{align*} b_4 = \sum_{j = 1}^{17} \left(Y_j \cdot w_j\right) \Bigg/ \sum_{j = 1}^{17} w_j \end{align*}$$

where $\{w_j\}$ is a suitable weighting vector to account for the enhancement of the signal power from the left to the right part of the data segment, as a result of the topology of the BBH waveform.

Considering a smoothed curve with p = 5 and $\tau = 0.7\cdot\alpha + \mu$, b₅ is the length of the shortest cluster of zeros.

Considering a smoothed curve with p = 5 and $\tau = 0.7\cdot\alpha + \mu$, b₆ is the number of the ones within the curve.