Sparse Representation of Gravitational Sound

Gravitational Sound clips produced by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and the Massachusetts Institute of Technology (MIT) are considered within the particular context of data reduction. It is shown that these types of signals can be approximated at high quality using much less elementary components than those required within the standard orthogonal basis framework. Furthermore, a measure a local sparsity is shown to render meaningful information about the variation of a signal along time, by generating a set of local sparsity values which is much smaller than the dimension of the signal. This point is stressed further by recourse to a more complex signal, generated by Milde Science Communication to divulge Gravitational Sound in the form a ring tone.


Introduction
In 1905 Henri Poincaré first suggested that accelerated masses in a relativistic field should produce gravitational waves [1]. The idea was magisterially pursued by Einstein via his celebrated theory of general relativity. In 1918 he published his famous quadrupole formula resulting from the calculation of the effect of gravitational waves [2]. A century later, the LIGO Scientific Collaboration and Virgo Collaboration published a paper about the gravitational radiation they had detected on September 2015 [3]. Ever since scientists believe to have entered in a new era of astronomy, whereby the universe will be studied by 'its sound' [4][5][6][7][8]. Gravitational Sound (GS) signals will then be here scrutinized with advanced techniques.
In the signal processing field, the problem of finding a sparse approximation for a signal consists in expressing the signal as a superposition of as few elementary components as possible, without significantly affecting the quality of the reconstruction. In signal processing applications the approximation is carried out on a signal partition, i.e., by dividing the signal into small pieces and constructing the approximation for each of those pieces of data. Traditional techniques would carry out the task using an orthogonal basis. However, enormous improvements in sparsity can be achieved using an adequate over-complete 'dictionary' and an appropriate mathematics method. For the most part, these methods are based on minimization of the l 1 -norm [9] or are greedy pursuit strategies [10][11][12][13][14][15][16][17], the latter being much more effective in practice.
Sparse signal representation of sound signals is a valuable tool for a number of auditory tasks [18,19]. Moreover, the emerging theory of compressive sensing [20][21][22] has enhanced the concept of sparsity by asserting that the number of measurements needed for accurate representation of a signal informational content decreases if the sparsity of the representation improves. Hence, when some GS tones made with the observed Gravitation Wave (GW) were released, we felt motivated to produce a sparse approximation of those clips.
We simply analyze the GS tones from a processing viewpoint, regardless on how and why they have been generated. We consider a) a short tone made with the chirp gw151226 that has been detected, b) the theoretical simulated theoretical GS, iota 20 10000 4 4 90 h, and c) the Black Hole Billiards ring tone, which is a more complex signal produced by superposition with an ad hoc independent percussive sound. The ensuing results are certainly interesting. If, in the future, GS signals are to be generated at large scale (as astronomical images have been produced [23, 24]), it is important to have tools for all kinds of processing of those signals.
The central goal of this Communication is to present evidences of the significant gain in sparsity achieved if a GS signal is approximated with high quality outside the orthogonal basis framework. For demonstration purposes we have made available the MATLAB routines for implementation of the method.

Some Preliminary Considerations
The traditional frequency decomposition of a signal given by N sample points, f Finding a sparse solution is the goal of sparse approximation techniques.
The problem of the sparse approximation of a signal, outside the orthogonal basis setting, consists in using elements of a redundant set, called a dictionary, for constructing an approximation involving a number of elementary components which is significantly smaller than the signal dimension. For signals whose structure varies with time, sparsity performs better when the approximation is carried out on a signal partition. In order to give precise definitions we introduce at this point the notational usual conventions: R and C represent the sets of real and complex and numbers, respectively.
Boldface fonts are used to indicate Euclidean vectors and standard mathematical fonts to indicate components, e.g., d ∈ C N is a vector of N-components d(i) ∈ C N , i = 1, . . . , N. The operation ·, · indicates the Euclidean inner product and · the induced norm, i.e. d 2 = d, d , with the usual inner product definition: A partition of a signal f ∈ R N is represented as a set of disjoint pieces, f q ∈ R N b , q = 1, . . . , Q, henceforth to be called 'blocks', which, without loss of generality, are assumed to be all of the same size and such that QN b = N. Denoting byĴ the concatenation operator, the signal f ∈ R N is 'assembled' from the blocks as f =Ĵ Q q=1 f q . This operation implies that the first N 1 components of the vector f are given by the vector f 1 , the next N 2 components by the vector f 2 , and so on.
A dictionary for R N b is an over-complete set of (normalized to unity) elements D = {d n ∈ R N b ; d n = 1} M n=1 , which are called atoms.

Sparse Signal Approximation
Given a signal partition f q ∈ R N b , q = 1, . . . , Q and a dictionary D, the k q -term approximation for each block is given by an atomic decomposition of the form The approximation to the whole signal is then obtained simply by joining the approximation for the

The Method
The problem of finding the minimum number of K terms such that f −f K < ρ, for a given tolerance parameter ρ, is an NP-hard problem [12]. In practical applications, one looks for tractable sparse solutions. For this purpose we consider the Optimized Hierarchical Block Wise (HBW) version [25] of the Optimized Orthogonal Matching Pursuit (OOMP) [13] approach. This entails that, in addition to selecting the dictionary atoms for the approximation of each block, the blocks are ranked for their sequential stepwise approximation. As a consequence, the approach is optimized in the sense of minimizing, at each iteration step, the norm of the total residual error f − f K [25]. As will be illustrated in Sec. 3.3, when approximating a signal with pronounced amplitude variations the sparsity result achieved by this strategy is remarkable superior to that arising when the approximation of each block is completed at once, i.e., when the ranking of blocks is omitted. The OHBW-OOMP method is implemented using the steps indicated below.
OHBW-OOMP Algorithm 1) For q = 1, . . . , Q initialize the algorithm by setting: r 0 q = f q , f 0 q = 0, Γ q = ∅ k q = 0, and selecting the 'potential' first atom for the atomic decomposition of every block q as the one corresponding to the indexes ℓ q 1 such that 2) Use the OHBW criterion for selecting the block to upgrade the atomic decomposition by adding one atom If 3) Calculate Upgrade the set Γ q ⋆ ← Γ q ⋆ ∪ ℓ k q ⋆ +1 and increase k q ⋆ ← k q ⋆ + 1.
4) Select a new potential atom for the atomic decomposition of block q ⋆ , using the OOMP criterion, i.e., choose ℓ q k q ⋆ +1 such that including, for numerical accuracy, the re-orthogonalizing step: 6) Check if, for a given K and ρ either the condition Q q=1 k q = K +1 or f −f K < ρ has been met. If that is the case, for q = 1, . . . , Q compute the coefficients c kq (n) = b kq n , f q , n = 1, . . . , k q . Otherwise repeat steps 2) -5).
Remark 1: For all the values of q, the OOMP criterion (6) in the algorithm above ensures that, fixing the set of previously selected atoms, the atom corresponding to the indexes given by (6) minimizes the local residual norm f q − f kq q [13]. Moreover, the OHBW-OOMP criterion (3), for choosing the block to upgrade the approximation, ensures the minimization of the total residual norm [25]. Let us recall that the OOMP approach optimizes the Orthogonal Matching Pursuit (OMP) one [11]. The latter is also an optimization of the plain Matching Pursuit (MP) method [10](see the discussion in [13]).

The Dictionary
The degree of success in achieving high sparsity using a dictionary approach depends on both, the suitability of the mathematical method for finding a tractable sparse solution and the dictionary itself.
As in the case of melodic music [25,26], we found the trigonometric dictionary D T , which is the union of the dictionaries D C and D S given below, to be an appropriate dictionary for approximating these GS signals. In the above sets w c (n) and w s (n), n = 1, . . . , M are normalization factors.
We also found that sparsity may benefit by the inclusion of a dictionary which is constructed by translation of the prototype atoms, p 1 , p 2 and p 3 in Fig. 1. Denoting by D P 1 , D P 2 and D P 3 the dictionaries arising by translations of the atoms p 1 , p 2 , and p 3 , respectively, the dictionary D P is built as D P = D P 1 ∪ D P 2 ∪ D P 3 . The whole mixed dictionary is then D M = D T ∪ D P , with D T = D C ∪ D S . Interestingly enough, the dictionary D M happens to be a sub-dictionary of a larger dictionary proposed in [27] for producing sparse representations of astronomical images, the difference being that, in this case, sparsity does not improve in a significant way by further enlarging the dictionary. From a computational viewpoint the particularity of the sub-dictionaries D C and D S is that the inner product with all its elements can be evaluated via FFT. This possibility reduces the complexity of the numerical calculations when the partition unit N b is large [25,26]. Also, the inner products with the atoms of the dictionaries D P 2 and D P 3 can be effectively implemented, all at once, via a convolution operation.
Note: The MATLAB routine implementing the OHBW-OOMP approach, dedicated to the dictionary introduced in this section, has been made available on [28].

The Processing
We process now the three signals we are considering here:  The quality of an approximation is measured by the Signal to Noise Ratio (SNR) which is defined as SNR = 10 log 10 f 2 f − f K 2 = 10 log 10 The sparsity of the whole representation is measured by the Sparsity Ratio (SR) defined as SR = N K , where K is the total number of coefficients in the signal representation as defined above.

Audio representation of the chirp gw151226
This clip, made with the detected short chirp gw151226, is plotted in the left graph of Fig.2. The graph on the right is its classic spectrogram. When an orthogonal basis for approximating these The graph on the right is its classic spectrogram. When processing the signal with DCT the best sparsity result when the approximation is carried out block by block up to the same error is SR=4.2, for SNR=40dB, and corresponds to N b = 16384. However, with N b = 2048 the OHBW version for  Fig. 4 represents the difference between the signal and its approximation, up to SNR=40dB. It is worth commenting that, if with the same dictionary, the approximation were carried out without ranking the blocks, i.e., approximating each block at once up to the same SNR, the value of SR would be only 6.7. This example highlights the importance of adopting the OHBW strategy for constructing the signal approximation, when the signal amplitude varies significantly along the domain of definition.

The Role of Local Sparsity
The SR is a global measure of sparsity indicating the number of elementary components contained in the whole signal. An interesting description of a the signal variation is rendered by a local measure of sparsity. For this we consider the local sparsity ratio sr(q) = N b kq , q = 1, . . . , Q where, as defined above, k q is the number of coefficients in the decomposition of the q-block and N b the size of the block.  Black Hole Billiards ring tone. The lines in the right graph discriminate the inverse local sparsity ratio produced with atoms in the dictionary D P (blue), in the dictionary D T (red) and in the whole dictionary D M (black).
Since the Black Hole Billiards ring tone is a more complex signal, due to the superposition of the artificial sound, the information given by the local sparsity ratio is richer than in the previous cases. Notice for instance that the dark line in the left graph of Fig. 6 clearly indicates the offsets in the percussive part of the clip which has been superimposed to the GS chirp. Moreover this line, joining 512 points of inverse local sparsity ratio, also roughly follows the signal variation envelop.
The graph on the right discriminates the local sparsity measure corresponding to atoms in the trigonometric component of the dictionary, and those in the dictionary D P . From bottom to top the first line (blue) represents the inverse local sparsity values corresponding to atoms in D P and the next line (red) to atoms in D T . The top line (black) corresponds to atoms in the mixed dictionary D M for facilitating the visual comparison. In this clip 20% of atoms are from dictionary D P and, as indicated by the blue line in the right graph of Fig. 6, a significant contribution of those atoms takes place within the blocks where the rapid rise of the GS tone takes place (c.f. spectrogram in Fig. 4).

Conclusions
We have here advanced an effective technique for the numerical representation of Gravitational Sound clips produced by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and the Massachusetts Institute of Technology (MIT). Our technique is inscribed within the particular context of sparse representation and data reduction. We laid out a detailed procedure to this effect and were able to show that these types of signals can be approximated with high quality using significantly fewer elementary components than those required within the standard orthogonal basis framework.