Nonorthogonal wavelet transformation for reconstructing gravitational wave signals

Detections of gravitational-wave signals from compact binary coalescences have enabled us to study extreme astrophysical phenomena and explore fundamental physics. A crucial requisite for these studies is to have accurate signal models with characteristic morphologies, which have been challenging for many decades, and researchers are still endeavoring to incorporate important physics. Therefore, morphology-independent methods have been developed for identifying a signal and its reconstruction. The reconstructed signal allows us to test the agreement between the observed signal and the waveform posterior samples from parameter estimation. These methods model observed signals using a nearly orthogonal wavelet basis in the frame of continuous wavelet transformation. Here, we propose log-uniform scales to construct the wavelets, which are are highly redundant (nonorthogonal) compared to the conventional octave scales but more efficient for reconstructing the signals at high frequencies. And we introduce a semi-model-dependent reconstruction method using the posterior samples of the events, where we model the signal using Gabor-Morlet wavelets with log-uniform scales. We demonstrate the ability to detect deviation using a numerical simulation of an eccentric binary black hole merger, where the signal in the data does not belong to the search template waveform manifold. Finally, we apply this method to each binary black hole merger event in GWTC-1. We have found that the signal produced by the GW150914 event has 96% agreement with the waveform posterior samples. As the detector sensitivity improves and the detected population of black hole mergers grows, we expect the proposed method will provide even stronger tests.


I. INTRODUCTION
Detecting gravitational waves (GWs) from the merging compact binaries consists of black holes and neutron stars has ushered in a new era of observational astronomy [1,2].So far, the network of Advanced LIGO and Advanced Virgo Detectors has observed more than 90 of these events [3][4][5][6][7][8].Observations of compact binary mergers allow us to design a range of novel tests of general relativity (GR) in the strong-field regime [9][10][11].A current limitation on tests of beyond-GR physics with compact binary coalescences is the lack of understanding of the strong-field merger regime in nearly all modified theories of gravity.Thereby, we generally perform the tests assuming GR to be the null hypothesis or evaluate the consistency of the data with predictions from the theory.At the same time, the latter test can reveal if any unknown binary parameter influences the observed signal.
The continuous wavelet transformation (CWT) based time-frequency analysis has played a crucial role in identifying the gravitational wave events and detecting the deviation from general relativity.The most potential gravitational wave sources, such as compact binary mergers, supernovae explosions, produce non-stationary signals.The CWT is a powerful analysis tool that allows us to obtain a time-frequency localized projection * soumen.roy@nikhef.nl;roy.soumen61@gmail.com of a non-stationary signal.This procedure calculates the wavelet coefficients by performing a scalar product between the signal and wavelet basis.These coefficients are used not only for time-frequency analysis but also to reconstruct the original time-domain signal.The existing time-frequency-based methods such as coherent wave burst (cWB) [12][13][14], BayesWave [15,16], and X-pipeline [17,18] are designed for searching unmodeled gravitational-wave bursts and reconstructing the signal from the coherent response of a network of detectors.The previous studies reported satisfactory agreement between the reconstructed signal and the estimated waveform of binary black hole mergers for the louder events [3,19,20].
The conventional burst search method identifies a cluster of pixels from the time-frequency map of a detector's strain data if each of those pixels has more power than one expects from the noise alone, called triggers.If the triggers from a network of detectors are coincident, then the method claims a detection of a signal [21].The cWB method first constructs a multi-resolution time-frequency map of the strain data using Wilson-Daubechies-Meyer wavelets [22] and uses the same strategy to detect the events.After that, the method employs a coherent framework of the constrained maximum likelihood approach to yield the event properties: sky location, wave polarization, and signal reconstruction [13,23].Whereas the BayesWave method employs a framework of Bayesian statistics without relying on any prior assumption of waveform morphology [15].This method reconstruct the gravitational wave signals and instrumental noise, where both of them are characterized as a superposition of Gabor-Morlet wavelets or chrirplets [24].The number of time-localized wavelets and their parameters are determined via Reversible Jump Markov Chain Monte Carlo algorithm.After that, whether the likelihood of the event is favoured to a true gravitational wave signal or instrumental glitch is determined by using Bayesian model selection strategy.Finally, the posterior samples are used to reconstruct the signal as a superposition of the wavelets.
A wavelet can be regarded as a time-frequency localized bandpass filter.The bandwidth of which is determined by the properties of that wavelet.The reconstructed signal has an inevitable noise contribution that passes through those wavelet filters.The match between the reconstructed signal and the estimated waveform increase with the signal-to-noise ratio (SNR).However, we would never achieve absolute agreement (unit match) even if our predicted theory is complete because of the additional noise that passes through the coherent wavelets.
In this paper, we propose a quasi model-dependent method for reconstruction of signal from noisy data.We identify the essential wavelets using the binary black holes merger waveforms estimated using the standard Bayesian analysis [25][26][27][28].The unmodelled reconstruction methods identify the wavelets that are coherent across the network of detectors.The collection of those coherent wavelets in the time-frequency plane form a cluster.Similarly, the essential wavelets form a cluster in the time-frequency plane which can be collectively looked as a single patch with an irregular boundary.
The patch area determines how much noise persists in the reconstructed signal.The wavelets are discretely placed over the time-scale plane, where the patch area is determined by the placement method and sparseness.The dyadic grid placement is considered to be the most efficient method, leading to the construction of an orthonormal wavelet basis.This placement assumes an output of octave decomposition (power-of-two logarithmic) to construct the grid.The sparseness of the grid imposes a limit on the acceptable loss in signal characteristics.We target to achieve maximum sparseness and still have an adequate representation.A highly sparse grid might provide a sufficient signal representation, but the area covered by the essential wavelets would be larger than the case of a denser grid.On the contrary, an overdense grid can not reduce the patch area beyond a specific limit and also increases the computational cost.
In this paper, we propose a log-uniform scale to place the wavelets over the time-scale domain.Further, we derive an inverse wavelet transformation formula for loguniform scale.The wavelets with such scales are highly redundant (i.e., over complete), but they can further reduce the area covered by the essential wavelets.As a result, We can reduce the noise contribution in the recon-structed signal, which leads to a more accurate signal reconstruction.Moreover, the essential wavelets with loguniform scale can provide a more appropriate signal representation than the octave scale for the high-frequency linear chirp signals.It implies that the log-uniform scale is more efficient for reconstructing the signals from low mass binary systems.
This paper is organized as follows: Section II describes the reconstruction methods using the continuous wavelet transformation with octave and log-uniform scales; Section III demonstrates the performance of the reconstruction methods for the linear chirp signals and gravitational wave signals in simulated Gaussian noise; Section IV describes the efficiency of the wavelet reconstruction method when a signal in the data does not belong to the manifold of search template waveform by injecting an eccentric waveform; Section V exhibits results of our analysis for binary black hole merger events observed by the Advanced LIGO and Virgo during the first and second observation runs.Finally, we conclude in Section VI.

II. RECONSTRUCTION FORMULAE FOR NON-ORTHOGONAL WAVELET
The CWT of a continuous signal x(t) is a linear mapping onto a time-scale space of wavelets: where, ψ * is the complex conjugate of the shifted and scaled version of the time-localized mother wavelet ψ, in which a and b are the scale and time-shift parameters, respectively.Therefore, the CWT provides a time-scale representation of the signal by performing a sliding crosscorrelation with a continuous family of wavelets.The quantity 1/ √ a outside the integral is an energy normalized factor.It assures that each wavelet has the same energy, whatever the value of scaling and shift.The function ψ(t) must satisfy a set of mathematical criteria to be a wavelet, it must have finite energy, and must hold the admissibility condition, where, ψ(ω) is the Fourier transform of ψ(t), ω = 2πf is the angular frequency, and C ψ is called the admissibility constant.The above condition implies that ψ(ω) is endeavoured to zero faster than ω and must not have zero frequency component, ψ(0) = 0.If the wavelet function satisfy both the criteria, then the inverse wavelet transform can be described as a superposition of the dual wavelets: An alternative approach of inversion formula was found by Morlet, we can even choose a a completely different wavelet function, Dirac δ-function δ ((b − t)/a) instead of the analysing wavelet, and that leads to a single integral inverse formula [29], where C δ is the admissibility constant of the δ-function that can be computed using the Eq.( 3).In this work, we use this single integral inversion formula for reconstructing the gravitational wave signals.
The wavelet function can be regarded as a impulse response of a bandpass filter.The associated frequency of the wavelet can be treated as frequency value in the timefrequency domain, which is known as pseudo-frequency (f p ).It depends on the central frequency (f c ) and the scale parameter of the wavelet, This equation can used to represent the CWT in timefrequency frame as like the short-time Fourier transform.One of the most commonly used mother wavelet is the Morlet wavelet [30], which is a complex non-orthogonal wavelet.As the complex wavelet can separate the phase and amplitude components associated with the signal, it is more suited to determine the instantaneous frequency.The Morlet wavelet consists of an harmonic oscillation with Gaussian window, where ) −1/2 the normalization constant.The second term in the bracket is a correction to preserve the admissibility condition.It ensures that the zero-frequency component vanishes.In practice, we ignore this term as it is approximately zero for the values of f c 0. Note that the above equation does not contain the time-shift and scale parameters.When we analyze a signal, we replace the variable t with (t − b)/a.
The central frequency f c is a crucial parameter in the time-frequency analysis to regulate the tradeoff between temporal precession and spectral resolution.It controls the number-of-cycles of the wavelet without modifying the shape of Gaussian window.A large value of central frequency leads to an increased temporal precision at the cost of decreased spectral precision and vice versa for a smaller value of central frequency.It is impossible to achieve simultaneously good precision in time as well as in frequency.A rule of thumb can be proposed for analyzing the gravitational waves from compact binary mergers.The signals from low-mass systems spend several tens of seconds or more than a hundred seconds in the bandwidth of Advanced LIGO-like detectors, which look like a long-duration chirp signal.A larger value of central frequency is convenient to achieve an overall good timefrequency precision for those systems.On the other hand, a smaller value of central frequency is suitable for high mass systems as their signals are short-duration bursts with a tiny chirp.

A. Choice of wavelet scales: octave
The CWT generally suffers due to the redundancy at large scales, where the neighboring wavelets are highly correlated.Once a wavelet function is chosen, it is important to choose a tightest set of scales that forms an orthonormal wavelet set.The key mathematical criteria to choose a set of discrete wavelets is that every function f ∈ L 2 (R) 1 must be fully expressed as a superposition of those wavelets i.e., the set is complete in L 2 (R).The conventional approach of choosing the scales is an output of octave decomposition [31,32], where N is the total number of samples in x(t), δt is sampling interval, a 0 = 2δt f c is the smallest resolvable scale, and δj is the spacing between the discrete scales.
For Morlet wavelet, an adequate scale resolution can be achieved for the values of δj 0.5.The quantity J that determines the highest value of the scale is associated with the Nyquist frequency.The octave scale leads to a simpler implementation of inverse discrete wavelet transformation as the quantity da/a = a 0 δj ln( 2) is a constant.Now, the inverse formula Eq. ( 5) yields,

B. Choice of log-uniform scale
For non-orthogonal wavelet analysis, one can adopt an arbitrary set of scales to build up a more complete picture if that set satisfy the completeness criteria.We propose a log-uniform scale, The quantity δ = δf p /f c stands for the inverse of spacing between two discrete wavelets and smallest resolvable scale a 0 = f c /δf p .In this work, we choose δf p = 1/T , where T is duration of the time-series.The log-uniform scale leads to an uniformly spaced pseudo-frequency.In this formalism, Morlet wavelet transformation is equivalent to the scaled Gabor transformation as the pseudofrequencies are uniformly spaced.For log-uniform scale da/a 2 = δ , thus, the discrete inverse formula Eq. ( 5) yields, C. Identifying the essential wavelets The conventional method for identifying essential wavelets is a straightforward approach, where one discards the low-magnitude wavelet coefficients by applying a threshold, which is known as wavelet thresholding.[33,34].This approach aims to remove the noise from data without affecting the basic features of the signal.In this work, we propose a slightly different approach.We employ the modelled waveforms that are obtained from the standard Bayesian parameter estimation analysis.We compute wavelet coefficients for a waveform and set a threshold on coefficient value to select the essential wavelets.The threshold is chosen such that the resultant power from the essential wavelets is equal to a percentage of total power of the spectrogram.We calculate the threshold (E * ) for a given value of fractional power loss (R * ) by solving an equation, We shall call that R * is the spectral loss parameter.As the number of essential wavelets is considerably fewer than the total number of wavelets, a significant amount of noise can be removed while preserving the basic features of the signal.

D. Occupied area by the essential wavelets
Reconstructing a signal using the essential wavelets trades with the sensitivity, by which we are able to remove the noise from data.However, the reconstructed 1.An illustration of the area covered by the cluster of essential wavelets is produced using a gravitational waveform of a nonspinning equal mass binary system with a chirp mass of 30M .The areas bounded by the log-uniform wavelets and octave scale wavelets are 13 and 15, respectively.For this computation, we have assumed aLIGO [35] noise curve with a fixed lower cutoff frequency of 20Hz, the waveform is generated using SEOBNRv4 model [36].We have considered R * = 0.05 for determining the essential wavelets.
signal would always contain a nominal amount of inseparable noise, which passes through the wavelet filters.As a wavelet is regarded as a time-localized bandpass filter, it can be seen as a patch on time-frequency plane.The area (time-frequency bandwidth) of that patch can be used to determine the amount noise released through that wavelet.Therefore, the total amount of inseparable noise can be estimated by calculating the total area covered by the essential wavelets.
The area on the time-frequency plane covered by a cluster of wavelets {(b i , a i )} n0 i=1 is related to their placement.For octave scale wavelet with a central frequency of f c , the covered area is: where, f pi is the i th pseudo frequency, f pi = f c /a i .Whereas, the time-frequency area governed by a cluster of wavelets {(b i , a i )} n0 i=1 with log-uniform scale is:  12).We used those values to choose the essential wavelets.
Fig. 1 shows that log-uniform scale occupy a smaller area than the octave scale.The figure indicates that the log-uniform scale allows ∼ 13% less noise in the reconstructed signal than the octave scale for a high mass system, can provide a nearly identical signal representation.

III. PERFORMANCE OF THE RECONSTRUCTION PROCEDURES
In this section, we demonstrate the performance of the wavelet based signal reconstruction and compared between the choice of Octave and log-uniform scale.First, we exhibit that the essential wavelets are adequate to represent a chirp signal containing a broad range of frequency.Second, we carry out injection analysis to evaluate the performance for real data analysis.The gravitational wave signals are drawn from the binary black merger and added to a simulated Gaussian noise.

A. Reconstruction of chirp signal without noise
We consider a simple chirp signal with a constant amplitude and phase is upto a quadratic order in time (t), where f 0 is the starting frequency of the signal and the time range is chosen to be between 0 and 2 s.Therefore, the ending frequency of the chirp is 3f 0 .This type of chirp signal is adequate to exhibit the performance of the reconstruction method over a broad range of frequencies.We consider five different cases to identify the essential wavelets.For each case, the total power contained in the essential wavelets equals a fraction of signal power.
It is inevitable to have an overall amplitude loss since the number of essential wavelets is a small subset of the set of wavelets representing the whole time-scale space.However, for a given value of R * , we can set the overall amplitude by looking at the amplitude loss when determining the essential wavelets.An alternative approach is to define normalized minimum squared error (NMSE) for a reconstructed signal x rec (t), where, the symbol • denotes the norm.We use this above equation to quantify the adequateness of the essential wavelets to characterize a linear chirp signal.In Fig. 2, we illustrate the accuracy of the wavelet reconstruction of the linear chirp signals.The solid and dashed lines correspond to the log-uniform scale and octave scale, respectively.We consider five different cases of fractional loss in spectrogram power to select the essential wavelets used for reconstruction.The reconstruction of high-frequency chirps using log-uniform scale wavelet is more accurate than the octave scale and vice versa for the low-frequency chirps.

B. Reconstruction of gravitational wave signal in simulated noise
We estimate the performance of the reconstruction methods for gravitational wave signals from compact binary mergers of equal mass nonspinning black holes.For each case, an identical signal is injected in many noise realizations of stationary Gaussian distribution weighted by Advanced LIGO zero-detuned high-power (aLIGO) design sensitivity [35].The waveforms are generated for an wide range of chirp mass between 10M and 40M using SEOBNRv4 model [36] with a fixed lower cutoff frequency of 20Hz.To determine the reconstruction accuracy with signal-to-noise ratio (SNR) of the injection (injected SNR), we choose a set of values between 5 and 50.The injected SNR (ρ inj ) of a waveform (h inj ) is defined as: where, hinj (f ) denotes the Fourier transform of h inj (t) and S n (f ) denotes the one-sided detector noise power spectral density.We reconstruct the signal from noisy data using the log-uniform scale wavelet transform, where the essential wavelets are selected using the Eq.(12).To quantify the signal reconstruction accuracy, we calculate the match between the injected signal h inj and reconstructed signal h rec , which is defined as a inner product between two normalized waveforms ( ĥ * = h * / h * | h * ): where the term with angular brackets represents the following inner product, where hinj (f ) and hrec (f ) denote the Fourier transform of h inj (t) and h rec (t), respectively.The match between two waveforms varies between -1 and 1, depending on their correlation but not their overall amplitudes.A match value of 1 indicates a perfect positive correlation such that two waveforms change with equal proportion, and 0 means no correlation.The equal proportion changes with reverse direction indicate perfect negative correlation, for which match value is -1.To quantify the agreement between the wavelet reconstruction for a network of detectors, we compute the network match (M net ) as given in [37], where k represents the k th detector.In this paper, we consider three detectors configuration of Hanford (H1), Livingston (L1), and Virgo (V1).We consider aLIGO design sensitivity for H1 and L1 detectors [35,38,39], and advanced Virgo design sensitivity for V1 [40,41].
Please note that we demonstrate the injection analysis in this section assuming the single detector with the aLIGO design sensitivity.
A whiten Gaussian noise is considered to be distributed normally over time-scale domain.This implies that the noise energy is equally distributed over time-scale domain.The amount of noise retains in the reconstructed signal is determined by the spectral loss parameter R * .A higher value of R * can further reduce the inseparable noise in the reconstructed signal.At the same time, the essential wavelets would not be able to represent the complete signal characteristics.It implies that the spectral loss parameter plays a role that imposes a limitation on achieving the maximum overlap between the injected and reconstructed signal.Therefore, we want to optimize this parameter.For a given source parameter and injected SNR, we perform the injection analysis with a set of R * values.The noise averaged match for a given injected signal h inj ( λ, ρ inj ) and R * can be defined as where, h rec [n i ] represents the reconstructed signal from the data of i th noise realization n i .To determine the essential wavelets for a given value of R * , we compute the wavelet coefficients of whitened waveform and plug in to Eq. (12).Fig. 3 shows the noise averaged match as a function of R * for a set of injected SNR values and the number of essential wavelets as a function of R * .If we fix the spectral loss parameter and increase the injected SNR, the signal contribution to each essential wavelet increases.At the same time, the average noise contribution remains the same since the number of essential wavelets and their properties does not change.This leads to an improvement in the reconstruction accuracy.A fixed SNR curve in the Fig. 3 indicates that a smaller value of R * allows many nonessential wavelets in the analysis.The signal contribution to those wavelets is trivial, but the inseparable noise in the reconstructed signal increases.On the other hand, a higher value of R * discards many moderate essential wavelets.Thus, the choice of R * is a tradeoff between the signal and noise contributions to the wavelets.To find an optimum value R, we maximize the quantity M (R * ) noise over R * .
FIG. 4. The figures shows the mismatch between the injected gravitational wave signal and reconstructed signal, where the filter wavelets were constructed using the log-uniform scale as defined in Eq. (10).For producing the injections, we have assumed the nonspinning equal mass binary black holes with a wide range of chirp mass as labelled in x-axis.The legend of the figures stand for the optimal SNR of the injected signal.We have injected an identical signal into many noise realization to produce the distribution of mismatch for each chirp mass.We have assumed the aLIGO noise curve with a fixed lower cutoff frequency of 20Hz.The waveforms are generated using the SEOBNRv4 model.Fig. 4 demonstrates the mismatch (1 − M ) in reconstruction for nonspinning binary systems with equal component masses.The solid curve represents the median of the mismatch distribution, and the shaded region shows the ±σ width of that distribution.For the case of a fixed optimal SNR, we can see that mismatch substantially decreases with an increase of chirp mass.It is intuitively expected: the number of essential wavelets for high chirp mass systems is fewer than low chirp mass systems.The waveform of a binary system with high chirp mass can be characterized using a few wavelets as the waveform is short.In contrast, the waveform of a low chirp mass system is longer, for which one requires a large number of wavelets to represent the signal.
As the injection chirp mass increases, the number of essential wavelets decreases, and they cover a smaller area over the time-frequency plane.The signal contribution to the essential wavelet coefficients increases with the injection chirp mass for a fixed SNR case.At the same time, the noise contribution decreases due to the smaller area.That explains why the reconstruction signal for a high chirp mass system is more accurate than a low chirp mass.

IV. IDENTIFYING THE DEVIATION
As of yet, we discussed the performance of wavelet reconstruction where injected waveforms were used for determining the essential wavelets.That results tell us the efficiency of the semi-model dependent wavelet reconstruction method when a signal in the data belongs to the search template waveform manifold.However, a signal in real data may not belong to that manifold.The deviations could arise due to: the influence of unknown binary parameters (such as eccentricity if the waveform model only considers the circular binary), missing physics in the waveform model, deviation from GR, or noise artifacts.We demonstrate one of such cases, an injection waveform simulation for an eccentric BBH merger.Since recent template-based analysis by LIGO-Virgo collaborations used quasi-circular waveform [4], we consider a recently developed IMRPhenomXP model to generate the template waveform [42].This model includes the effect of orbital precession but does not consider orbital eccentricity.
The inclusion of eccentricity in the injected signal leads to a deviation from the search template waveform manifold.To demonstrate this, we consider a numericalrelativity simulation of eccentric IMR waveform picked from Simulating eXtreme Spacetimes (SXS) catalog [43].This system is a nonspinning binary with mass ratio of 1.22.The numerical simulation is performed for an eccentric BBH system, where eccentricity evolves with time.The reference eccentricity is e ref = 0.194 as measured at a reference orbital frequency of M f 0 = 0.0137.In our analysis, we scale the simulation for a total mass of 60 M , which is suitable for ground based detectors.
We keep the system in the face-on configuration where the inclination angle is 0 • .As the component masses are nearly equal and zero inclination, the contribution of higher-order modes to the injected signal is negligible.This configuration assures us that the deviation enters only due to the orbital eccentricity.We inject the signal assuming a three detector configuration of Advanced LIGO and Advanced Virgo [35,[38][39][40][41].We set the source's sky position and distance such that the injected SNR in H1, L1 and V1 are 25, 20, and 15, respectively.Since the average result of our wavelet reconstruction over many Gaussian noise realizations with an identical injected signal leads to the case of zero noise, we consider zero noise realization to construct the data stream.We analyze the data using the standard Bayesian parameter estimation library LALInferenceNest [27], a Bayesian inference nested sampling code implemented in the LIGO Algorithm Library (LALSuite) [44].We use the python-based package PESummary [45] to process the data from parameter estimation analysis and generate the template waveform in the detector frame.
In order to reconstruct the signal from data stream, we follow these steps: 1.For each posterior sample, we generate the CBC template waveform in the detector frame.Since the GR allows only two polarization states, referred to as the plus (h + ) and cross (h × ), the time-domain response h I (t) of a given detector I is determined by the antenna response functions (F + I and F × I ) of those polarizations [46], where, α and δ refers to the source sky location in terms of right ascension and declination, ψ is the polarization angle, D L is the luminosity distance to the source, ι is the inclination angle of the binary plane, λ represents the set of intrinsic parameters of the binary system, and ∆t I (≡ ∆t I (α, δ, t)) is the travel time of the signal from geocenter to the detector.
2. Whiten the waveform weighted by the noise amplitude spectral density such that the norm of the whitened waveform is equal to its optimal SNR.
3. Compute the wavelet coefficients for each whiten template waveform using the CWT as shown in Eq. ( 1),, where wavelets are constructed using the log-uniform scale.
4. Determine the essential wavelets using Eq. ( 12), where spectral loss parameter R * is obtained from Eq. ( 21).We inject the best fit template (corresponds to maximum likelihood sample) waveform in many simulated noise realizations and estimate the value of R * for each detector's data.We have found the value of R * for H1, L1, and V1 is 2.5%, 4%, and 5%, respectively.Please note that we estimate the spectral loss parameter for the maximum likelihood sample and use that value for all the posterior samples.
5. Use those essential wavelets to reconstruct the signal from detector strain.
The above-described procedure is used to obtain a reconstructed signal for each posterior sample.These wavelet-based reconstructed signals are very similar, and their 90% interval is very thin and looks like a line.Therefore, we illustrate the median of reconstructed signals at every time index.2Fig. 5 illustrates the results from wavelet reconstruction, 90% credible region of LALInference template waveform, and injected waveform.We can see the amplitude and phase of the reconstructed waveform are approximately consistent with the injected waveform.On the contrary, the LALInference waveform is out of phase over a bit of the region, which indicates a significant deviation from the search template waveform manifold.We also compute three different matches using the reconstructed signal (h rec ), LALInference template waveform (h bif ) obtained from posterior samples as described in step 1, and injected waveform (h inj ): M (h inj , h bif ), M (h bif , h rec ), and M (h inj , h rec ).Table I summarizes the match comparison.It signifies the wavelet-based reconstructed waveform is more faithful than the LALInference.However, we can not compute the M (h rec , h inj ) and M (h bif , h inj ) for an actual event case as it is infeasible to know the true signal in data.We propose to compute the residual SNR (ρ res ) obtained by subtracting each template waveform of LALInference posterior samples from the corresponding wavelet-based reconstructed waveform.LVC commonly uses this procedure in the test of GR with BBH events [9][10][11]47].We report the median of residual SNR in the last column of Table I.

V. ANALYSIS OF EVENTS IN GWTC-1
We apply the proposed method to each binary black hole event in the first gravitational-wave transient catalog GWTC-1 [3] to reconstruct the signals from individual detectors.In this analysis, we use the on-source data from the Gravitational Wave Open Science Center [48,49], released for GWTC-1 [3].For determining the essential wavelets, we use the posterior samples of source properties obtained from Bilby's [50,51] reanalysis of GWTC-1.Parameter estimation analysis was performed using the IMRPhenomPv2 [52,53] waveform model, and PSD was estimated using the BayesLine algorithm [54].
In order to reconstruct the signal from data using the posterior samples, we follow the steps described in Sec- tion IV.Fig. 6 shows the results of GW150914 [1]: CBC template waveform from posterior sample and wavelet reconstruction.The agreement for H1 data is better than L1 data as the reconstruction accuracy increases with SNR.For H1 data, the most probable value (mode) of the match values is 0.975, and the maximum is 0.982.We also compute the network match values and its mode 0.962.It implies an excellent agreement between the GR template waveform and the observed data.We have seen that the reconstructed signal does remain almost identical even when the essential wavelets are selected using different posterior samples.Therefore, we consider the case of the maximum likelihood sample only to illustrate the signal is time-domain.Note that match values are computed between every posterior waveform and the corresponding wavelet-based reconstructed signal; we call this on-source match.
Our reconstruction method trades with sensitivity for identifying the essential wavelets and removing the noise from on-source data.Consequently, the technique cannot discern the early inspiral or late ringdown part of the CBC waveforms where the signals are weaker.In the time-frequency domain, the early inspiral part spreads over time, whereas the late-ringdown part of the signal over frequency direction.In Fig. 6, it is visible that the early inspiral part of the L1 signal (before 0.25 seconds) fades out.However, the H1 signal is still present and consistent with the template waveform because of the higher SNR in H1.At the same time, the late ringdown part of both the signals is disappeared.In a time-frequency frame, the ringdown part of a signal spreads out over frequency despite having a fixed frequency because of its exponential decay term, leading to a Lorentzian spread along the frequency direction.
Further, we perform the reconstruction analysis on the remaining events of GWTC-1 [55][56][57][58][59] and report the results in Table II.In order to understand the reconstruction efficiency depending on the SNR and chirp mass, we also reported them in the same Table .These values are obtained from the posterior samples by Bilby [51].The solid vertical line in Fig. 7 shows the 90% interval of the on-source network match values, and the solid circle marks their mode value.We found the best agreement with GR for the GW150914 event, the minimum for GW170608, albeit the observed SNR from the latter event was higher than the nominal threshold in both the detectors.The match value depends not only on SNR but also on the time-frequency area covered by the essential wavelets of the signal.In section III B, we have seen that the match value increases with injection chirp mass while the injection SNR is kept at a fixed value.GW170608 event has the lowest chirp mass in the catalog, for which its essential wavelets occupy the largest area over the time-frequency plane.
The vertical dashed line in Fig. 7 shows the 90% interval of expected match that is determined by injecting the maximum likelihood Bilby waveform (h * bif ) in many Gaussian noise realizations, similar study is shown in Fig. 4. We have reconstructed the signal from each realization and computed the match with h * bif .We see that there is significant overlap between the distribution of expected match and on-source match, except for GW151226 and GW170608.In particular, the distributions are far from each other for GW170608, for which a further study is worthy.We will focus on it in future work.
For GW170729, the match V1 data is negative, which is unexpected.The reconstructed signal and template waveform are in nearly antiphase for most of the posterior samples as shown in Fig. 8, leading to a negative match value.It could be due to the noise artifacts.The antiphase could occur due to poor constraint of event time in the detector frame.The total width of event time distribution for H1, L1, and V1 are 19, 18, and 63 ms, respectively.The spread for V1 data is significantly larger than the other two.It is impossible to draw an appropriate conclusion about whether the deviation presents or not when the SNR is very low.Thereby, we also compute the network match excluding the contribution of V1.The median network match for two detectors is 0.88, and the mode is 0.89.An extensive study on agreement between the CBC template waveform and the BayesWave reconstruction also found the similar results [60], where the contribution of V1 data was not considered for this study.We excluded the contribution of V1 for GW170729 to produce Fig. 7.As a further investigation, we project the H1 and L1 signal over V1 based on the sky location and polarization angle of parameter estimation samples.We follow these steps to obtain the projected signal: (a) perform the inverse transformation of whiten procedure over H1 and L1 reconstructed signals and apply timeshift based on the sky location to obtain these signals in geocentric coordinate, (b) obtain the two GW polarizations based on the sky location and polarization angle, and (c) project on the frame of V1 detector as shown in the bottom panel (blue dashed) of Fig. 8.We can see that the projected signal and CBC template waveform are in phase, and the median match is 0.64.It implies that the noise probably plays a role in V1 reconstructed signal being antiphase.

VI. CONCLUSION
This paper presents a wavelet-based semi-modeldependent method for reconstructing the gravitational wave signals produced from compact binary mergers.We have employed the framework of continuous wavelet transformation, where Morlet wavelets represent the signals.The semi-model-dependent approach determines the essential wavelets using the posterior samples from parameter estimation to reconstruct the signals from the data.
In general, the wavelets are constructed using an output of octave scale, which provides the tightest set of wavelets.Such wavelets yield a nearly orthogonal wavelet basis.In this paper, we have proposed a log-uniform scale for constructing the wavelets.Such wavelets are highly redundant, i.e., non-orthogonal wavelets.However, this new scale is more efficient for representing the linear chirp signals at high frequencies than the octave scale.As the wavelets are regarded as time localized bandpass filters, a reconstructed signal always contains a nominal amount of noise that passes through the essential wavelets.The amount of noise depends on the area covered by the essential wavelets over the time-frequency plane.We have shown that the essential wavelets with the log-uniform scale occupy a smaller space than the octave scale, which enables us to reconstruct the weak signals better.
We have conducted injection analysis to evaluate the reconstruction efficiency for the signals from binary black hole mergers by computing the match between injected waveform and the reconstructed signal.The reconstruction accuracy increases with the SNR such that the mismatch is approximately proportional to SNR squared, 1 − M 1/SNR 2 .Also, the reconstruction accuracy strongly depends on chirp mass.We have seen that the mismatch at a fixed SNR decreases as chirp mass increases.
We have demonstrated the ability to detect the devi-ation where the injected waveform is outside the region of space enclosed by the search template waveform manifold.We have performed the parameter estimation analysis using LALInference by injecting an electric BBH waveform, where the template waveforms are generated by a quasi-circular waveform model IMRPhenomXP.After that, we studied the match between the injected waveform and the reconstructed waveform.As reported in Table I, the LALInference waveform agrees less with the injected waveform than the wavelet reconstruction.Also, the LALInference waveform is out of phase over a bit of the region.In comparison, the amplitude and phase of the reconstructed waveform are nearly consistent with the injected waveform.
We have applied this wavelet-based reconstruction analysis to each binary black hole merger event in GWTC-1.We have seen a satisfactory agreement between the reconstructed signal and the estimated theoretical waveform.As the gravitational wave catalog grows, we expect the wavelet-based semi-model-dependent reconstruction method to provide a more precise view of the agreement between the observed data and the waveform model.
There are many avenues for future work and extensions of this method: combining the octave and log-uniform scales to have an optimal method for reconstructing the signals from compact binaries; applying this new method for reconstructing the signal from binary neutron star mergers; developing appropriate statistics for detecting the deviations and investigating the match study when instrumental glitch presents in on-source data; applying this method on GWTC-2 and GWTC-3 [4,5].

FIG. 2 .
FIG. 2. Performance of octave and log-uniform scale wavelets reconstruction for the linear chirp signals.The signals are generated using Eq.(15).The y-axis represents the normalized minimum squared error in the reconstructed signal.The solid line and dashed line correspond to the log-uniform scale and octave scale, respectively.The bottom and top x-axes represent the starting and ending frequency of the chirp, respectively.The legend represents a set of spectral loss parameters (R * ) as described in Eq. (12).We used those values to choose the essential wavelets.

FIG. 3 .
FIG.3.The noise averaged match ( M noise ) as a function of spectral loss parameter (R * ) is produced using a gravitational waveform of a nonspinning equal mass binary system with a chirp mass of 30M .The solid dot on the curve indicates its maxima.The right y-axis shows the relative number of essential wavelets ( NW ) as a function of R * .The unit value of NW corresponds to a 99.9% match between the original and reconstructed signal for the zero noise case.These results are produced assuming the aLIGO noise curve with a fixed lower cutoff frequency of 20Hz.The waveforms are generated using the SEOBNRv4 model.

FIG. 5 .
FIG. 5.An illustration of wavelet reconstruction (solid red line) and LALInference (cyan band) with IMRPhenomXP waveform model, obtained from an injection analysis using an NR eccentric waveform SXS:BBH:0323 picked from SXS catalog with a total mass of 60M .The dashed black line represents the injected waveform.For this computation, we have used the three detector configuration of Advanced LIGO and Advanced Virgo.

FIG. 6 .
FIG.6.Log-uniform scale wavelet reconstruction results of GW150914 event, obtained using the on-source data of H1 and L1 detectors, where the times (in seconds) are shown relative to a reference time 1126259462.0.The solid black line in the top-left panel shows the median reconstructed signal obtained from H1 data and the solid blue line in the bottom-left panel for L1 data.The orange band in these panels is produced using the CBC template waveforms from parameter estimation samples.We demonstrate the reconstructed signals and CBC template waveforms in units of the standard deviation of the noise, which implies the norm is SNR.The right panel shows the histogram of the match between every posterior waveform and the corresponding wavelet-based reconstructed signal.

FIG. 7 .
FIG. 7. Network match between the CBC template waveform by Bayesian parameter estimation and wavelet-based reconstructed waveform.The vertical solid line represents the 90% credible interval of the on-source match.The dot over each line is the mode of the on-source match distribution.The vertical dashed line indicates the 90% interval of the expected match under the assumption of Gaussian noise, which is similar to Fig. 4.

TABLE I .
Comparison between the LALInference template waveform and log-uniform wavelet reconstruction for an eccentric numerical relativity waveform SXS:BBH:0323 picked from SXS catalogue with a total mass of 60M .The numerical simulation is performed for a nonspinning, nearly equal-mass eccentric BBH system, where the measured reference eccentricity is e ref = 0.194 at a reference orbital frequency of M f0 = 0.0137.The quantity M denotes the median of the match values.

TABLE II .
List of match values of GWTC-1 events obtained by computing the match between the reconstructed waveform and Bayesian inference template waveform.We report the chirp mass and the SNR to indicate the efficiency of the reconstruction as demonstrated in section III B. We report the median (denoted by tilde) value of the distribution.Dashes (-) correspond to detector not included in the analysis.