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.

Introduction

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 method of extracting signals from a given time series has important applications in many fields of science, ranging from neuroscience to astrophysics[2-8]. We are interested in the application of this approach for the detection of GWs[1,9-20].

To this end, since the first direct detection of GWs by the Laser Interferometer Gravitational-wave Observatory (LIGO) on 14 September 2015, more than fifty-two “confirmed events” (from the first, second, and third observation runs O1, O2, O3a, and O3b ) have also been detected[21-24,24-30] In this paper, we have analyzed fifteen publicly available GW events in GWTC-1 and GWTC-2 catalogues from O1, O2, and O3a runs (see Table 1).

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.

Event	Delay (ms)		P %	SNR	R
	LIGO	Physical delay
GW150914[46]	6.9-0.4+0.5 – L-first	7.3-0.5+0.3– L-first	P_H=0.99, P_L=0.79,	ρ_H=12.1,ρ_L=14.3	0.92-0.007+0.007
GW151012[61]	0.6-0.6+0.6	0.5-0.3+0.5 – L-first	P_H=0.78, P_L=0.99,	ρ_H=19.6,ρ_L=10.2	0.68-0.005+0.005
GW151226[21]	1.1-0.3+0.3 – L-first	1.2-0.5+0.7 – L-first	P_H=0.99, P_L=0.97,	ρ_H=11.9,ρ_L=10.6	0.33-0.008+0.008
GW170104[22]	3.0-0.5+0.4 – H-first	3.2-0.2+0.5 – H-first	P_H=0.83, P_L=0.81,	ρ_H=21.8,ρ_L=28.5	0.86-0.013+0.013
GW170608[23]	7 – H-first	6.8-0.5+0.2 – H-first	P_H=0.99, P_L=0.99,	ρ_H=5.9,ρ_L=3.0	0.72-0.004+0.004
GW170729[24]	n/a	1.8-0.9+1.0 – L-first	P_H=0.99, P_L=0.99,	ρ_H=13.7,ρ_L=17.7	0.60-0.023+0.023
GW170809[24]	n/a	9.5-0.5+0.5 – L-first	P_H=0.99, P_L=0.88,	ρ_H=7.9,ρ_L=9.8	0.63-0.023+0.023
GW170814[25]	8 – L-first	7.8-0.5+0.8 – L-first	P_H=0.99, P_L=0.99,	ρ_H=8.8,ρ_L=5.4	0.76-0.011+0.011
GW170817[16]	n/a	n/a	P_H=10^-3, P_L=0.92,	ρ_H=32.3,ρ_L=2.4	0.01-0.10+0.10
GW170818[24]	n/a	4.8-0.8+.5 – L-first	P_H=0.99, P_L=0.70,	ρ_H=4.6,ρ_L=21.0	0.70-0.023+0.023
GW170823[24]	n/a	1.5-0.5+0.7 – H-first	P_H=0.71, P_L=0.99,	ρ_H=21.3,ρ_L=9.8	0.55-0.022+0.022
GW190412[26]	n/a	4.0-0.2+0.3 – L-first	P_H=0.98, P_L=0.93,	ρ_H=7.0,ρ_L=21.8	0.79-0.002+0.002
GW190425[27]	n/a	n/a	P_V=10^-3, P_L=0.50,	ρ_V=32.1,ρ_L=6.6	0.01-0.13+0.13
GW190521[29]	n/a	2.0-0.4+0.3 – L-first	P_H=0.99, P_L=0.99,	ρ_H=6.5,ρ_L=4.8	0.95-0.002+0.002
GW190814[28]	n/a	3.4-0.6+0.5 – H-first	P_H=0.94, P_L=0.99,	ρ_H=14.3,ρ_L=8.4	0.49-0.062+0.062

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-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 [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-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.

A template-free approach

The direct detection based on laser interferometry by LIGO relies on amplitudes and phases of the strain data obtained from two identical setups at Hanford and Livingston (denoted henceforth by “H” and “L”, respectively) and Virgo detectors (“V”) with a standard template-dependent analysis[46]. It is based on the correlation of a set of limited GW template bank and the empirical data, utilizing the optimally matched filter method[53]. From matching with templates, one obtains the source properties, and from the time delay between the arrival of the gravitational wave to the detectors, one obtains limited information on the events’ angular position, which has led to the emission of the GW.

Our event-waveform detection approach is based on the statistical comparison in the time-frequency domain of background noise and the portion of data that includes the event (“Materials and methods”). The method consists two pre-processing steps and a spectral subtraction noise reduction algorithm that performs well for extracting any non-noise features from a noisy time series. The steps are: (i) First, we use a high-pass filter to filter out frequencies below . We also use notch filtering to eliminate certain high amplitude spectral lines, potentially disrupting the attempt to search for GW event-waveforms. These narrow resonances (high amplitude spectral lines) are caused by different sources, including harmonics of electrical power, violin modes due to mechanical resonances of the mirror suspensions, and/or spectral lines produced by calibration[54]. (ii) We whiten the raw data, which 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. This helped the frequency content in data appears equally, and even the smallest contributions could be observed. Time windows of length 8 s was selected because it was relatively large enough compared to GW events duration to avoid undesired boundary effects. While it was comparatively small enough to increase the weight of all frequency bands, including the frequency range of GW events, to have equal weight in the power spectrum[54]. (iii) Finally, we employ a generalized template-independent method to suppress the background noise, which is a combination of a two-step decision-directed noise reduction method[6,6,8] using Wiener filtering as its gain function and a noise estimation method based on recursive averaging algorithms analyzing GW events in the presence of uncertainty[7,55] (“Materials and methods”). We refer to the final data as processed data. (iv) From analyzing of the processed data, we provide values for the physical time delays between the arrival of gravitational waves to the detectors, using instantaneous frequencies of extracted waveforms, derived from Hilbert spectrum[3,56,57] (“Materials and methods”).

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 , with . 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 as the new time series (processed data) obtained from , as well as and 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 (in which the template of GW150914, y(t), is imposed), in left upper panel, and extracted waveform, i.e. denoting the processed waveform, and denoting the extracted waveform from our method with in the left lower panel. The cross-correlation coefficient of extracted waveform with 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 , has about statistical significance. We note that the one can interpret the parameter as the distance of GW source from the detectors (“Materials and methods”).

Figure 1

Extraction of injected signal into background noise with different signal to noise ratio (SNR). (a) Template of event on 14 September 2015 added to the raw data measured in the H-detector with . (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 . 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, , 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.

From the spectrum of the processed background noise (i.e. ), and the spectrum of the segment of data that includes the event (i.e. ), one can define the signal to noise ratio (SNR) as , where and represent the jth spectral component of the time segment i of the processed waveform and the processed noise respectively[6], and is the mean operator. The size of time windows are about 0.2 s and we consider an averaged weighted spectrum over 10 windows with the same size to obtain (“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 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 that the value of 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 , where the detection algorithm provides the signal to noise ratio with the cross-correlation coefficient between the extracted waveform 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 , see Table 1.

To reduce the false trigger rate, for a given , we estimate the joint probability (i.e. and ) in shifted time-delay between the detectors, by assuming absence of correlated noise between detectors. This gives the rate of false alert for to be 1 alert per 200 days of running time[60]. The generalization of false alert rate for N-point (N-detector) joint probability 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”).

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 signal to noise ratio of processed background noise. Each SNR is estimated for the 32 frequency bins between Hz for a 10 s segments of time series in O2 and O3a run (see Supp. Info. for the GPS time of analyzed data). Inset the PDF is plotted on a log-linear scale compared to the right-tail of estimated and Rayleigh PDFs.

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 with a time-lag where L-detector received the signal first. Also, the cross-correlation coefficients of the extracted waveforms in H and L detectors from our method with the templates of event GW150914, reported by LIGO, are and , respectively.

Figure 3

Extracted waveforms for the event GW150914 in H and L and instantaneous frequency derived by Hilbert spectrum. Extracted waveforms for the event GW150914 in units of the standard deviation of the processed background noise are given. The (Yellow) and (Orange ) confidence intervals were obtained via bootstrapping. The cross-correlation coefficient between two extracted waveform has value in time lag L-first. Cross-correlation coefficients of extracted waveforms in H and L with processed template of event GW150914 are and , respectively. For extracted waveforms in H and L, instantaneous frequencies derived by the Hilbert spectrum is plotted. Color bar indicating the instantaneous power of the signal at each point in the waveforms. We set the time to zero at the event GPS time reported by LIGO.

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 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.

Figure 4

Extracted waveforms of binary black-hole mergers in observing runs O1, O2 and O3a. Extracted waveforms of binary black-hole mergers in O1, O2, and O3a observing runs. They are labeled by the event name and detector that received the signal. The strains are given in units of the standard deviation of the processed background noise. For extracted waveforms in H, L, and V, the Hilbert spectrum’s instantaneous frequencies are plotted. The cross-correlation coefficients (R) of extracted waveforms are given in 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.

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 , and and calculate the mean increment in time lags (for instance for H and L) for the condition of , where . The value of that minimize (global minimum) , is the physical time-delay. In our analysis, we didn’t face the situation that 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.

Figure 5

Probability distribution functions of normalized background and its whitened noises before the event 29 July 2017. In lower panel, probability distribution functions (PDF) of normalized background noise estimated from data points is plotted in log-linear scale. In upper panel, the PDF of its whitened time series is depicted which has a Gaussian.

Figure 6

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.

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 non-stationary 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 cross-correlation 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.

Template-free methods

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 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 . 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 an iterative Wiener filtering as the gain function and a noise estimation method based on recursive averaging algorithms considering GW events presence uncertainty[6-8,55].

For convenience, we first explain the decision-directed noise reduction method and the 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 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 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 aswhere 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 and 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 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 . However, we can estimate the a-priori SNR and 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 from previous frame .

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 of the noisy time series fluctuates very rapidly over time. As the result, we use a first-order recursive version of the estimator called decision-directed approach. Given that the a priori and the a posteriori SNRs can be estimated as:andwhere P[.] denotes the half-wave rectification and 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 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 as a trade-off between noise reduction and the event distortion.

Then, using the Eq. (3), we are able to obtain the estimated clean GW waveform . 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 . 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 and the background noise 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 [7].

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.

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.

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 time-frequency 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 denotes the conditional probability of GW event being absent in time-frequency segment i, j given the noisy time series spectrum collected in LHV network. Similarly, denotes the conditional probability of GW event being present in time-frequency segment i, j given the noisy time series spectrum collected in LHV network. We can compute our conditional probabilities in Eq. (9), using Bayesian rules as follows:where is the ratio of the probabilities of GW event absence or presence, and is the likelihood ratio. Similarly for we have:

Substituting (10) and (10) into (9) reads:

When GW event is absent in frequency bin j in detector, the term is approximately equals to the PSD of the noisy time series, . Alternatively, we can approximate using the PSD of the noise in previous frames, [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 is the smoothing factor[7]. Hence, under Gaussian assumption, can be computed as:where and are the and the that for detector can be estimated using decision-directed 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 using iFFT. We find that due to the non-stationarity of the time series under study, the noise frame of 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-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–13). Then, we substitute the estimated expected value of from this step into Eq. (7) and estimate the apriori SNR for time frame and frequency component. Finally, we substitute into Eq. (8) which is the Wiener gain function and estimate the extracted signal.

Interpretation of parameter

The total strain indicated by , the first term n(t) represents the noise of detector and y(t) is the signal. Since gravitational wave follows the wave equation as , the solution of this equation depends on the distance from the source (in Minkowski Space) as . So the coefficient depends to r as . Now, we can write the time dependent x(t) as where r is the distance of GW source and is a given distance from the observer and is the signal of a GW source at the distance of . In our analysis varying within range of to means that we allow the distance of the GW source to change by fifty times. In other words, then . For instance, for the event ”GW150914” located at distance of , the time series at any arbitrary distance of r is related to as . For the maximum peak of strain of in the template, using means that distance of the source is located at few Gpc (taking into account the cosmological calculation) with the strain of .

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 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, , where the frequency band is divided in p-bins. We employ time-frequency discriminator for gravitational wave detection as follows. In each frequency bin the normalized SNR is and we define and . Provided that the SNRs in each bin have Gaussian distribution, it is shown that , where is the incomplete gamma function[59]. High and low values of are indicators of chirp and spurious signals. Intuitively, this stems from the fact that in chirp signal in the chosen frequency band, all SNRs have finite values and their dispersion are small, however in spurious signal, some bins have high SNR and in other bins have smaller, then there is the high value of dispersion.

To check the Gaussianity of SNRs in each frequency bin, we estimate the SNR for a random set of 560 captures from Hanford detectors during the O2 (from GPS Time 1164556817 to GPS Time 1187733618) and O3a (from GPS Time 1238166018 to GPS Time 1253977218) run[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 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 frequency bins between [30 - 512 Hz].

We find that our defined SNR in each bin follows from the Gaussian distribution, which is needed to employ the time-frequency discriminator test (see Fig. 7).

Figure 7

The results show PDF of normalized SNRs for each frequency bin. The total time duration of analysed time series is about 637 h. In each subgraph a Gaussian PDF with unit variance is plotted for comparison.

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., for instance, for L and H.

Hilbert transform and Hilbert spectrum

Let us consider extracted waveforms in LIGO and Virgo detectors 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 aswhere P is the Cauchy principal value of the integral. We define analytical signal , where and therefore local phase is given by . From the local phase one can calculate the local frequency f(t) via its time derivative .

As an example, consider a chirp waveform aswhere the time dependent frequency is given by , is the initial phase, 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 . The waveform is plotted for the parameters as , and . We repeat the Hilbert transformation for initial phase of , and notice that the time-dependent frequency obtained via Hilbert transformations with and are the same.

Figure 8

A chirp waveform and its Hilbert transform. A sinusoidal exponential chirp waveform and its theoretical time-dependent frequency as well as estimated frequency from Hilbert transform of discretized waveform. The waveform is plotted for the parameters are , with , and .

To estimate frequency-time relationship for the extracted waveform of , 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 for the event GW150914 in H-detector, using the EMD method.

Figure 9

Intrinsic mode functions of extracted waveform of event GW150914 in H. The intrinsic mode functions (IMF1 to IMF4) and residual (res.), for the extracted waveform of the event GW150914 in H from empirical mode decomposition.

To estimate local frequency or frequency-time relationship, we then apply Hilbert transform to each IMF components of the extracted waveform , excluding the final residuals, where one can express it as the real part () of the sum of the Hilbert transform of all the IMF components[68],where and 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 Hilbert spectrum. The inherent characteristics of a nonlinear and/or non-stationary waveform can be identified from the Hilbert spectrum. Local frequencies of all extracted waveforms in this work are estimated from the Hilbert spectrum. One can verify that the Hilbert spectrum and the local time-frequency, such as a sinusoidal, exponential chirp waveform, provide the same frequency-time relationship with different initial phases .

Supplementary Information 1.