A template-free approach for waveform extraction of gravitational wave events

We develop a general data-driven and template-free method for the extraction of event waveforms in the presence of background noise. Recent gravitational-wave observations provide one of the significant scientific areas requiring data analysis and waveform extraction capability. We use our method to find the waveforms for the reported events from the first, second, and third LIGO observation runs (O1, O2, and O3). Using the instantaneous frequencies derived by the Hilbert transform of the extracted waveforms, we provide the physical time delays between the arrivals of gravitational waves to the detectors.

We develop a general data-driven and template-free method for the extraction of event waveforms in the presence of background noise. Recent gravitational-wave observations provide one of the significant scientific areas requiring data analysis and waveform extraction capability. We use our method to find the waveforms for the reported events from the first, second, and third LIGO observation runs (O1, O2, and O3). Using the instantaneous frequencies derived by the Hilbert transform of the extracted waveforms, we provide the physical time delays between the arrivals of gravitational waves to the detectors.
An important prediction of General Relativity is the existence of Gravitational Waves (GWs), which act as ripples of space-time 1 . LIGO and VIRGO have observationally confirmed this prediction after a century. GWs detection has initiated a new era of research in gravity and has opened a new window to observe and study the Universe. One of the challenging issues in the GW events' extraction is the low amplitude of the ripples and the presence of significant background noise. Therefore, the techniques of signal detection and extraction of its waveform in presence of background noise become paramount.
The experimental efforts to detect GWs go back to Weber's cylindrical bar detector in the 1960s 10 . However, the efforts were not fruitful due to the weak effect of GWs on the vibration of Weber bar. On the other hand, the indirect evidence for the existence of GWs was provided in the astronomical data concerning the energy loss and hence the change in the rotational frequency of the Hulse-Taylor binary pulsar 11,13 . The energy loss matched well with the prediction of general relativity due to the gravitational waves' emission. Thanks to the great progress in the conceptual design of suspended mirrors and optomechanical technology in the Michelson-Morley type experiment of LIGO 12 , the direct detection of GW signals has become possible.
The GWs detection opens new observational windows involving astrophysical events 14 to the large scale structures of the Universe [31][32][33][34][35][36][37][38][39] . Moreover, GWs can place constraints on modified gravity theories 15,40,41 . Also, the detection of the GWs from merging binary neutron stars 16 opens a new area of multi-messenger astrophysics 17 . These observations can also enable us to investigate the classical and quantum properties of Black Holes (BH) 18,42 . There is even a possibility that GW candidates observed by LIGO could be due to the merger of primordial black holes which are candidates for dark matter 43,44 in the possible mass range of 20M ⊙ ≤ M pbh ≤ 100M ⊙ 43,44 .
The basis of GW-astronomy's diverse applications is detection, analysis, and interpretation of the GW data. There are several template-dependent and template-free methods, which have been developed in recent years for the detection of gravitational waves [45][46][47][48][49][50][51][52] . In this work, we develop a data-driven and template-free approach to extract event waveforms from gravitational-wave strain time series available from LIGO and Virgo detectors. Furthermore, we use the Hilbert-Huang transform and Hilbert spectrum to extract the time series's instantaneous frequency. The Hilbert spectrum enables us to find the physical time delays of the events in the GW detectors.

Results
Let us demonstrate our approach's applicability to detect the waveform of an event imposed into the background noise with different signal to noise ratio, using the introduced method. We consider a background time series n(t) (the data measured in the either H-or L-or V-detectors) as the stochastic background noise, where its mean value is subtracted. The next step is to superimpose a GW waveform template to the background noise. We take a typical time-series y(t) of one of the generated GW templates by LIGO (for instance, here we use the template for the event observed on 14 September 2015). For superimposed process to generate simulated data from the two time series of n(t) and y(t), we generate a new set of time series x α (t) = n(t) + α y(t) , with α ∈ [0, 1] . One can interpret the α factor as the strength of the GW template time series compared to the noisy background. In processing the simulated data, we apply our three steps introduced in previous section. We introduce x α (t) as the new time series (processed data) obtained from x α (t) , as well as ỹ(t) and ñ(t) obtained from the time series of y(t) and n(t) (with a different segment of the background data), respectively.
In Fig. 1, we plot the time series x α (t) (in which the template of GW150914, y(t), is imposed), in left upper panel, and extracted waveform, i.e. αỹ T (t) denoting the processed waveform, and ỹ(t) denoting the extracted waveform from our method with α = 0.005 in the left lower panel. The cross-correlation coefficient of extracted waveform with ỹ(t) is about 0.99 (right panel of Fig. 1). The strains are given in units of the standard deviation of the processed background noise, which provides us the statistical significance of an extracted waveform in each instant. For instance, the maximum amplitude of waveform for α = 0.005 , has about 12 σ statistical significance. We note that the one can interpret the parameter α as the distance of GW source from the detectors ("Materials and methods").
From the spectrum of the processed background noise (i.e. ñ(t) ), and the spectrum of the segment of data that includes the event (i.e. ỹ(t) ), one can define the signal to noise ratio (SNR) where ỹ(i, j) and ñ(i, j) represent the jth spectral component of the time segment i of the processed waveform ỹ(t) and the processed noise ñ(t) respectively 6 , and E[· · · ] is the mean operator. The size of time windows are about 0.2 s and we consider an averaged weighted spectrum σ 2 (i, j) over 10 windows with the same size to obtain E[|ñ(i, j)| 2 ] ("Materials and methods").
We define the integrated signal to noise ratio, ρ in effective frequency bands of 30-500 Hz, using thirty equal size-frequency bins 58 . This allows us to employ χ 2 time-frequency discriminator for gravitational wave detection, which enables us to reject the spurious events from the real ones 59 . The method applies to each data set and provides the probability P that the value of χ 2 would be obtained from the chirp signal 59 , see Table 1 and "Materials and methods". Our obtained integrated ρ depends on α as shown in Fig. 1c. As an example, we simulate the data with GW strains of α ≈ 5 × 10 −4 , where the detection algorithm provides the signal to noise ratio ρ H ≃ 58 with the cross-correlation coefficient ≈ 0.99 between the extracted waveform ỹ(t) and the processed one. The correlation coefficients of extracted waveforms and processed templates for different values of α are given in Fig. 1c. These results demonstrate our approach's high sensibility for the detection of events and the extraction of their waveforms in the time series. The estimated values for ρ using our approach for events reported by LIGO are ∼ 2.4 − 32.3 , see Table 1. www.nature.com/scientificreports/ To reduce the false trigger rate, for a given ρ * , we estimate the joint probability p([ρ H , ρ L ] > ρ * ) (i.e. ρ H > ρ * and ρ L > ρ * ) in shifted time-delay between the detectors, by assuming absence of correlated noise between detectors. This gives the rate of false alert for ρ * > 4 to be 1 alert per 200 days of running time 60 . The generalization of false alert rate for N-point (N-detector) joint probability p can be found in 60 . By analysing extensive background noise of about 637h of data, we find that the ρ has Rayleigh probability distribution function (see Fig. 2 and details are given in "Materials and methods").
We apply our method to the GW time series in the LIGO public database. In Fig. 3, we plot the extracted waveforms for GW150914, binary black-hole mergers in H and L detectors. The cross-correlation coefficient between the two extracted waveforms has the value ≃ −0.92 with a time-lag ≃ 7.3 +0. 3 −0.5 ms where L-detector received the (b) Extracted waveform and processed template are also shown. We find a very high similarity between extracted and processed waveforms. The correlation coefficient between two waveforms is ≃ 0.99 . The extracted waveform is in units of the standard deviation of the processed background noise. We applied our proposed steps to extract the injected template out of background noise. We called this an extracted waveform. Also we applied steps of our method only on the simulated template of GW150914 without adding it into background noise (which we called it processed template). (c) Signal to noise ratio (SNR) and cross correlation coefficients of extracted waveform with processed template for different α s. Error bars indicate the standard error of the mean averaged over 100 ensembles of background noise for each α . We did not preset the time of the injected waveform in 100 ensembles. In this paper SNR is defined as the ratio of extracted waveform and the mean value of processed background noise, ρ(i, j) = |ŷ(i, j)| 2 /E[|n(i, j)| 2 ] , summed over a certain frequency band. We note that due to the stochastic behavior of the background noise, different frequency bands in the spectrum will have effectively different SNRs in each time frame. Owing to the nature of our noise estimation approach, the noise will have a nonuniform effect on the SNR ratio of the extracted waveform specially when α is small 6 . However, as indicated in (c) when α increases, this effect will decrease and the relation between the SNR for the extracted waveform and the value of α becomes linear, with slops about 9300 and 5300 for L and H, respectively.

Figure 2.
Probability distribution function (PDF) of signal to noise ratio ρ for processed background noise. The PDF of signal to noise ratio ρ for the processed background noise as well as Rayleigh distribution (red) with variance 0.67. The total time duration of the analyzed time series is about 637h. The PDF is estimated from ∼ 74 × 10 6 signal to noise ratio ρ of processed background noise. We perform a similar analysis for the other fourteen events (ten events from the first and second observing runs of O1 and O2 and four events in the third observing run O3a). Our results find thirteen events out of the fifteen reported ones with cross-correlation coefficients between extracted waveforms obtained from our approach and simulated waveforms generated by IMRPhenomD, IMRPhenomPv3HM, IMRPhenomPv2 for GW190814, GW190521, and other remaining GWs, respectively, larger than 0.48. Their extracted waveforms are depicted in Fig. 4. The values of P and SNR of extracted waveforms are reported in Table 1. We excluded waveforms of the GW170817 and the GW190425 from Fig. 4, since our method was not able to extract clean chirp waveforms for these events.
Finally, we determine the time delays between the arrival of the signals to the detectors. In our method, the physical time delay is estimated from instantaneous frequencies derived from the Hilbert spectrum of the extracted waveforms (see "Materials and methods"). We determine f H (t) , f L (t) and f V (t) and calculate the mean increment in time lags (for instance for H and L) for the condition of |τ | < 10 ms , where The value of τ that minimize (global minimum) C(τ ) , is the physical time-delay. In our analysis, we didn't face the situation that C(τ ) of instantaneous frequencies cross twice or having two minima with the same depths.
In Figs. 3 and 4, we plot the instantaneous frequency derived from the Hilbert spectrum for all events in the first and second observing runs of O1, O2 and O3a ("Materials and methods").

Discussion
In summary, our aim in this work has been to introduce a new approach for the extraction of waveforms from GW events. Moreover, we estimate the physical time delays between the arrivals of gravitational waves at the GW detectors. The standard whitening procedure (in window size of 8 s) allows us to use the standard methods of false-alert rate estimation (see 60 and the references therein). The signal-to-noise ratio distribution is shown to be a Rayleigh distribution for the whitened data; see Fig. 5. We also provide the cross-correlation coefficients between the extracted waveforms in Fig. 4 and the simulated waveforms generated by PyCBC library for the events in the first, second, and third observing runs of O1, O2, and O3a, given in Fig. 6 62 . Moreover, the absolute value of the cross-correlation coefficients (R) of the extracted waveforms for H and L detectors are presented in Table 1.
The physical time delay between the amplitudes of the two detectors' signals is a complicated function of the orientation of the source, the orbital plane of the binary system, and the antenna patterns of the detectors. However, the Hilbert-Huang transform and Hilbert spectrum provide the instantaneous frequency of the extracted waveform independent of these geometric details. Therefore we are calling the estimated time delay as "physical time delay". In Methods, we generalized our approach to bivariate time series in which the extracted waveforms are estimated using the information provided by other detectors. We employed our generalized approach in the extraction of all waveforms in this study ("Materials and methods"). Although, we demonstrate the efficiency and applicability of our method in GW detection, however, this approach has the potential for extraction of event waveforms in many fields of science, ranging from neuroscience, physics to biology.

Materials and methods
Contemporary approaches to extract GWs. Contemporary approaches to extract GWs can be divided into two major categories.
Template-dependant methods. This seeks to find GW events from the statistical comparison between LHV datasets and a range of simulated GW templates in time or frequency domains. The most prominent candidate method in this regard is using match filtering technique that tries to match the time-series from detectors to a wide range of pre-generated GW template waveforms, under certain physically meaningful criteria, through an optimal match filter design 63 . This technique is well developed to search for a GW event included in the nonstationary stochastic background noise. Another technique introduced in Ref. 64 is to first find the GW events and then tries to extract a clean coalescence signal from noisy data. It is also accompanied with some presumptions related to a common gravitational waveform behavior. There is also another approach based on crosscorrelation of detectors in the time domain 65 . In this work, a simple measure consisting of the cross-correlation of detectors in short time intervals is introduced. The problem regarding this technique is that in contrast to the matched filter approach which is mainly based on cross-correlation in the frequency domain, here the correlation is measured in the time domain, thus demanding extra care. It means that the techniques used to improve the SNR should be carried out carefully to give rise to the possibly existing signal while suppressing the noise. This task alone, as is discussed and demonstrated in our work, needs to be implemented with much care since it is very sensitive to the true estimation of the non-stationary, non-Gaussian noise statistics. And also considering the abundance of instrumental glitches in the detectors with similar structures, the chances of detecting a glitch as a transient can be considerable. Possible correlations of the noise between detectors can be misleading too.  Table 1. We set the time to zero at the event GPS time reported by LIGO. We excluded waveforms of the GW170817 and the GW190425, since our method was not able to extract clean chirp waveforms for these events. We are using the information in bivariate time series in which the extracted waveform from one time series is estimated using the information provided by another detector. Therefore in any case that in given time series due to antenna pattern, local noises, etc. it possesses low SNR, this influence the extracted waveform from other time series. On the other hand these methods try to find the presence of GW events without any prior assumption about GW templates. In the wavelet reconstruction method, although it is true that wavelet reconstruction is able to obtain clean GW event-waveforms, it highly depends on the selection of the appropriate mother wavelet function which could limit the attempts for general noise suppression problems in real-world 63 The advantage of our proposed method compared to other common algorithms in the field is that our method can detect clean shape GW merger events without any prior information about the presence of the event. We start with the estimation of the noise variance in the Fourier domain from the signal after common pre-possessing steps. Then it updates the noise variance using an add-overlap method and tries to remove the noise based on this estimation. We also generalized this approach to calculate the noise statistics using a network of arrays (i.e. by updating the noise variance with a conditional joint probability bayesian estimator described in Eqs. (9)-(11), see below). Although this method works well to detect binary black hole mergers, we could not detect reliable results for binary neutron stars. The reason is due to the merger time duration in binary neutron stars that raise a challenge for the time constant that we consider for our add-overlap step. We showed that the results are consistent with events reported by LIGO (see cross-correlation coefficients presented in Fig. 6). Moreover, due to the stochastic characteristic of the background noise detected in different detectors and the different angles of arrival of events in detectors, we estimate time delays between detectors in the frequency domain using the Hilbert-Huang transform. This helps us find correct physical time delays for each event compared to time delays estimated from cross-correlation coefficients in the time domain. After deriving of instantaneous frequencies using the Hilbert-Huang transform, we are looking for chirp waveforms where the power increases monotonically in frequency and time.
At first we split the signal into 8 s windows to apply prepossessing steps. Then, we apply an add-overlap method on 0.06 s windows in the Fourier domain to estimate expected values of the noise during the noise estimation step. Finally, we reconstruct the signal in the Fourier domain from summation over all 0.06 s windows and estimate the Wiener gain function for the larger windows (8-s windows) to extract the noise from the raw data.
We did not obtain a clean extracted signal for binary neutron star mergers. Thus, we pick the frequency band 30 − 500Hz for our estimation to focus on binary black hole mergers.
Details of the event-waveform detection method. In this section, we explain in detail the method introduced in the main body of the manuscript.
Our event-waveform detection approach is based on the statistical comparison in the time-frequency domain of background noise and a segment of the data that may includes the GW event. The steps are: (i) First, we use a high-pass filter to filter out frequencies below 30 Hz . We also use notch filtering to eliminate certain high amplitude spectral lines, potentially disrupting the attempt to search for GW event-waveforms 54 .
(ii) We whiten the raw data that is a technique commonly used in astrophysics and cosmology. Whitening is equivalent to flattening the spectrum of a given signal, allowing all the bands to participate equally in the power spectrum of a given time series. We applied whitening over time windows of length 8 s.
(iii) Finally, we employ a template-independent method to suppress the background noise. This includes our generalized implementation of a two-step decision-directed noise reduction method proposed in 55 and Table 1. Results of our approach for detection of fifteen GW candidates. n/a time delays in the LIGO column means the corresponding time delay for that GW event is not reported by LIGO, at least we were not able to find it. n/a in the physical delay column means we didn't find a statistical meaningful time delay to report for that GW event. www.nature.com/scientificreports/ an iterative Wiener filtering as the gain function and a noise estimation method based on recursive averaging algorithms considering GW events presence uncertainty [6][7][8]55 . For convenience, we first explain the decision-directed noise reduction method and the a priori SNR used in Wiener filtering gain function described in Ref. 6 . Then we introduce the noise power spectrum density (PSD) estimation method used in our noise reduction algorithm. Let x(t) = y(t) + n(t) be the noisy time series that we obtain from pre-processing steps, which consists the GW waveform denoted by y(t) and the stochastic background noise denoted by n(t). We assume that y(t) and  . Cross correlation coefficients of the extracted waveforms for the events in the first, second and third observing runs O1, O2 and O3a by LIGO with waveforms generated by PyCBC library. Cross-correlation coefficients of the extracted waveforms for the event reported by LIGO with simulated waveforms generated by PyCBC library, where their waveforms are given in Fig. 4. For all events, A refers to the extracted waveform in H or L and B refers to simulated waveform for each detector 62 . Correlation coefficients between waveforms extracted from different detectors depend on two general parameters. One is the antenna pattern or the spatial angle between the detector and the GW which receives there. The other one is the SNR ratio of the recorded GW in each detector. www.nature.com/scientificreports/ n(t) are statistically independent, and both are zero mean time series. It is clear that in practice, we only have access to the noisy time series x(t) collected during LHV observations, and the presence of a GW event y(t) is unknown. However, we can estimate the noise variance during time frames in which no GW events are reported. We note that the nature of noise in advanced LHV time series is highly non-stationary and has a nonuniform effect so that different frequency bands have effectively different SNRs 54 . As a result, in order to be able to conduct a statistical comparison between consecutive time frames in each frequency bands, we continue our analysis in time-frequency domain using short-time Fourier transform (STFT) over windowed data segments as where index k indicates the corresponded detector in LIGO-VIRGO network array, X(i, j), Y(i, j), and N(i, j) represent the i-th time frame and j-th spectral components over time-frequency domain of the noisy time series, GW event-waveform, and background noise, respectively. Owing to the assumed independence of the GW waveform and the stochastic background noise, the periodogram of the noisy time series is approximately equal to the sum of the periodograms of GW signal and noise, respectively, that is: Once we obtain the magnitude and the spectrum of the desired GW waveform Y(i, j), we are able to reconstruct the final GW event-waveform in the time domain y(t), using inverse fast Fourier transform (iFFT). In practice, no direct solution for the spectral estimation exists. As the result, most common signal enhancement techniques require the evaluation of two parameters as the a posteriori and the a priori SNRs for the noisy components and GWs spectral components, respectively. Then, an estimate of Y(i, j) is subsequently obtained by applying a spectral gain function as: where ξ(i, j) and γ (i, j) are known to be a priori and a posteriori SNRs in i-th time frame and j-th frequency band, respectively, defined as: Here E[.] is the expectation value operator and q k shows the signal absence probability of the time series related to kth detector in time-frequency domain, which is set to be zero here, meaning that we assume the existence of signal in all short time-frequency segments that we are analyzing. The function G is the spectral gain function that should be applied to each short-time frame of the spectral component of the signal to obtain the spectral component of the clean signal. The choice of the gain function determines the gain behavior that determines the level of noise reduction in this method. We applied a Wiener gain function as described in Eq. (8) see below, and Ref. 6 .
We note that, one can't estimate Eq. (3) directly because we don't know the value of Y k (i, j) . However, we can estimate the a-priori SNR and Y k (i, j) respectively using the assumption that a-priori SNR for frame i can be estimated from Eq. (7) that is a decision-directed step (as known as first-order recursive function) to update ξ(i, j) from previous frame i − 1.
Substituting Eq. (10) into Eq. (9) yields: By tracking the noise power spectrum density (PSD), we are able to estimate the a priori and a posteriori and SNRs. In order to be able to bridge the broadest peak in the GW signal, we choose relatively wide enough window length − 0.06 s-for the estimation. However, the periodogram |X(i, j)| 2 of the noisy time series fluctuates very rapidly over time. As the result, we use a first-order recursive version of the γ (i, j) estimator called decisiondirected approach. Given that the a priori and the a posteriori SNRs can be estimated as: and where P[.] denotes the half-wave rectification and Ŷ (i − 1, j) is the estimated waveform spectrum at previous time frame in each frequency bands. We used half-wave rectification during the a priori and the a posteriori estimation using harmonic regeneration noise reduction-HRNR-method as proposed in Ref. 6 . Using HRNR to estimator SNRs has shown to be more accurate in reducing background noise while it helped us to prevent the distortion 6 . The control parameter 0 < < 1 is the weighting factor in decision-directed approach which determines the smoothness of the estimation. also is a smoothing factor that determines the smoothness of the a-priori estimator. Higher values of give smoother results. We choose = 0.95 as a trade-off between noise reduction and the event distortion. (1) Scientific Reports | (2021) 11:20507 | https://doi.org/10.1038/s41598-021-98821-z www.nature.com/scientificreports/ Then, using the Eq. (3), we are able to obtain the estimated clean GW waveform Ŷ (i, j) . The gain function we used in this paper is based on Wiener filtering as 6 : So far, we assumed that an estimation of the noise spectrum, N(i, j) is available and the noise PSD is given by σ 2 n = E[|N(i, j)| 2 ] . However, in general, we only have access to the noisy time series which yields having only the noisy signal spectrum X(i, j) to be known, while the PSDs of the GW waveform |Y (i, j)| 2 and the background noise σ 2 n are still to be unknown. Thus, both a posteriori SNR and a priori SNRs have to be estimated. In literature, there are different approaches to estimate the noise PSD σ 27 .
(i) Minimum statistics methods, assume that the noisy time series's power in each frequency band often decays to the noise power level. Thus, one can estimate the noise power spectrum by searching for minimum values in the noisy time series PSD, at each segment over the time-frequency domain. (ii) Histogram based methods, assume that the most frequent energy levels in the power spectrum distribution at each segment over the time-frequency domain correspond to the noise level. As a result, the noise PSD is estimated by finding maximum values in the power spectrum distribution. (iii) The time-recursive averaging methods are developed based on the assumption that the presence of noise has a nonuniform effect on the spectrum of a given time series so that different frequency bands have effectively different SNRs. Consequently, one can estimate the noise spectrum in each time-frequency segments whenever SNR is very low. This leads to a noise estimation approach based on a weighted average statistical comparison between the past estimates and the present noisy time series spectrum 7,55 . In this manner, the weighted average variable updates adaptively based on the noise and the GW event's absence or presence.
Here we use a generalized version of the Improved Minima Controlled Recursive Averaging (IMCRA) algorithm to obtain more accurate estimator due to the non-stationarity nature of the background noise 55 . In general, the presence or absence of a GW event or equivalently the detection of such event in jth frequency bin, in timefrequency domain of the data obtained from kth detector in LVH network array, is cast as a detection problem which yields in two following hypotheses: The noise PSD in kth detector can be estimated in terms of minimum mean square error as: where P(H j 0 |X k (i, j)) denotes the conditional probability of GW event being absent in time-frequency segment i, j given the noisy time series spectrum X k (i, j) collected in LHV network. Similarly, P(H j 1 |X k (i, j)) denotes the conditional probability of GW event being present in time-frequency segment i, j given the noisy time series spectrum X k (i, j) collected in LHV network. We can compute our conditional probabilities in Eq. (9), using Bayesian rules as follows: is the ratio of the a priori probabilities of GW event absence or presence, and is the likelihood ratio. Similarly for P(H j 1 |X k (i, j)) we have: Substituting (10) and (10) into (9) reads: When GW event is absent in frequency bin j in kth detector, the term E[σ 2 n,k (i, j)|H j 0 ] is approximately equals to the PSD of the noisy time series, |X k (i, j)| 2 . Alternatively, we can approximate E[σ 2 n,k (i, j)|H j 1 ] using the PSD of www.nature.com/scientificreports/ the noise in previous frames, σ 2 n,k (i − 1, j) 7,55 . Substituting these approximations into (11) gives the PSD estimate of the noise as: Equation (12) has the form of a time-recursive noise estimator where α k (i, j) = r k � k (i,j) 1+r k � k (i,j) is the smoothing factor 7 . Hence, under Gaussian assumption, � k (i, j) can be computed as: where ξ k (i, j) and γ k (i, j) are the a priori and the a posteriori that for kth detector can be estimated using decisiondirected estimator given by Eq. (6) and (7) 6,7 .
As we obtain spectral components of the desired GW waveform Y, we can reconstruct the denoised GW event-waveform in the time domain ỹ(t) using iFFT. We find that due to the non-stationarity of the time series under study, the noise frame of ∼ 10 s duration adjacent to the mainframe in which GW events received provide good results for the average weighted smoothing factor. This leads to an acceptable estimation of the averaged noise spectrum and also avoid the spectral artifacts that may appear locally in time. We also demonstrate the data in the unit of its standard deviations of a normalized window comprising the processed noise and the processed GW event's time windows, which in a way can be considered as an illustration of the signal contribution in noisy signal.
In summary, our method is based on the assumption that the power of the noisy signal in each frequency band decays to the power level of the noise as described in Refs. [6][7][8] . Therefore, we can estimate the noise level in each frequency band by tacking the minimum of the power in the noisy signal. This method was first developed to estimate the presence of speech signals under noisy conditions. Here, we use this idea and try to generalize it to obtain the clean GW event waveforms without having any prior information about the existence of the GW event in the analysed data. Moreover, we use the Wiener gain function as described in Ref. 6 and combine it with a noise estimator algorithm described in Ref. 55 to obtain better results for the power spectrum of the noise during our estimation. Our steps after pre-possessing include noise estimation as described in Eqs. (9)(10)(11)(12)(13). Then, we substitute the estimated expected value of |N(i, j)| 2 from this step into Eq. (7) and estimate the apriori SNR for ith time frame and jth frequency component. Finally, we substitute ξ(i, j) into Eq. (8) which is the Wiener gain function and estimate the extracted signal.
Interpretation of parameter α. The total strain indicated by x(t) = n(t) + αy(t) , the first term n(t) represents the noise of detector and y(t) is the signal. Since gravitational wave follows the wave equation as h µν = 0 , the solution of this equation depends on the distance from the source (in Minkowski Space) as h µν ∝ 1/r . So the coefficient α depends to r as α ∝ 1/r . Now, we can write the time dependent x(t) as x(t; r) = n(t) + (r 0 /r)y 0 (t) where r is the distance of GW source and r 0 is a given distance from the observer and y 0 (t) is the signal of a GW source at the distance of r 0 . In our analysis varying α within range of 10 −3 to 5 × 10 −3 means that we allow the distance of the GW source to change by fifty times. In other words, α 1 /α 2 = 50 then r 2 /r 1 = 1/50 . For instance, for the event "GW150914" located at distance of r 0 = 800Mpc , the time series at any arbitrary distance of r is related to α as r = 800Mpc/α . For the maximum peak of strain of 4 × 10 −20 in the template, using α = 0.005 means that distance of the source is located at few Gpc (taking into account the cosmological calculation) with the strain of 2 × 10 −22 .

Cross correlation coefficients of extracted waveforms with those inferred by LIGO.
To check the similarity of our extracted waveforms with those inferred by LIGO, we have estimated the Pearson correlation coefficients of the waveforms with the best fits among the waveforms simulated based on posterior BBH parameters reported by LIGO 30,58,66,67 . To this end, we have used the PyCBC package 62 to generate the time domain waveforms for each detector, using the posterior source files reported by LIGO for each event, including the parameters: primary and secondary masses, cartesian spin elements for each object, distance, inclination, declination, right ascension and polarization. The approximates used here are "IMRPhenomPv2" for all events except for GW190814 and GW190521, for which "IMRPhenomD" and "IMRPhenomPv3HM" is used, respectively. We note that by best fit, we are pointing to the one simulated waveform that returns the largest correlation coefficient with our extracted waveforms, which is chosen among thousands of posterior sets, typically. In Fig. 6, the absolute values of the correlation coefficient between all extracted waveforms and their associated simulated waveforms that match the extractions the best, are reported for both (L and H) detectors. We note that the correlation coefficients are evaluated in the 0.2s time-windows.
The χ 2 time-frequency discriminator. In this paper signal to noise ratio (SNR) is defined as the ratio of extracted waveform and the mean value of processed background noise, ρ(i, j) = |ŷ(i, j)| 2 /E[|n(i, j)| 2 ] , where the frequency band is divided in p-bins. We employ χ 2 time-frequency discriminator for gravitational wave detection as follows. In each frequency bin the normalized SNR is z i and we define z = p i=1 z i and Provided that the SNRs in each bin have Gaussian distribution, it is shown that is the incomplete gamma function 59 . High and low values of P χ 2 <χ 2 0 are (12)  30,58 . The duration of each capture was 4096 s and data points were sampled at 4KHz. The total time duration was 637h. A list of the name of time series can be found at Supp. Mat. To estimate SNRs, we first applied a highpass filter with a frequency pass > 30 Hz. We assumed the first 10 s of each time series as the initial noise to start our noise reduction algorithm. Finally, we estimated SNRs of each time series in p = 16 frequency bins between .
We find that our defined SNR in each bin follows from the Gaussian distribution, which is needed to employ the χ 2 time-frequency discriminator test (see Fig. 7).
Time delays between the arrivals of gravitational waves to the detectors.. We provide the value for the time delays between the arrival of gravitational waves to the detectors from analyzing the processed data. We use the Hilbert spectrum for the physically meaningful time-lag for the two detectors, i.e., < 10 ms for instance, for L and H.
Hilbert transform and Hilbert spectrum. Let us consider extracted waveforms in LIGO and Virgo detectors x(t) :=ỹ(t) and determine their local phases using for instance Hilbert transform or marked events method 3 . To determine the local "phase" of time series, we apply the Hilbert transform to process x as where P is the Cauchy principal value of the integral. We define analytical signal z(t) = x(t) + iy(t) = a(t)e iφ(t) , where a(t) = [x(t) 2 + y(t) 2 ] 1/2 and therefore local phase is given by φ(t) = tan −1 (y(t)/x(t)) . From the local phase one can calculate the local frequency f(t) via its time derivative f (t) = 1 2π d dt φ(t).
As an example, consider a chirp waveform as where the time dependent frequency is given by f (t) = f 0 k t , φ 0 is the initial phase, f 0 is the starting frequency (at t=0), and k is the rate of exponential change in frequency. In Fig. 8, we plot the waveform as well as the time dependent theoretical frequency and also the one estimated from Hilbert transform of discretized with dt = 10 −5 . The waveform is plotted for the parameters as k = 2.8 , f 0 = 25 and φ 0 = 0 . We repeat the Hilbert transformation for initial phase of φ 0 = π/3 , and notice that the time-dependent frequency obtained via Hilbert transformations with φ 0 = 0 and φ 0 = π/3 are the same.
x(t) := h(t) = sin φ 0 + 2πf 0 k t − 1 ln(k) , www.nature.com/scientificreports/ To estimate frequency-time relationship for the extracted waveform of ỹ(t) , we can apply Empirical Mode Decomposition (EMD) (see [49] 68 for details). The aim is to decompose the original signal into a hierarchy of Intrinsic Mode Functions (IMF) that separate the signal's different frequency components by counting the maxima and zero crossings. In Fig. 9, we plot the IMFs of the extracted waveform ỹ(t) for the event GW150914 in H-detector, using the EMD method.
To estimate local frequency or frequency-time relationship, we then apply Hilbert transform to each IMF components of the extracted waveform ỹ(t) , excluding the final residuals, where one can express it as the real part ( R ) of the sum of the Hilbert transform of all the IMF components 68 , where a j (t) and ω j (t) are the amplitude and the frequency of the jth IMF component, respectively. Thus, the amplitude is a function of time and frequency. The frequency-time relationship of the amplitude is known as the ỹ(t) = R