Autoregressive Search of Gravitational Waves: Denoising

Abstract

Because of the small strain amplitudes of gravitational-wave (GW) signals, unveiling them in the presence of detector/environmental noise is challenging. For visualizing the signals and extracting its waveform for a comparison with theoretical prediction, a frequency-domain whitening process is commonly adopted for filtering the data. In this work, we propose an alternative template-free framework based on autoregressive modeling for denoising the GW data and extracting the waveform. We have tested our framework on extracting the injected signals from the simulated data as well as a series of known compact binary coalescence (CBC) events from the LIGO data. Comparing with the conventional whitening procedure, our methodology generally yields improved cross-correlation and reduced root mean square errors with respect to the signal model.

I Introduction

The existence of gravitational wave (GW) is one of the most remarkable predictions of general relativity (GR)[1, 2]. GW is a tidal acceleration that propagates in spacetime at the speed of light. According to the Einstein field equation, it requires stress at an order of $c^{2}/8\pi G\sim 10^{43}$ N m${}^{-2}$ to produce a unit of curvature. Therefore, the amplitude of GW is expected to be very small. And it requires catastrophic phenomena involving compact objects to produce such tiny ripples in the spacetime (e.g. binary black hole mergers).

It is the small amplitude of GW that makes the detection challenging. The first compelling evidence for the existence of GW came indirectly from the long-term pulsar timing of the Hulse-Taylor binary PSR B1913+16 [3]. The behavior of this binary (e.g. decay of orbital period) is fully consistent with the prediction by GR as the system loses its orbital energy in GW.

Thanks to the improved sensitivity, on 14 September 2015, the advanced Laser Interferometer Gravitational-wave Observatory (LIGO) has directly detected a GW event, GW150914, from a binary black hole (BBH) coalescence for the first time [4]. This has opened the possibility of exploring our Universe without limiting to the window of electromagnetic radiation. Two years later, the era of multi-messenger astronomy was highlighted by the discovery of GW event GW170817 resulted from the merger of two neutron stars (NSs) [5, 6, 7], which was found to be associated with the $\gamma-$ray burst GRB 170817A [8]. This marks the first case that both GW and electromagnetic radiation were detected from the same astrophysical object.

Currently, in the Gravitational Wave Transient Catalog (GWTC) maintained by LIGO/Virgo/KAGRA collaboration [10, 11, 12, 13], there are 93 GW transient events so far have been confidently detected (i.e. probability of origin from an astrophysical source $p_{\rm astro}>0.5$). These include 89 from BBH coalescence, 2 from NS-NS mergers, and 2 from BH-NS mergers. Apart from these confident events, there are $>20$ marginal candidates.

For further advancing GW astronomy, while enhancing the instrumental sensitivity is vitally important [e.g. 14], progress can also be achieved by improving the methodology of data processing and analysis. Currently, the standard search method for CBC events in the GW community is matched filtering [cf. 15], which is done by cross-correlating a template of known waveform and the interferometer output at different time delays to produce a filtered output. With the signal-to-noise ratio (SNR) as the ratio of the value of the filtered output to the corresponding value root mean square value for the noise, it can be proved that a matched filter comprises the ratio of the template of the actual waveform to the spectral noise density of the interferometer can optimize the SNR under several assumptions [16].

Although the technique of matched filtering has unveiled a considerable population of GW events as aforementioned, it has a number of limitations. For the matched filter to have optimal performance, the data have to fulfill the assumptions of wide-sense stationarity (WSS) and zero-means, which are generally not satisfied in the raw interferometric data. And most importantly, the construction of matched filter requires knowledge of waveform for the expected signal. However, the forms of GW signal from many possible sources are poorly modeled (e.g. highly eccentric BH binaries) or even unknown (e.g. fast radio bursts). In such cases, matched filtering technique cannot be employed. Even for the cases that the waveform can be determined such as circular BH binaries, this technique still requires the construction of a large template bank to cover a sufficiently large parameter space. A search over this extensive template bank by brute-force is computationally expensive.

Furthermore, even though the technique of matched filtering is capable to detect the GW signals from CBC events with known waveform, it does not enable one to visualize the signal directly. For visualizing the GW signal from the data and extracting its waveform for a comparison with the prediction by numerical relativity [e.g. Fig. 1 in 4], one must filter the raw time series with a bandpass filter for removing the data out of the detectors’ most sensitive frequency band as well as apply the frequency-domain whitening process for suppressing the colored noises at low frequencies and the spectral lines resulted from instrumental/background effects. Frequency-domain whitening is a procedure to equalize the spectrum through dividing the Fourier coefficients by the estimate of the amplitude spectral density of the noise. While this is considered as a standard procedure and has been adopted for noise suppression in many works [e.g. 17, 18], it is still important to explore alternative techniques for de-noising and compare their performance with that of conventional whitening filter. For example, Tsukada et al. have proposed a time-domain whitening filter for optimizing latency in the CBC data analysis pipeline [19].

In this paper, we explore the feasibility of a template-free method based on autoregressive modeling in filtering GW data. In Section II, we provide an overview of the methodology. In Section III, we will demonstrate the feasibility of our framework by a series of experiments. And we will summarize our results and provide an outlook for further development in Section IV.

II Methodology

II.1 Autoregression

Time series data that we acquire in nature can be affected by various random processes and exhibit stochastic behaviors. Due to the high sensitivity of GW detectors, the raw data are typically corrupted with the various kinds of noise [e.g. 20], in which the assumptions for the matched filter to attain the optimal performance such as stationarity are generally not fulfilled. In our proposed framework, we adopted an autoregressive (AR) approach in developing an efficient time-domain noise filtering scheme without any a priori knowledge on the noise.

In a recent astronomical application of AR modeling, Caceres et al. have developed a methodology of the autoregressive planet search (ARPS) for treating a wide variety of stochastic processes so as to improve the search of transit signals by exoplanets in the residuals after noise reduction [21]. In exoplanet search, people are looking for small dips in the light curve resulting from a transit submerged by the much larger brightness variability of the parent stars. And the aperiodic colored noise in the photometric data is notoriously difficult to treat [22]. With a procedure based on AR, Caceres et al. have demonstrated the stellar variability can be identified and removed.

We notice that the aforementioned challenge is shared by the GW astronomy, namely searching for the small strain amplitude of GW signals in the presence of instrumental/environmental noise with amplitude orders of magnitude larger. This comparison has motivated us to explore whether the AR technique can also be applied in extracting GW signals.

AR modeling can be applied to any dynamical system whose status in the present time has a dependence on its past status (i.e. autocorrelated behavior). The simplest model AR($p$) can be built by regression with the estimate at time $t$, $\hat{x}_{t}$, being modeled by the linear combination of past values $x_{t-i}$ plus a random noise term:

$$\hat{x}_{t}=\sum^{p}_{i=1}a_{i}x_{t-i}+\epsilon_{t}$$ (1)

where $p$ is the order of AR model (i.e. the number of lags in the model), $a_{i}$ are the model parameters, $x_{t-i}$ is the $i$-th past data, and $\epsilon_{t}$ is the noise term distributed as a Gaussian with zero mean and unknown variance.

In the application of ARPS, the best-fit AR model can be treated as a good estimator of stochastic noise. By subtracting the model from the raw data, the residuals can be obtained as follow:

$$x_{r}=x-\hat{x}$$ (2)

where $\hat{x}$ is the best-fit AR model on data $x$. From $x_{r}$, we can investigate whether there is any astrophysically-interesting signals can be extracted [see 21, 23].

II.2 Autoregressive integrated moving-average model (ARIMA)

However, for an AR model to provide a legitimate description of the data, the time series is assumed to be stationary. For the time series with systematic trends, such data cannot be treated by the AR model. For converting a non-stationary time series into a stationary one, the differencing operation is found to be efficient (e.g. $x_{t}^{{}^{\prime}}=x_{t}-x_{t-1}$ where $x_{t}^{{}^{\prime}}$ is the differenced series obtained from the change between consecutive values in the original time series). With the backshift operator $B$ defined as $Bx_{t}=x_{t-1}$, the aforementioned process can be described as $(1-B)x_{t}$. This is known as first-order differencing. To generalize the process to a higher orders, the operation can be modified as

$$x_{d}=(1-B)^{d}x_{t}$$ (3)

which is commonly referred to as Integrated process (I) of $d^{\rm th}$ order. The output $x_{d}$ can be modeled as a stationary time series.

While AR model uses past values in a time series to predict the current value, a Moving Average model of order $q$, MA($q$), predicts the current value by a linear combination of past error terms

$$\hat{x}_{t}=\sum^{q}_{j=1}b_{i}\epsilon_{t-j}+\epsilon_{t}$$ (4)

where $\epsilon_{t-j}$ is the error term for the $j-$th time step in the time series $x_{t}$ and $b_{i}$ are the model parameter.

Combining Equation 1, Equation 3 and, Equation 4, an AutoRegressive Integrated Moving-Average model (ARIMA) can be constructed as:

$$(1-B)^{d}x_{t}=\sum^{p}_{i=1}a_{i}x_{t-i}+\sum^{q}_{j=1}b_{i}\epsilon_{t-j}+\epsilon_{t}$$ (5)

with $a_{i}$ and $b_{i}$ estimated simultaneously for the whole time series by any optimization method (e.g. maximum likelihood). On the other hand, the orders of the model (e.g. $p$, $q$, $d$) can be determined by the procedures of model selection.

By applying ARIMA modeling to the light curves of 156717 stars as observed by NASA’s Kepler satellite, Caceres et al. shows that the brightness fluctuations of the parent stars can be effectively reduced. Subsequent searches of transit-like signals from the ARIMA residuals resulted in a recovery of a significant fraction of confirmed exoplanets [23].

We started our experiment by testing whether the existing code auto.arima from the R package forecast is capable to fit the ARIMA model to GW data with the orders and the parameters of the model estimated automatically. auto.arima was used in constructing ARPS pipeline [23].

We have tested whether auto.arima can extract the waveform of GW150914 from LIGO data obtained by both detectors in Hanford (H) and Livingston (L) with a sampling frequency of 4 kHz. The data were downloaded from the event portal managed by the GW Open Science Center (GWOSC). However, since the strain amplitude is at an order of $h\lesssim 10^{-18}$, the modeling apparently suffered from an underflow problem.

In order to circumvent the underflow, we have normalized the data to the order of unity. In this case, auto.arima yields the orders of (p,q,d)=(4,0,1). By subtracting the resultant model from the original data, we have obtained the residuals. However, even with the aid of a low-pass filter, we are unable to identify any waveform which resembles that of GW150914 from the residuals. In view of such an undesirable behavior, instead of employing auto.arima as in [21, 23], we are going to develop our own algorithm which is more suitable for reducing noise in GW data.

Figure 1: The structure of our proposed framework seqARIMA.

Figure 2: Illustrations for the effects of each stage in seqARIMA by using the LIGO-H data of GW150914.

In our proposed framework, we break the noise reduction process into a sequence of procedures as shown in Figure 1. Hereafter we refer it as sequential ARIMA model (seqARIMA), which consists of four stages: integrated process, autoregressive process, moving-average process, and bandpass filtering.

In this section, we take the LIGO data of GW150914 (Hanford, 32 s, 4 kHz sampling) to demonstrate the performance of seqARIMA. The effects of each stage in our procedure are illustrated in Figure 2 and described in the following subsections (i.e. subsection II.3-subsection II.6).

II.3 Integrated process

As a first step of our proposed framework, integrated process plays an essential role for ensuring the stationarity of a given GW time series.

For demonstrating the procedure, we start with the raw LIGO-H data of GW150914 (top panel in Figure 2) and with the parameter corresponding to the order of differencing initialized as d=0. We employ the Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test [26], which is a standard test for stationarity [27], to examine whether the raw data exhibits trend and non-zero mean level. Figure 2 clearly shows that the raw data is non-stationary.

Instead of ensuring global stationarity on the entire input time series, we consider local or segmented stationarity within time windows which are comparable to the signal duration from the typical BBH coalescence. A segment within a time series refers to a sequence of data points collected or recorded at a given time interval. Most time series segmentation algorithms can be classified into three primary categories: sliding windows, top-down, and bottom-up approaches [28]. Lovrić et al. have demonstrated that the process of segmenting time series into a limited number of homogeneous segments aids in the extraction of time segments with similar observations [29]. These techniques process the input time series and return a piecewise linear representation (PLR). Our method, however, divides time series into non-overlapping times series of equal length for the best performance. Since the total length of the time series is much longer (32 s) than the duration of the typical BBH coalescence signal (i.e. $\lesssim$ 0.5 seconds), the whole time series is divided into 64 segments and the KPSS test is applied on each segment with the length 0.5 s (i.e. $t_{\rm s}=0.5$). If the $p$-values from KPSS test on all those segments are greater than or equal to our predefined threshold ($p$-value = 0.1), the given data is determined as a stationary time series. The threshold is chosen to be larger than the conventional value of $p$-value = 0.05 so as to reduce the false negatives.

Otherwise, if there is any segment exhibits non-stationary behavior, we modify the parameter $d$ as $d+1$ and apply differencing by Equation 3. Such process will be iterated until $x_{d}$ satisfies the stationary condition by passing the KPSS test (See algorithm 1). In the second panel of Figure 2, we show an optimal differencing model for our test data with $d=2$.

Since the result of hypothesis test can be influenced by the volume of data used, we further test the robustness of algorithm 1 by running KPSS tess on different $t_{s}$. In the aforementioned experiment, we took $t_{s}=0.5$ s which gives $\sim 2000$ data points for a sampling rate of 4 kHz. For investigating the possible impact of non-stationarity detection by the length of data segment, we have re-run algorithm 1 on the same data by varying $t_{s}$ from 0.25 s to 0.75 s. And we found that all cases yield the same optimal differencing model with $d=2$. In view of this, we conclude that the results from KPSS test and hence algorithm 1 is robust and our adopted segment length of $t_{s}=0.5$ s is sufficient.

Input : $x_{t}$, the input time series,

$t_{\rm t}$, the time length of $x_{t}$,

$t_{\rm s}=0.5$, the length of each segment,

$c_{p}=0.1$, the threshold for the acceptance of p-values

Output : $x_{d}$

Initialization : $n_{\rm s}=\lceil t_{\rm t}/t_{\rm s}\rceil$, the number of segments,

$\mathbf{S}=\{s_{1},...,s_{n_{\rm s}}\}$, the consecutive segments

1 $d=0$

2 while true do

3       $x_{d}=(1-B)^{d}x_{t}$ for $i=1$ to $n_{\rm s}$ do

4             Perform KPSS test on $s_{i}$

5             Compute the $p$-values for the level stationarity $p_{L}$-value and the trend stationarity $p_{T}$-value with $s_{i}$

6            

7      if $p_{L}\geq c_{p}$ and $p_{T}\geq c_{p}$ $\forall s_{i\in\{1,2,...,n_{s}\}}$ then

8             return $x_{d}$

9            

10      else

11             $d=d+1$

12            

13      

Algorithm 1 Integrated process (Sec. II C)

II.4 Autoregressive process

Once the stationary data set is obtained as an output $x_{d}$ of the integrated process, we build the AR model of $x_{d}$. A set of AR models can be produced by:

$$\hat{x}_{d,p}=\sum^{p}_{i=1}a_{i}x_{d,t-i}+\epsilon_{t}$$ (6)

where $p\in\{1,...,p_{\rm max}\}$.

For constructing a set of candidate models, we need to fix the upper-bound of $p$ which is set by the hyper-parameter $p_{\rm max}$. The optimal AR model is determined by model selection based on Akaike Information Criteria [30, AIC], which is defined as $AIC=2p+n\ln(\hat{\sigma_{p}}^{2})$ where $n$ is the sample size and $\hat{\sigma_{p}}^{2}$ is the maximum likelihood estimator for the variance of the noise term. For each model in the set, AIC is calculated. And the model that attains the optimal AIC will be selected.

We have examined the effect with different values of $p_{\rm max}$ by subtracting the corresponding optimal AR model from the data (cf. Equation 2). The results are shown in Figure 4 which show that $p_{\rm max}=8192$ (corresponding to $\sim$ 2 seconds for 4 kHz sampling frequency) can lead to recognizable waveform in the residuals. We have examined whether the result can be further improved by setting $p_{\rm max}$ at higher values. However, among all the experiments presented in this work, we found that the optimal $p$ selected by AIC all converged below $8192$ even with $p_{\rm max}$ set at higher values. In view of this, we fixed $p_{\rm max}=8192$ as the hyper-parameter throughout this work for an efficient computation.

In our framework, the model parameters (i.e. AR coefficients $\{a_{i}\}$) are estimated by Burg method, which fits the model to $x_{d}$ for minimizing the sums of squares of forward and backward linear prediction errors [31, 32]. The function ar.burg from the R package stats is adopted in our experiment.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: (upper panels) PSDs of the pure noise data, AR residuals and the whitened data. (lower panels) Q-Q plots for comparing the distributions of the raw/processed data (red lines) and with a Gaussian distribution (black lines).

Since AR is the major component of our procedure, before we apply it to the data with a CBC signal embeded, we have first investigated its performance on the pure noise data and examined whether the processed data can satisfy the requirements of stationarity and normality [cf. 33]. For this test, we have used both simulated and real noise data. We started by generating 100 simulated noise data of 32 s from sampling the updated Advanced LIGO sensitivity design curve . And we have processed them with algorithm 1 and algorithm 2. For comparison, we have also separately processed the simulated noise with the standard whitening. To quantify the difference between the distribution of the data from normality, we have run the Anderson-Darling (A-D) test [35]. Taking the $p-$value of 0.05 as the benchmark for rejecting the null hypothesis, all the simulated data fail to pass the A-D test which yield a mean $p-$value of $\sim 10^{-16}$. This suggests they are all significantly different from a Gaussian distribution. After subtracting the noise data from the AR models, $\sim 70\%$ of these samples become conform with normality (yield a $p-$value $>0.05$). And we found that whitening results in a similar fraction that pass A-D test.

On the other hand, all these simulated noise data are found to pass KPSS test and do not demonstrate any non-stationarity. In order to search for the non-stationary noise data for the experiment, we have searched over the LIGO data. We have chosen 24 time segments which do not encompass any confirmed GW events and yield a $p-$value smaller than 0.05 in the KPSS test. Also all these segments do not conform with normality which yield a mean $p-$value of $\sim 10^{-14}$ in the A-D test. After whitening (or processing with our method), all 24 processed pure noise data pass the KPSS test (all yield $p-$value $>0.1$). Also, both methods result in similar fraction ($\sim 65\%$) for passing the normality test. Therefore, we conclude that both our method and whitening have a comparable performance on the pure noise in attaining normality and stationarity. As an example, we compare the power spectral density (PSD) for one of our real noise sample with those of AR residuals and whitened data in Figure 3. These plots also demonstrate the capability of a AR model in line removal. In the low panels, we have also constructed the quantile-quantile (Q-Q) plots for comparing the distribution of the data with the normal distribution. They clearly show that both AR residuals and the whitened data distribute as a Gaussian.

For our test data which encompasses the transient signal from GW150914, an optimal AR order of $p=7931$ is obtained. Before passing to the next stage of processing, we obtain the residual time series $x_{r}$ by subtracting the optimal model from $x_{d}$ (cf. Equation 2, algorithm 2). The waveform of AR residual data is shown in the third panel in Figure 2, in which the modulation resulting from the BBH coalescence starts emerging. To further suppress the random fluctuation in $x_{r}$, we proceed to the next stage (see below).

Input : $x_{d}$, the input time series,

$p_{\rm max}=8192$: the maximum number of the order of the autoregressive process

Output : $x_{r}$

1 Estimate the parameters $\phi_{p}$ of an AR($p$) model on $x_{d}$ with the variance $\sigma^{2}_{p}$ calculated using Burg’s method and

2 $\{\phi_{p},p_{\rm opt}\}=\mathop{\arg\min}\limits_{\phi_{p},p_{\rm opt}}{AIC}$ for $p\in\{1,...,p_{\rm max}\}$

3 with $AIC=2p+n\ln(\hat{\sigma_{p}}^{2})$

4 Obtain time series $\hat{x}_{p_{\rm opt}}$ from the optimal model $\{\phi_{p},p_{\rm opt}\}$

5 return the residual $x_{r}=x_{d}-\hat{x}_{p_{\rm opt}}$

Algorithm 2 Autoregressive process (Sec. II D)

Figure 4: Effects of AR with different $p_{\rm max}$ on the difference data of $d=2$. These time series are the residuals obtained by subtracting the fitted AR model from the data (i.e. $x_{r}$ outputted by algorithm 2).

II.5 Moving-average process

Different from the conventional ARIMA model (Equation 5) in which MA performs the regression with the past forecast errors $\epsilon_{i}$. In our framework, MA refers to the method of estimating the trend in the residuals which is taken as a form of low-pass finite impulse response filter. The process is expressed as follows:

$$x_{\rm ma}=\frac{1}{q}\sum_{j=-k}^{k}x_{t+j}$$ (7)

where $q$ is the order of MA and $k=(q-1)/2$. Since we consider a two-sided (centered) MA, if an even value of order $q$ is specified, two MAs of $k$ rounded-down and rounded-up will be averaged [cf. 36].

For choosing the value of $q$, a model with small $q$ might have the signal remain buried by the random fluctuations in $x_{r}$. On the other hand, a large $q$ can smear out the signal. Since the GW from the BBH coalescence has its frequency varying, a MA model with a fixed $q$ can suffer from the aforementioned trade-off. Another problem of a single MA model is found when large values of $q$ are adopted. In the Figure 5, we show the power spectral density (PSD) of the output from the MA model with different $q$ annotated as blue lines in each panel. Within our concerned frequency band (32$-$512 Hz), We found that power spectral leakage starts appearing with $q\geq 7$, which can be possibly resulted from over-smoothing.

To overcome the aforementioned problems, rather than using only a single MA($q$), we adopted the method of Ensemble of Averages (EoA) [cf. 37] which combines MAs from a range of $q$ ($q\in\{$1…$q_{\rm max}\}$). It aggregates a number of MAs with an ensemble of moving averages. In our work, we utilize EoA and demonstrate that using median as the collector function can be a very effective filter (algorithm 3).

Figure 6 shows the outputs of EoA for different choices of $q_{\rm max}$. Empirically, we found that $q_{\rm max}=20$ with a median collector function gives a desirable result. It can eliminate most random fluctuations and retains the signal fidelity as all three stages of coalescence (i.e., inspiral, merger, and ringdown) can be clearly visualized. Furthermore, with EoA, the problem of spectral power leakage in PSD is resolved (annotated as red lines in Figure 5).

Input : $x_{r}$,

$q\in\{q_{\rm min},...,q_{\rm max}\}$, the order of the moving average

where $q_{\rm min}=1$ and $q_{\rm max}=20$ by default

Output : $x_{\rm EoA}$

1 for $q=q_{\rm min}$ to $q_{\rm max}$ do

2       $x_{q}=\frac{1}{q}\sum^{k}_{j=-k}x_{r+j}$

3$\hat{x}=median(\{x_{q_{\rm min}},...,x_{q_{\rm max}}\})$

4 return $\hat{x}$

Algorithm 3 Moving-average process (Sec. II E)

Figure 5: PSDs of the outputs from EoA with different $q_{\rm max}$ (red lines). These are compared with the corresponding result from MA with single value of $q$ (blue lines). Cases of different $q_{\rm max}$ and $q$ are shown.

Figure 6: EoA with a median collector on the output of AR stage. Collector function aggregates MAs from $q=1$ to $q=q_{\rm max}$.

II.6 Bandpass filtering

In the last step of our framework, we have $x_{\rm EoA}$ bandpass filtered in the frequency range of $32-512$ Hz for removing noise out of this band (e.g. seismic noise at low frequencies and photon shot noise at high frequencies). We adopt the Finite Impulse Response (FIR) filter by using the functions filtfilt and fir1 from the R package signal. The bandpass filtered signal of GW150914 is shown at the bottom panel in Figure 2.

In comparison with the output from the previous step (i.e. $x_{\rm EoA}$), no significant improvement can be found as a result of bandpass, which indicates that seqARIMA has already efficiently suppressed the noise in the raw data.

Although the bandpass does not appear to be necessary in the case of GW150914, we keep it in our framework to ensure all unwanted modulations outside this band are removed for the sake of comparing with the whitening results.

III Experimental results

III.1 Simulated data

For comparing the performance between the frequency-domain whitening filter and seqARIMA in noise reduction as well as waveform visualization, we have carried out a series of experiments. We started by simulating clean waveforms of a BBH coalescence at different luminosity distance $d_{L}$ by the code get_td_waveform from pycbc with the model SEOBNRv4_opt.

We have considered $d_{L}$ in a range from 200-4000 Mpc with a step size of $\Delta d_{L}$=200 Mpc. In order to analyse the denoising performance for a variety of waveform, for each $d_{L}$, we have generated 100 waveform of randomly sampled individual component masses $m_{1}$ and $m_{2}$. For the other parameters such as dimensionless spin and eccentricity, the default values of get_td_waveform are adopted (See ). These waveform are defined as the signals $s$.

For the sampling of waveform parameters, we have firstly fitted the distributions of $m_{1}$ and $m_{2}$ from all the 81 confirmed BBH CBC events with the R package gamlss. Among all the distribution functions available in gamlss, generalized Beta distributions of second kind provides the best description in accordance with AIC. And we sampled $m_{1}$ and $m_{2}$ from these best-fitted distributions.

For $d_{L}$, we have generated 100 noise data of 32 s, which is defined as $n$. They are sampled from a PSD simulated by aLIGOZeroDetHighPower in pycbc from LALsimulation with low_freq_cutoff of 15 Hz. Each of them are generated with different random seed. The preparation of the simulated data was finished by injecting $s$ in a random time location of $n$.

This simulated dataset allows us to compare the performance of seqARIMA and whitening filter in extracting the injected signal at varying $d_{L}$. In both methods of seqARIMA and whitening, the same bandpass filter of $32-512$ Hz were applied. For whitening, we have adopted a segment length of 4 s and the an overlap percentage of $50\%$ in all experiments. Such choices of whitening parameters follow the standards given by the pycbc documentation

For quantifying the fidelity of the extracted signal, we computed the cross-correlation functions (CCFs) defined as,

$${\rm CCF}(t^{\prime})=\sum_{t=-\infty}^{\infty}s(t)\hat{s}(t-t^{\prime})$$ (8)

where $s$ is the simulated waveform and $\hat{s}$ is the denoised data. Then we obtained the maximum values of $|{\rm CCF}|$, ${\rm CCF_{max}}$, as the metric of measuring the similarity between the $s$ and $\hat{s}$. In order to evaluate the noise reduction performance, we also computed the root-mean-square errors (RMSEs) defined as,

$${\rm RMSE}=\sqrt{\frac{\sum^{N}(s-\hat{s})^{2}}{N}}$$ (9)

where $N$ is the length of data, which reflects how the noise is suppressed in the whole time series. For each $d_{L}$, we have re-sampled $n$ with 100 different random seeds and computed the median and the 95% confidence interval of ${\rm CCF}_{\rm max}$ and RMSE from this sample.

In Figure 7, the results are shown for $d_{L}=$400, 2000, 3000, and 4000 Mpc. For a visual comparison of the similarity of the extracted signal and the injected waveform $s$, we have also overlaid $s$ in all the panels of Figure 7 as the red solid curves.

In the left panels of Figure 8, we compare how ${\rm CCF}_{\rm max}$ and RMSE vary with $d_{L}$ in both schemes. The error bars represent 95% confidence intervals calculated from 100 simulated waveform with randomly sampled $m_{1}$ and $m_{2}$ as well as different random seeds for generating the noise. Comparing the extracted signals by these two methods, we found that those obtained by seqARIMA generally have a larger degree of similarity with $s$ and lower level of noise. Although whitening process attains better results for small distance ($d_{L}\leq 600$ Mpc), seqARIMA has shown advantage in de-noising for increasing $d_{L}$ (i.e. larger ${\rm CCF}_{\rm max}$ and reduced RMSE).

In the right panels of Figure 8, we show the fractional improvements in both metrics as yielded by seqARIMA at different $d_{L}$. Comparing with the whitening results at $d_{L}=4000$ Mpc, seqARIMA has improved ${\rm CCF}_{\rm max}$ by $\sim 42~{}$% and suppressed RMSE by $\sim 23~{}$%.

Figure 7: The comparison of extracted waveforms from the simulated data (top panel) by whitening (middle panel) and seqARIMA denoising (bottom panel) for $d_{L}=$400, 2000, 3000, and 4000 Mpc, from top to bottom, respectively. In each panel, we have overlaid the injected signal (red lines) on the data (black lines). For the sake of comparison, the strain amplitudes are normalized in each panel.

Figure 8: (left panel) The comparison of ${\rm CCF_{max}}$ and RMSE between the injected signal and the waveform extracted from seqARIMA (red triangles) and the whitening process (blue circles) with varying $d_{L}$. The error bars represent 95% confidence errors calculated from 100 sampled parameters and 100 different random seeds. (right panel) Fractional difference of ${\rm CCF_{max}}$ and RMSE resulted from seqARIMA denoising with respect to the corresponding metrics resulted from whitening as a function of $d_{L}$.

III.2 LIGO data

To demonstrate the capability of seqARIMA in handling real data, we have attempted to extract the signals from a number of known GW events from the LIGO data. In this test, we have chosen all the events in GWTC-1 as observed during the first and second observation runs (O1 and O2) in 2015-2017 [44] plus two additional interesting events. All the data with a length of 4096 s with 4 kHz sampling frequency are obtained from GWOSC. Except for the NS-NS merger GW170817, we windowed the 4096 s data with a frame of 32 s for all the events. For GW170817, because of its much longer timescale, we apply a window of 50 s instead.

Event Name	Dataset	RMSE		CCF${}_{\rm max}$
		H1	L1	H1	L1
GW150914	Whiten	0.034	0.0288	0.628	0.488
	seqARIMA	0.024	0.0321	0.646	0.607
		-29.4 %	+11.4 %	+2.83 %	+24.4 %
GW151012	Whiten	0.0318	0.0326	0.127	0.0922
	seqARIMA	0.0304	0.0311	0.201	0.177
		-4.39 %	-4.71 %	+57.9 %	+92.1 %
GW151226	Whiten	0.0119	0.0119	0.0763	0.0736
	seqARIMA	0.0114	0.0117	0.145	0.109
		-4.05 %	-1.81 %	+90.6 %	+47.6 %
GW170104	Whiten	0.0289	0.0283	0.212	0.212
	seqARIMA	0.0262	0.0264	0.298	0.293
		-9.27 %	-6.71 %	+40.3 %	+38.2 %
GW170608	Whiten	0.0145	0.0144	0.0926	0.105
	seqARIMA	0.0138	0.0143	0.18	0.143
		-4.92 %	-1.23 %	+94 %	+37 %
GW170729	Whiten	0.0375	0.0391	0.278	0.379
	seqARIMA	0.033	0.0309	0.443	0.51
		-12.2 %	-20.9 %	+59.7 %	+34.5 %
GW170809	Whiten	0.0328	0.0309	0.246	0.305
	seqARIMA	0.0315	0.0304	0.367	0.35
		-3.75%	-1.5 %	+49.4 %	+15 %
GW170814	Whiten	0.0309	0.027	0.303	0.402
	seqARIMA	0.0255	0.0276	0.467	0.546
		-17.6 %	+2.0 %	+54 %	+35.8 %
GW170817	Whiten	0.00836	0.00828	0.0719	0.106
	seqARIMA	0.00816	0.00793	0.113	0.163
		-2.38 %	-4.16 %	+57.1 %	+54.1 %
GW170818	Whiten	0.035	0.0335	0.119	0.273
	seqARIMA	0.0325	0.0319	0.218	0.359
		-7.19 %	-4.99 %	+82.9 %	+31.2 %
GW170823	Whiten	0.0336	0.0349	0.306	0.3
	seqARIMA	0.0352	0.0337	0.434	0.362
		+4.78 %	-3.51 %	+41.7 %	+20.7 %
GW190814	Whiten	0.0118	0.0117	0.0951	0.109
	seqARIMA	0.0113	0.0112	0.171	0.182
		-4.68 %	-4.62 %	+79.7 %	+67.2 %
GW200105	Whiten	-	0.0121	-	0.0316
	seqARIMA	-	0.0115	-	0.174
		-	-5.59 %	-	+449 %

Table 1: The comparison of RMSE and CCF${}_{\rm max}$ yielded by whitening and seqARIMA on 13 confirmed CBC events as observed by LIGO with reference to the model waveform as specified in the literature. The third row of each event show the percentage change resulted from seqARIMA with respect to whitening. For GW200105, LIGO-H was not operational during this event as hence there is no data available [45].

III.2.1 GWTC-1 events

All 11 events in GWTC-1 can be well extracted by seqARIMA. In Figure 9, we show the spectrograms/oscillograms of the extracted signals from three representative cases, GW150914, GW151012, and GW170817, as detected by both observatories in Hanford (H: left panels) and Livingston (L: right panels). For the results of other GWTC-1 events, we have put them in the Appendix (Figure 11).

GW150914 is the first case that a GW signal was directly detected [4]. Its high signal-to-noise (SNR) of 26 has put it among the strongest signals of BBH merger detected so far. In Section II, we have already used this case for illustrating the feasibility of seqARIMA, in which we demonstrate that the signal of GW150914 can be clearly recovered. In the top row of Figure 9, we have produced the spectrograms of this event with Q-transform for visualizing how the frequency of the signal varies over the entire process. The characteristic sweeping chirp can be clearly seen in the spectrograms.

GW151012 is the BBH merger detected with a SNR of 10, which puts it as the weakest signal in GWTC-1 [44]. Its low significance as found from the initial discovery in O1 did not make it as a confirmed detection. And hence it was firstly considered as a candidate which was named as LVT151012 [46]. With a more detailed analysis, it was found to meet the criteria of a confident detection and was subsequently re-named as GW151012. In the second row of Figure 9, we show the spectrograms of the signals of GW151012 as extracted by seqARIMA. The chirp-like feature can be seen from the denoised data though it is not as clear as in the case of GW150914 because of its low significance.

The GW signal from the event GW170817 is resulted from a merging NS-NS binary, which is the first GW event that has the counterpart detected across the whole electromagnetic spectrum [5, 6, 7, 8]. It is associated with a short $\gamma-$ray burst GRB170817A, detected by Fermi Gamma-ray Burst Monitor (GBM) 1.7 s after the coalescence [8]. It has provided a long-sought evidence for the link between NS-NS mergers and short $\gamma-$ray bursts. Unlike BBHs, the inspiral time of GW170817 is much longer. Therefore, we take this event as a test for the capability of our framework in handling a signal with a longer timescale. Apart from adopting a wider window in the analysis, since LIGO-L data of GW170817 suffered from the transient noise (or glitch) at the GPS time of 1187008881.389 (around 1.1 s before the coalescence), we have used the data after noise subtraction following the glitch model described in [5]. The spectrograms of GW170817 resulted from seqARIMA denoising are shown in the bottom panels of Figure 9. The inspiral and the merging process over $\sim 30$ s can be clearly visualized in both data.

Figure 9: The spectrograms and oscillograms of seqARIMA-denoised LIGO data of GW150914, GW151012, and GW170817 as selected from GWTC-1 (Results for the other GWTC-1 events are shown in Figure 11). The color scale for the normalized power of the spectrogram is given at the upper-left corner of each panel. The reference epochs $t_{0}$ for each case are the event time reported in GWOSC. For GW170817, the glitch-removed LIGO-L data is used.

III.2.2 GW190814 & GW200105_162426

Apart from reproducing the GWTC-1 events, we have further tested our framework on two additional sources: GW190814 & GW200105_162426. These events were chosen because their inferred properties are somewhat different from those 11 events in GWTC-1.

GW190814 was detected in the third observing run (O3) with a SNR of 25 [48]. Parameter estimation suggests that the masses of the compact objects in their progenitor binary are highly unequal. While one component has its mass estimated as $\sim 23M_{\odot}$ which is consistent with a stellar BH, the mass of the other one is likely lying in a range of $\sim 2.5-3M_{\odot}$ which put it in a mass gap of being either a very massive NS or a low-mass BH. In the Fig. 1 of [48], we notice that the timescale of GW190814 is $\sim 2-3$ s long which is different from those of GWTC-1 sources. Therefore, we have included GW190814 in our test.

For the same reason, we have also included GW200105_162426 (hereafter GW200105) in our experiment. It was detected by a single detector (LIGO-L) during O3 with an SNR of $\sim 14$ [45]. It is estimated to have component masses of $\sim 8.9M_{\odot}$ and $\sim 1.9M_{\odot}$ which makes it likely an NS-BH binary. The signal of GW200105 shows a track of excess power with increasing frequency over $\sim 3$ s in the spectrogram [see Fig. 1 in 45].

In Figure 10, we show the spectrograms of these two sources produced in our framework. The tracks of the signals in both cases are clearly visible. In comparing the spectrogram of GW200105 resulted from seqARIMA and the one obtained from spectral whitening as shown in Fig. 1 of [45], we found that our result can attain a higher clarity which shows the inspiraling stage has a duration up to $\sim 6$ s.

In order to compare the performance of signal extraction by whitening and seqARIMA, we computed the ${\rm CCF}_{\rm max}$ and RMSE resulted from both schemes with reference to the waveforms generated by pycbc with the model of SEOBNRv4_opt for BBHs and IMRPhenomPv2 for GW170817 (BNS), GW190814 (Mass-gap), and GW200105 (NSBH) according to the parameters given in the corresponding literature. The results are summarized in Table 1. For comparing RMSE between seqARIMA and whitening, we have seen general improvement in most cases. However, there are a few cases that the noise reduction resulted from seqARIMA are worse than that from whitening. The most notable one is from the LIGO-L data of GW150914. This might suggest that for the events with SNR sufficiently large as in the case of GW150914, seqARIMA may not have the advantage over the conventional whitening. This is also reflected by the non-monotonic behavior for the small values of $d_{L}$ in Figure 8. On the other hand, in terms of ${\rm CCF}_{\rm max}$ (i.e. the similarity between the extracted signals and the model), seqARIMA has shown improvement in all our tested cases.

Figure 10: The spectrograms and oscillograms of seqARIMA-denoised data of GW190814 (LIGO-H, L) and GW200105_162426 (LIGO-L).

IV Summary & Future Prospects

In this work, we have proposed a novel de-noising technique in processing GW data, which is based on autoregressive modeling. By coupling with other techniques (i.e. integrated process, EoA), we have developed a framework we refer as seqARIMA pipeline (cf. Figure 1). The effects of each component in the pipeline have been investigated (see Figure 2 and subsection II.3-subsection II.6 for details). We have tested the performance of our proposed framework with a series of experiments.

We have examined the ability of seqARIMA pipeline in extracting the simulated GW signal with varying waveform and distance. By comparing the noise-subtracted time series and the injected signal (Figure 7), we have computed CCF${}_{\rm max}$ and RMSE resulted from both seqARIMA and whitening process. At larger distance, we found that seqARIMA can attain a higher CCF${}_{\rm max}$ and lower RMSE than those resulted from whitening (Figure 8).

We have also applied our method in extracting a number of known GW events from the LIGO data. All 11 events cataloged in GWTC-1 can be well recovered by seqARIMA (Figure 9 & Figure 11). We have further tested the method in two additional sources GW190814 (mass-gap object) and GW200105 (NS-BH merger), which have the timescale of their GW signals different from those in GWTC-1. We showed their signals can also be successfully extracted (Figure 10).

We have further compared the CCF${}_{\rm max}$ and RMSE resulted from both seqARIMA and whitening by comparing the noise-subtracted time series of these events with the model waveforms generated in accordance with the parameters specified in the corresponding literature (see Table 1). We found that seqARIMA generally yields improvement over whitening in terms of these performance metrics.

We have demonstrated that seqARIMA can enhance the noise suppression and therefore it is capable to provide an alternative to the conventional frequency-domain whitening process. For further improving the denoising performance, seqARIMA can be coupled with deep learning. Many recent studies have investigated the feasibility of denoising the GW data with deep neural network and showed that this can significantly suppress the noise and recover the signal [e.g. 49, 50, 51, 52]. We notice that these recent studies remain using whitening as a preprocessing procedure. Therefore, it will be encouraging to explore whether combining seqARIMA with these machine-learning based architectures can boost the denoising performance to a further extent. By substituting whitening with our proposed method, dedicated studies can also explore whether parameter estimation can also be benefited from seqARIMA.

We can also consider the feasibility of incorporating seqARIMA into a template-free low-latency detection pipeline. Since whitening can be a dominant source for the latency, it is desirable to reduce the computational cost in this stage [e.g. 19]. However, the conventional frequency-domain whitening process do not have many degree of freedom for improving the computational efficiency. On the other hand, the complexity of seqARIMA can be controlled by the hyper-parameters $p_{\rm max}$ and $q_{\rm max}$, which gives the flexibility of this process. For example, in trading off the fidelity of the extracted signal, a low $p_{\rm max}$ can result in a more efficient modeling. Therefore, one can examine whether seqARIMA can be adopted in a candidate identification pipeline. With the improved noise subtraction, the signal from a CBC or burst event can possibly be identified as a cluster of bright pixels in the spectrograms [e.g. 53, 54] which allows a GW event candidate to be detected without a priori knowledge of its waveform. A quantitative analysis on the execution speed of our proposed framework will be important for examining the capability of rapid real-time processing.

Appendix A Spectrograms and oscillograms of GWTC-1 events produced by our framework

Figure 11: The spectrograms and oscillograms of seqARIMA-denoised LIGO data of GWTC-1 sources which are not shown in Figure 9.

Figure 12: The spectrograms and oscillograms of seqARIMA-denoised LIGO data of GWTC-1 sources which are not shown in Figure 9.