SOAP: A generalised application of the Viterbi algorithm to searches for continuous gravitational-wave signals

TL;DR

This paper introduces a Viterbi algorithm-based search method for continuous gravitational waves that is model-agnostic, computationally efficient, and capable of detecting signals with unknown frequency evolution, demonstrated on real and simulated data.

Contribution

The authors develop a novel, computationally efficient Viterbi-based search method for continuous gravitational waves that is less sensitive but more versatile than traditional template-matching approaches.

Findings

Achieves 95% efficiency with 1% false alarm rate in various datasets.

Depth sensitivities of approximately 33, 10, and 13 Hz^{-1/2} in different data sets.

SNRs of about 60, 72, and 74 in the tested datasets.

Abstract

All-sky and wide parameter space searches for continuous gravitational waves are generally template-matching schemes which test a bank of signal waveforms against data from a gravitational wave detector. Such searches can offer optimal sensitivity for a given computing cost and signal model, but are highly-tuned to specific signal types and are computationally expensive, even for semi-coherent searches. We have developed a search method based on the well-known Viterbi algorithm which is model-agnostic and has a computational cost several orders of magnitude lower than template methods, with a modest reduction in sensitivity. In particular, this method can search for signals which have an unknown frequency evolution. We test the algorithm on three simulated and real data sets: gapless Gaussian noise, Gaussian noise with gaps and real data from the final run of initial LIGO (S6). We show that at 95% efficiency, with a 1% false alarm rate, the algorithm has a depth sensitivity of $\sim 33$, $10$ and $13$ ,Hz$^{-1/2}$ with corresponding SNRs of $\sim 60$, $72$ and $74$ in these datasets. we discuss the use of this algorithm for detecting a wide range of quasi-monochromatic gravitational wave signals and instrumental lines.

Tables

Table 1. Table 1 : Table shows the ranges of the search parameters for the gapless Gaussian injections search Sec.  III.1 and the optimised parameters. There are 10 parameter values spaced linearly between the limits.

	$\tau$	$w_S$	$w_L$	$p\left(M_L\right)/p\left(M_N\right)$
limits	[1.0,1.3]	[0.1,5.0]	None	0.0
optimised	1.1	2.06	None	0.0

Table 2. Table 2 : Table shows the ranges of the search parameters for the S6-like injections into Gaussian noise with gaps Sec.  III.3 and the optimised parameters. The parameters were spaced linearly between the limits.

	$\tau$	$w_S$	$w_L$	$p\left(M_L\right)/p\left(M_N\right)$
limits	[1.0,1.4]	[0.1,5.0]	None	0.0
optimised	1.1	2.06	None	0.0

Table 3. Table 3 : Table shows the ranges of the search parameters for the S6 search Sec.  III.3 and the optimised parameters. The parameters, τ 𝜏 \tau , w L subscript 𝑤 L w_{\rm{L}} and w S subscript 𝑤 S w_{\rm{S}} were distributed in log space between the limits and p ​ ( M L ) / p ​ ( M N ) 𝑝 subscript 𝑀 𝐿 𝑝 subscript 𝑀 𝑁 p(M_{L})/p(M_{N}) is distributed uniformly.

	$\tau$	$w_S$	$w_L$	$p\left(M_L\right)/p\left(M_N\right)$
limits	[1.0,1.1]	[0.1,5.0]	[0.1,6.0]	[0.0,1.0]
optimised	1.00000001	4.0	5.0	0.0387

Equations

$$p(\mathbf{\nu }\mid D)=\frac{p(\mathbf{\nu })p(D\mid \mathbf{\nu })}{p(D)},$$

$$p(\mathbf{\nu })$$

$$=p({\nu }_{N-1}\mid {\nu }_{N-2})p({\nu }_{N-2}\mid {\nu }_{N-3})\dots p({\nu }_1\mid {\nu }_0)p({\nu }_0)$$

$$=p({\nu }_0)\prod _{j=1}^{N-1}p({\nu }_j\mid {\nu }_{j-1}),$$

$$p(D\mid \mathbf{\nu })=\prod _{j=0}^{N-1}p({\mathbf{x}}_j\mid {\nu }_j),$$

$$p(\mathbf{\nu }|D)=\frac{p({\nu }_0)p({\mathbf{x}}_0|{\nu }_0)\prod _{j=1}^{N-1}p({\nu }_j|{\nu }_{j-1})p({\mathbf{x}}_j|{\nu }_j)}{\sum _S\left\{p({\nu }_0)p({\mathbf{x}}_0|{\nu }_0)\prod _{j=1}^{N-1}p({\nu }_j|{\nu }_{j-1})p({\mathbf{x}}_j|{\nu }_j)\right\}},$$

$$p(\hat{\mathbf{\nu }}|D)\propto \underset{\mathbf{\nu }}{\max }\left[p({\nu }_0)p({\mathbf{x}}_0|{\nu }_0)\prod _{j=1}^{N-1}p({\nu }_j|{\nu }_{j-1})p({\mathbf{x}}_j|{\nu }_j)\right].$$

$$\begin{split}\log p(\hat{\bm{\nu}}|D)=\max_{{\bm{\nu}}}{\biggl{\{}\log p(\nu_{0})+\log p({\bm{x}_{0}}|\nu_{0})}\\ \left.\sum_{j=1}^{N-1}\biggl{[}\log p(\nu_{j}|\nu_{j-1})+\log p({\bm{x}_{j}}|\nu_{j})\biggr{]}\right\}+\text{const}.\end{split}$$

$$p({\mathbf{x}}_j\mid {\nu }_{j,k},{\sigma }_{j,k},I)\propto \exp \left[C({\nu }_{j,k})\right].$$

$$C({\nu }_{j,k})\equiv C_{j,k}=\frac{1}{{\sigma }_{j,k}^2}\frac{1}{N_{\mathrm{s}}}{\left|\sum _{r=0}^{N_{\mathrm{s}}-1}{x}_{j,r}{\mathrm{e}}^{i{\nu }_{j,k}{t}_r}\right|}^2,$$

$$V_{j,k}=\underset{i}{\max }(C_{j,k}+T_i+V_{j-1,k+i}),$$

$$\begin{array}{r}p(\nu ,{\nu }^{(1)},{\nu }^{(2)}|D^{(1)},D^{(2)})\propto p(\nu )p({\nu }^{(1)},{\nu }^{(2)}|\nu )\\ p(D^{(1)}|{\nu }^{(1)})p(D^{(2)}|{\nu }^{(2)}),\end{array}$$

$$V_{j,k}=\underset{i,l,m}{\max }(C_{j,k+l}^{(1)}+C_{j,k+m}^{(2)}+T_{i,l,m}+V_{j-1,k+i}),$$

$$V_{j,k,s}=\underset{h}{\max }(C_{j,k}+T_{s,h}+V_{j-1,k+s,k+s+h}),$$

$$F_j=\sum _i^{N_s}{C}_{i,k}$$

$$\begin{array}{rl}p(D\mid \nu )& \propto p(x_1,x_2\dots x_n\mid \nu )\\ & \propto p(x_1\mid \nu )\dots p(x_n\mid \nu )\\ & \propto \exp \left(\sum _i{C}_{j,k}\right).\end{array}$$

$$p(F_j|{\nu }_j,M_{\text{N}},I)=\frac{1}{2^{d/2}\Gamma (d/2)}{F}_j^{d/2-1}\exp \left\{\frac{F_j}{2}\right\}$$

$$\begin{array}{r}p(F_j|{\nu }_j,\lambda ,M_{\text{S}},I)=\frac{1}{2}\exp \left\{-\frac{F_j+\lambda }{2}\right\}{\left(\frac{F_j}{\lambda }\right)}^{d/4-1/2}\\ I_{d/2-1}\left(\sqrt{\lambda F_j}\right).\end{array}$$

$$\begin{array}{r}p(F_j^{(1)}\mid {\nu }_j,M_{\mathrm{S}},I)={\int }_0^{\infty }p(\lambda ,w)p(F_j^{(1)}|{\nu }_j,\lambda ,M_{\text{S}},I)d\lambda .\end{array}$$

$$p(\lambda ,w)=\exp \left(\frac{-\lambda }{w}\right).$$

$$\begin{array}{r}p(M_{\mathrm{GL}}\mid F_j^{(1)},{\nu }_j,I)=p(M_{\mathrm{G}}\mid F_j^{(1)},{\nu }_j,I)\\ +p(M_{\mathrm{L}}\mid F_j^{(1)},{\nu }_j,I).\end{array}$$

$$\begin{array}{rl}{O}_{S/GL}^{(1)}(F_j^{(1)}\mid {\nu }_j)& =\frac{p(M_{\mathrm{S}}\mid F_j^{(1)},{\nu }_j)}{p(M_{\mathrm{GL}}\mid F_j^{(1)},{\nu }_j)}=\frac{p(M_{\mathrm{S}}\mid F_j^{(1)},{\nu }_j)}{p(M_{\mathrm{G}}\mid F_j^{(1)},{\nu }_j)+p(M_{\mathrm{L}}\mid F_j^{(1)},{\nu }_j)}\\ & =\frac{p(M_{\mathrm{S}})p(F_j^{(1)}\mid M_{\mathrm{S}},{\nu }_j)}{p(M_{\mathrm{G}})p(F_j^{(1)}\mid M_{\mathrm{G}},{\nu }_j)+p(M_{\mathrm{L}})p(F_j^{(1)}\mid M_{\mathrm{L}},{\nu }_j)}\\ & =\frac{p(F_j^{(1)}\mid M_{\mathrm{S}},{\nu }_j)p(M_{\mathrm{S}})/p(M_{\mathrm{G}})}{p(F_j^{(1)}\mid M_{\mathrm{G}},{\nu }_j)+p(F_j^{(1)}\mid M_{\mathrm{L}},{\nu }_j)p(M_{\mathrm{L}})/p(M_{\mathrm{G}})}\end{array}$$

$$\begin{array}{rl}\log \left[O_{S/GL}^{(1)}(F_j^{(1)})\right]& =\log \left[p(F_j^{(1)}\mid M_{\mathrm{S}})\right]\\ & -[\log (p(F_j^{(1)}\mid M_{\mathrm{G}})\\ & +p(F_j^{(1)}\mid M_{\mathrm{L}})p(M_{\mathrm{L}})/p(M_{\mathrm{G}}))]\end{array}$$

$$\begin{array}{r}{\hat{V}}_{i,j}=\underset{k,l,m}{\max }(T_{k,l,m}+b+V_{i-1,j+k}\\ +\log \left[O_{S/GL}^{(1)}\left(F_j^{(1)}\right)\right]),\end{array}$$

$$\begin{array}{r}p(F_j^{(1)},F_j^{(2)}\mid {\nu }_j,M_{\mathrm{S}},I)={\int }_0^{\infty }p(\lambda ,w_{{}_{\mathrm{S}}})\\ p(F_j^{(1)}|{\nu }_j,\lambda ,M_{\text{S}},I)p(F_j^{(2)}|{\nu }_j,\lambda ,M_{\text{S}},I)d\lambda .\end{array}$$

$$\begin{array}{r}p(F_j^{(1)},F_j^{(2)}\mid {\nu }_j,M_{\mathrm{L}},I)={\int }_0^{\infty }p(\lambda ,w_{{}_{\mathrm{L}}})\\ [p(F_j^{(1)}|{\nu }_j,M_{\mathrm{N}},I)p(F_j^{(2)}|{\nu }_j,\lambda ,M_{\mathrm{S}},I)\\ +p(F_j^{(1)}|{\nu }_j,\lambda ,M_{\mathrm{S}},I)p(F_j^{(2)}|{\nu }_j,M_{\mathrm{N}},I)]d\lambda ,\end{array}$$

$$\begin{array}{r}p(F_j^{(1)},F_j^{(2)}\mid {\nu }_j,\lambda ,M_{\mathrm{G}},I)=p(F_j^{(1)}\mid {\nu }_j,M_{\mathrm{G}},I)\\ p(F_j^{(2)}\mid {\nu }_j,M_{\mathrm{G}},I).\end{array}$$

$$\begin{array}{rl}\log \left[O_{S/GL}^{(2)}(F_j^{(1)},F_j^{(2)})\right]& =\log \left[p(F_j^{(1)},F_j^{(2)}\mid M_{\mathrm{S}})\right]\\ & -[\log (p(F_j^{(1)},F_j^{(2)}\mid M_{\mathrm{G}})\\ & +p(F_j^{(1)},F_j^{(2)}\mid M_{\mathrm{L}})p(M_{\mathrm{L}})/p(M_{\mathrm{G}}))]\end{array}$$

$$\begin{array}{r}{\hat{V}}_{i,j}=\underset{k,l,m}{\max }(T_{k,l,m}+b+V_{i-1,j+k}\\ +\log \left[O_{S/GL}^{(2)}\left(F_j^{(1)},F_j^{(2)}\right)\right]),\end{array}$$

$$\Delta f_{\mathrm{rot}}^{(1,2)}=\frac{({\mathbf{v}}^{(1)}-{\mathbf{v}}^{(2)})\cdot \hat{\mathbf{s}}}{c}{f}_0,$$

$$\Delta f_{\mathrm{rot}}\approx 6.5\times 10^{-7}{f}_0,$$

Full text

I Introduction

One of the main targets for current ground based gravitational-wave (GW) detectors, including Laser Interferometer Gravitational-wave Observatory (LIGO) Abbott et al. (2009); Aasi et al. (2015a) and Virgo Acernese et al. (2008, 2014), are sources of continuous gravitational waves. These are long-duration, quasi-monochromatic sinusoidal signals that are well-modelled by a Taylor series expansion in the signal phase. A likely source of such signals are rapidly spinning non axisymmetric neutron stars. A number of possible emission mechanisms are outlined in Prix (2009); Owen (2009).

These types of GWs are expected to give strain amplitudes that are significantly below the detector’s noise spectral density, and need sensitive search algorithms for detection. The most sensitive method is to use a coherent matched filter which requires knowledge of the waveform beforehand such that it can be coherently correlated with the data. This approach is used in searches for gravitational signals from known pulsars such as Dupuis and Woan (2005); Astone et al. (2010); Jaranowski et al. (1998a); Collaboration et al. (2017); Abbott et al. (2019a). For broad parameter space searches, where the parameters of the signal are unknown, a large number of template waveforms must be used to sufficiently cover the parameter space. This approach rapidly becomes computationally impractical as the search space grows, so semi-coherent search methods have been developed to deliver the maximum overall sensitivity for a given computational cost. Semi-coherent searches break the data up into sections of either time or frequency and perform a coherent analysis on these sections separately. These intermediate results can then be recombined incoherently in a number of different ways to form the final search result outlined in  Brady and Creighton (2000); Abbott et al. (2019b) and references therein.

The analysis that we present here is known as SOAP Ellis et al. (2006) and is based on the Viterbi algorithm Viterbi (1967). The algorithm models a process that has a discrete number of states at discrete time steps, and computes the set of states which gives the highest probability (suitably defined) given the data. Our implementation of SOAP is intended as a stand-alone search which is naturally non-parametric and has broad applications to both searches for known signal types and signals which have an unknown frequency evolution. The algorithm works in time-frequency plane, where our ‘states’ are represented by the time and frequency coordinates of a potential signal. We can then find the most probable set of frequencies a possible signal could have, i.e. we can find the most probable track in frequency as a function of time. This is not the first application of the Viterbi algorithm to GW data. Another variant of the algorithm Suvorova et al. (2016) has recently been used, amongst other applications, as part of a continuous wave (CW) search to track a pulsar with randomly wandering spin frequency Sun et al. (2018); Suvorova et al. (2017); Abbott et al. (2017, 2018); Sun and Melatos (2018).

In the next section we will describe the Viterbi algorithm and the basic SOAP implementation to GW time-frequency data. We then describe additional features to the algorithm, including the use of data from multiple detectors. As well as this we describe methods used to ignore instrumental effects in the data, such as incoherently summing data and a ‘line aware’ statistic. In the final section as well as a test of the computational cost of the search, we show results of a search performed on datasets of increasing complexity: Gaussian noise with no gaps (i.e., contiguous in time), Gaussian noise with gaps simulating real data more accurately, and finally real LIGO data taken during the sixth science run.

II The Viterbi Algorithm

The Viterbi algorithm is an efficient method for determining the most probable set of states (a single ‘track’ of steps on the time-frequency plane) in a Markov model dependent on data, where the model has a discrete number of states at each step. Rather than computing the probability of every possible track and selecting the most probable, the algorithm maximises this probability after every discrete step. As a result, a partial track which cannot ultimately be the most probable is rejected before the next step is calculated, and only a fraction of all possible tracks need to be computed to find the one that is most probable.

In this work we apply the Viterbi algorithm to a GW strain time-series to find the most probable track of a single variable-frequency signal in the noisy data. We divide the time series into $N$ equal-length and contiguous segments ${\bm{x}}_{j}$, defining the set $D\equiv\{{\bm{x}}_{j}\}$. The ‘states’ in the model correspond to the frequencies a signal could have in each segment. A ‘track’ is a list of such frequencies ${\bm{\nu}}\equiv\{\nu_{j}\}$, where $\nu_{j}$ is the frequency in the segment ${\bm{x}_{j}}$.

Our objective is to calculate the most probable track given the data, i.e., the track that maximises $p({\bm{\nu}}\mid D)$. Using Bayes theorem, this posterior probability can be written as

[TABLE]

where $p({\bm{\nu}})$ is the prior probability of the track, $p(D\mid{\bm{\nu}})$ is the likelihood of the track (i.e., the probability of the data given the track) and $p(D)$ is the model evidence (or marginalised likelihood).

The Viterbi algorithm treats the track as the result of a Markovian process, such that the current state depends only on the previous state. It is therefore useful to split the track’s prior into a set of transition probabilities such that

[TABLE]

where $p(\nu_{0})$ is the prior probability that the signal in the first time step has a frequency $\nu_{0}$ and $p(\nu_{j}\mid\nu_{j-1})$ is the prior ‘transition’ probability for $\nu_{j}$ given the frequency at the last step was $\nu_{j-1}$.

The noise in each of the segments can be treated as independent, so the likelihood component in Eq. 1 can be factorised as

[TABLE]

where $p({\bm{x}_{j}}\mid\nu_{j})$ is the likelihood of our signal having a frequency $\nu_{j}$ in the $j$th segment.

Using Eqs. 1,II and 3, the posterior probability is then

[TABLE]

where in the denominator we must sum over all possible tracks $S$. We require the specific track, or set of frequencies, $\hat{\bm{\nu}}$ that maximises the posterior probability, and which therefore maximises the numerator on the right-hand side of Eq. 4, i.e.,

[TABLE]

This track also maximises the log of the probability,

[TABLE]

The Viterbi algorithm finds the most probable track $\hat{\bm{\nu}}$ by calculating the quantities in Eq, 6 for each frequency at each time step. In the following sections we explain how this is achieved in practice.

II.1 The transition matrix

We define the ‘transition matrix’ $T$ as the matrix that stores the prior log-probabilities $\log p(\nu_{j}\mid\nu_{j-1})$. These transition probabilities depend only on the size and direction of the transition, and in our case correspond to a jump in frequency when moving from the $(j-1)$th to the $j$th state. It is within the transition matrix that we impose some loose model constraints. For example it is usual in the time-frequency plane for frequencies to only have discrete values (frequency bins) and a track might only be allowed to move by one bin in each time step, restricting it to a up, centre or down (UCD) transition or ‘jump’ or equivalently setting the size of the first dimension of the transition matrix $n_{1}=3$. We can also impose that the transition probabilities are independent of the current track location in frequency, i.e. $p(\nu_{j}\mid\nu_{j-1})=p(\nu_{j+k}\mid\nu_{j+k-1})$. This leads to the transition matrix containing only three numbers, corresponding to the three prior log-probabilities that the track was in the corresponding UCD frequency bin at the previous time step. These numbers are chosen to reflect the prior probability of a frequency deviation in the track and depend on the class of signals that one wishes to detect.

In later sections we will consider more complex situations in which the transition matrix describes the prior probability associated with sequences of even earlier transitions (‘memory’) and the case where there are multiple detectors. In these cases the number of dimensions of the transition matrix can grow substantially to account for the extra complexity of the problem.

II.2 Single detector

We will first consider the simple case of a single dataset $D$, generated by a single gravitational wave detector, and consider only a one-dimensional transition matrix, i.e. the transition matrix contains only three numbers refering the the track ‘jumps’. We will make use of discrete Fourier transforms so that frequencies, and hence the track frequencies, are also discrete. These frequencies will be indexed by $k$ and therefore $\nu_{j}\rightarrow\nu_{j,k}=k(j)\Delta f$ where $\Delta f=1/T$ is the frequency bin width for a segment of duration $T$.

The Viterbi algorithm determines the most probable track on the time-frequency plane by calculating the value of Eq. 6 for every discrete Fourier frequency, incrementally in time. In other words, at each time segment it finds the most probable earlier track which ends at each particular frequency. On reaching the final segment it can look back to identify the most probable track connecting segment 1 to segment $N$.

There are two main components to Eq. 6: the transition probabilities $p(\nu_{j}\mid\nu_{j-1})$ and the likelihoods $p({\bm{x}_{j}}\mid\nu_{j})$. The transition probabilities are pre-calculated and stored in a transition matrix according to Sec. II.1 above, and details of how these are calcuated are described in Sec. III. To calculate the likelihood we follow the approach of Bretthorst (1988) which gives, under the assumption of a single sinusoidal signal in additive Gaussian noise in data segment ${\bm{x}_{j}}$,

[TABLE]

where $C_{j,k}(\nu_{j,k})$ is the Schuster periodogram normalised to the noise variance at frequency $\nu_{j,k}$ of segment $j$. This is equivalent to the log-likelihood, and is defined as

[TABLE]

where $N_{\rm s}$ is the number of data points in each segment and $t_{r}$ is the time corresponding to $x_{j,r}$, the $r$th sample in the $j$th data segment. $\sigma_{j,k}^{2}$ is the noise variance and is calculated as an estimate of the noise power spectral density (PSD) in the $k$th sample and the $j$th data segment. The log-likelihoods of each segment can be calculated at discrete frequencies before running the algorithm by computing the power spectra for each segment from discrete Fourier transforms of the data. In the GW field these standard data forms are known as short Fourier transforms. It is worth noting at this point that it is also possible to write this as a likelihood ratio, and therefore write out detection statistic as a log-odds ratio, however, we will discuss this in more depth in Sec. II.6.

The Viterbi algorithm records two quantities for each frequency and time bin: The first, $V_{j,k}$, contains the value defined by Eq. 6, which is the log-probability of the most probable path ending in position $j,k$. The second, $A_{j,k}$, is the transition, or ‘jump’, used to achieve the most probable path. The algorithm can be divided into three main sections: initialisation, iteration and identification. These three sections are described in pseudo-code in Alg. 1 and a simple demonstration of the algorithm at work is shown in Fig. 1.

Initialisation

The two parts of Eq. 6, $\log p(\nu_{0})$ and $\log p({\bm{x}_{0}}\mid\nu_{0})$, must be computed before the main recursive part of the algorithm can start. Therefore, the initialisation section (lines 5–8) in Alg. 1 calculates the first column in the lower panel of Fig. 1. A priori, there is no preferred initial frequency, so we take the log-prior $\log p(\nu_{0,k})$ to be uniform over the complete frequency range. As a result, this is does not affect the maximisation for any jump, therefore, can be omitted from the calculation. We then use the pre-calculated log-likelihood values $C_{0,k}$ to fill the track probabilities $V_{0,k}$. There is no previous position to jump from in this case, so the transition probabilities are irrelevant and $A_{0,k}$ are set to zero.

Iteration

The main part of the calculation is the sum in Eq. 6. Lines 11–16 in Alg. 1 calculate the most probable tracks that end at each frequency bin for each segment by using

[TABLE]

where $i$ is the size and direction of the jump. For example, in Fig. 1 columns 1–4 are calculated in order using Eq. 9, where it maximises over three possible previous positions in frequency. These positions are the frequency bins UCD of the current position. The size and direction of the jump, $i$, which gives the maximum probability is then saved to $A_{j,k}$. These are shown in parentheses below the log-probabilities in Fig. 1 where UCD correspond to values of $i=[-1,0,1]$ respectively.

Identification

The final stage of the algorithm identifies the most probable track. This is done by initially finding the highest log-probability values in the final time segment, $\max_{k}(V_{N-1,k})$ (line 19 in Alg. 1). In the lower panel of Fig. 1 this is located at position $j,k=4,1$ with $V_{4,1}=22$. To find the track which corresponds to this, the values in $A_{jk}$ are followed backwards from this position (lines 20–21). For example, in Fig. 1 the final position is $j,k=4,1$ and $A_{j,k}=\rm{Center}=0$, this means that at the previous segment the most probable track was at position $j,k=4-1,1+0=3,1$. At this time $A_{3,1}=R=1$, therefore, the next track element is at $j,k=3-1,1+1=2,2$. This then continues until $j=0$ whereupon these retraced positions constitute the most probable track, highlighted in red in Fig. 1.

The most probable track is the one traced backwards from the highest probability final segment frequency position. However, tracks can also be traced back from any of the end-frequency positions, returning the most probable track conditional on a given final position. Such tracks should not be confused with the being equal to the second, third, fourth, etc. most probable tracks. Information regarding the rankings and properties of all possible tracks (excluding the most probable and conditionally most probable tracks) is lost during the maximisation procedures computed at each stage in the algorithm – a necessary consequence of the algorithm’s speed and efficiency.

II.3 Multiple detectors

If there are $Q$ detectors operating simultaneously we have $Q$ sets of data which can be combined appropriately to provide input to the Viterbi search described above. We must also modify the allowed transitions encoded within the transition matrix to take account of the extra prior constraints that are now available.

The received instantaneous frequency of a given astrophysical signal will be nearly the same for all ground-based GW detectors, and our algorithm should be sensitive to tracks that show this consistency in frequency. However there will be small differences between the frequencies measured at detectors that are not co-located, due to differential Doppler shifts caused by Earth rotation. As a result the signal could fall in different frequency bins at each detector.

To account for these small differences in signal tracks in each detector, we reference the observed tracks to a third (pseudo) detector located at the centre of the Earth which would be insensitive to Earth spin. The signal frequencies in each real detector are then allowed to vary within a certain number of frequency bins from the track in the reference detector. In the examples that follow, we only consider the possibilities that the track in each real detector is no more that one frequency bin away from the reference track. We can tune the length of the SFTs to ensure this is a valid assumption, this is explained in greater detail in III. As well as differences in signal frequency, due to antenna patterns and other effects, the measured signal amplitude may differ between the detectors Jaranowski et al. (1998b). In the following example we assume that the signal has the same amplitude in each detector, however, in Sec. II.6 we discuss the case where they differ.

We will now show how the algorithm in Sec. II.2 can be modified to handle a two-detector network (i.e., $Q=2$), however any number of detectors can easily be accommodated. In the two detector case the joint probability of two (real) tracks, $\nu^{(1)}$ and $\nu^{(2)}$, and the geocentric track $\nu$, given the data, is

[TABLE]

where $D^{(1)}$ and $D^{(2)}$ represent the data from the two detectors. The main difference between this and that described in Sec. II.2 is that the track probabilities $V_{j,k}$ are stored for the geocentric pseudo-detector. The main iterative calculation (defined for the single detector case in Eq. 9) now becomes

[TABLE]

where ${C}^{(1)}$ and ${C}^{(2)}$ refer to the log-likelihoods in detectors 1 and 2 respectively and the transition matrix $T$ is an $n_{1}\times n_{2}\times n_{3}$ matrix, where $n_{1}$ dimension refers to the jump from the previous time step, $n_{2}$ and $n_{3}$ refer to the relative frequency positions in each real detector. The transition matrix is now three-dimensional and holds the prior log-probabilities of $p(\nu)$ and $p(\nu^{(1)},\nu^{(2)}|\nu)$. We now need to maximise over three indices: $i,l$ and $m$. The index $i$ refers to the size and direction of the jump at the geocentre (as before). The indices $l$ and $m$ refer to the number of frequency bins by which the two real tracks deviate from the geocentre track. For example, if the most probable track in the geocentred detector is in bin $j,k=5,12$ and the values of $i,l,m=0,-1,1$, then detector 1 is in position $j,k={5,11}$ and detector 2 is in position $j,k={5,13}$ and the geocentred track was in the position $j,k={4,12}$ at the previous time step. As a result, the track at the geocentre is only affected by Doppler modulations from the Earth’s orbit whereas the tracks in the real detectors include Doppler modulations from the Earth’s spin.

At every time step the frequency bin position for each real detector is forced to be within $n_{l}$ or $n_{m}$ bins of the track in the geocentred detector, where $n_{l}$ and $n_{m}$ depend on how much each detector could possibly be Doppler shifted. As mentioned previously, we only consider the case where $n_{l}=1$ and $n_{m}=1$, allowing the track from each real detector to be at most one frequency bin away from the geocentred track position. While we tune the SFT length to keep this condition for different frequencies, it is also possible to tune the values of $n_{l}$ and $n_{m}$ to get a similar effect. The implementation of the multi-detector algorithm is similar to the single detector case described in Sec. II.2. However in the single detector case there is only a single variable to be maximised over for each time-frequency bin. This variable is the frequency jump from the position in the previous segment. For the multi-detector case there are at least three variables to be maximised over: the probability of the jump, $i$, at the geo-centre and the probability of the signal being in the surrounding positions in each on $Q$ real detectors, $l,m,\dots$. The values of $i,l,m,\dots$ are then saved to $A_{j,k}$ and are ultimately used to reconstruct the most probable consistent tracks in each real detector.

As in Sec. II.2, there are three main sections: Initialisation, iteration, and the identification. For the multi-detector case each element is modified as follows.

Initialisation

The first-row calculation (lines 5–8) in Alg. 1, are now modified to additionally maximise over the real detector track positions $l$ and $m$. For each time-frequency bin the maximum sum of the log-likelihoods is saved together with the frequency locations of the corresponding tracks in the real detectors. The index $i=0$ is kept constant as there is no previous position.

Iteration

To process the subsequent time segments, lines 13–14 in Alg. 1 are modified to account for two (or more) detectors. Line 13 of Alg. 1 is changed to calculate Eq. 11, the log-probability of a track at the geocentre ending in bin $j,k$ given that signal is in the real detector positions of $j,k+l$ and $j,k+m$. Line 14 is then modified so that $A_{j,k}$ stores the jump values, $i$, and the real detector positions, $l$ and $m$, which returned the highest probability.

Identification

The most probable track is identified in the same way as for the single detector case, first by finding the maximum value in the final time step of $V_{j,k}$ (line 19 in Alg. 1). The track at the geocentre can then be found by iteratively following the jump values stored in $A_{j,k}$ back from this position. The track in each of the real detectors is determined by using the values of $l$ and $m$ indices also stored in $A_{j,k}$ to find the relative position of the track in each real detector compared to the geocentre.

This method can be extended to more than two detectors by including additional datasets and expanding the corresponding number dimensions of the maximisation procedures in the iterative steps.

II.4 Memory

In this section we extend the basic Viterbi algorithm to improve its sensitivity to non-stochastic signals where there is some knowledge of its frequency evolution. We do this by including a form of ‘memory’ and this extension applies to both the single and multiple-detector cases. Rather than considering only the previous step in our decision-making process, we now include the previous $m+1$ steps and expand the transition matrix to include these values. A memory of $m=0$ therefore corresponds to the methods described in previous sections. With a non-zero memory the transition matrix can a-priori make certain sequences of jumps more probable and assign different prior probabilities for these jump sequences e.g., ‘up then centre’ may be less preferable to ‘centre then centre’. As a result we can increase the chance of the most probable track matching an expected astrophysical signal. In a single detector search with a memory of $m=1$, if we only allow UCD transitions, then for every frequency bin we save 3 values. These are proportional to the log-probabilities of a track coming from a UCD bin in the previous time step, where the maximisation is over the corresponding UCD bins two time steps back. Eq. 11 then is then modified to,

[TABLE]

where $s$ and $h$ refer to the UCD jumps at the time step $j-1$ and $j-2$ respectively. Similar to the previous two sections, the algorithm is split into three parts: initialisation, iteration, and the track identification:

Initialisation

The initialisation process needs to populate the first $m+1$ steps before the main iteration can start. At the first time step, the elements $V_{0,k,s}$ are set to the log-likelihoods $C_{0,k}$ as in Sec. II.2. There is no previous time step, so the element $s$ is not relevant. At the second time step, $V_{1,k,s}$ is calculated using Eq. 12, where there is no maximisation over $h$, it is assumed to be [math], or a center jump. As there is no data before $j=0$, the maximisation at this point will always return the jump which has the largest prior probability, which in this case is a center jump. Therefore, the maximisation returns the same value for all frequency bins and can be set to a center jump.

Iteration

For all following time steps the values for each element of $V_{j,k,s}$ in Eq. 12 are calculated. This quantity is proportional to the log-probability of the track ending in time-frequency bin $j,k$, which was in the previous position of $j-1,k+s$. The corresponding value of $h$ that maximised the log-probability of the track is recorded in $A_{j,k,s}$.

Identification

The most probable track is identified in a similar way to the non-memory cases, by finding the highest-valued last element, $V_{N-1,k,s}$. The values of $s$ and $h$ are then followed back to find the most probable track. As an example, let us assume the most probable track finishes in bin $j,k,s=10,5,0$, where the value of $m$ is $A_{10,5,0}=1=\rm{up}$. The previous position is then $j,k,s=10-1,5+s,m=10-1,5+0,1=9,5,1$ with a value $A_{9,5,1}=0=\text{Center}$, and the next track position is $j,k,s=9-1,5+1,0=8,6,0$ etc. The values of $j,k$ along this track describes most probable path.

The number of elements over which one must search increases rapidly with memory length, and has a strong impact on the computational cost of the analysis. For the single detector Viterbi approach the number of calculations made is $3\times N\times M$ if we only allow UCD jumps, where $N$ and $M$ are the number of time are frequency bins respectively. When memory is included this increases to $3^{m+1}\times N\times M$.

II.5 Summed input data

In this section we outline a method of incoherently-summing a set of SFTs to increase the signal-to-noise-ratio (SNR) of a signal in a segment. To be more precise, we actually sum the log-likelihoods, i.e. the quantity in Eq. 8. We can write the new summed set of data $F_{j}$ as,

[TABLE]

where $N_{s}$ is the number of SFTs to sum together and the log-likelihood $C(\nu_{i,k})$ is defined in Eq. 8. We can see this is possible by looking at Eq. 7, where we can use the product of likelihoods,

[TABLE]

If the data contains gaps where the detector was not observing, then we fill the gaps in the power spectrum with a constant value which is the expectation value of the log-likelihood. The procedure of filling in the gaps of the data is completed before any summing. Therefore, the data should have the same mean regardless of how much real data is in each sum. In the examples that follow, we sum the SFTs over the length of one day.

The main motivation for summing the data is to increase the SNR of a signal in the segments. The risk is that a signal can move between adjacent frequency bins during a day. To reduce this risk, we choose the frequency bin width such that it is more likely that a signal will be contained within a single frequency bin than cross a bin edge. In practise, to ensure that this is true, the segment or SFT length and the number of segments which are summed can be tuned for each search. As well as increasing the SNR, summing over one day should average out the antenna pattern. This means that the log-likelihood value in any bin should be more similar between detectors, however, there is still some variation due to the sky localisation and polarisation.

This also has two main effects on the transition matrix, the first is that as each segment of data is now one day long, a jump between frequency bins is far more likely, therefore, the transition matrix elements are modified to account for this. The second is that as the data is averaged over one day, the signal should remain is the same frequency bin between detectors, therefore, there is no longer a need for the multi-dimensional transition matrix described in Sec. II.3.

The volume of the data is also reduced by a factor of $1/N_{s}$, therefore, the time taken for the algorithm to run is also reduced by the same factor.

II.6 Line-aware statistic

The single-detector algorithm described in Sec. II.2 returns the most probable track of the loudest signal assumed to be in Gaussian noise. However, an astrophysical signal is not expected to have an amplitude which is orders of magnitude above the noise floor, but have an amplitude more similar to the noise. Therefore, a signal with a large amplitude is more likely to be of instrumental origin rather than astrophysical Coughlin et al. (2010); Aasi et al. (2015b); Covas et al. (2018).

We first consider the model of Gaussian noise with no signal present. Within a single summed segment, the likelihood of Gaussian noise at frequency $\nu$ is given by a $\chi^{2}$ distribution,

[TABLE]

where $F_{j}$ is the frequency domain power summed over sub-segments within a single day, as described in Sec. II.5 and $d$ is the number of degrees of freedom, equal to twice the total number of summed SFTs. $M_{\rm{N}}$ represents the model that the data is simply Gaussian noise. In the presence of a signal (model $M_{\text{S}}$), the power should follow a non central $\chi^{2}$ distribution in which the non-centrality parameter $\lambda$ is the square of the SNR, $(\lambda=\rho_{\rm{opt}}^{2})$, i.e.

[TABLE]

If a signal is present we therefore expect the SFT powers in the detector to follow Eq. 16. We can then determine the evidence for model $M_{\text{S}}$ by marginalising over $\lambda$,

[TABLE]

Here we set the prior on $\lambda$ to be an exponential distribution of width $w$, this is done somewhat arbitrarily as we expect the majority of signals to have a low SNR. This distribution follows,

[TABLE]

In this single-detector case, we expect an astrophysical signal to look very similar to that of a line other than its amplitude (or SNR). Therefore, we set the evidence for an astrophysical signal and an instrumental signal to follow Eq. 17, where the width $w$ different between the two models.

We then have three models, one for an astrophysical signal, one for an instrumental line and one for Gaussian noise.

The posterior probability of model $M_{\text{GL}}$, which contains the probability of Gaussian noise or Gaussian noise with a line (taken as mutually exclusive) is

[TABLE]

We can now find the posterior odds ratio for the presence of a signal over noise or a line,

[TABLE]

In practice it is convenient to use the log odds ratio,

[TABLE]

As we are only interested in the maximum of $\log\left[O^{(1)}_{\rm{S/GL}}(F^{(1)}_{j})\right]$, the factor $\log\left[p(M_{\rm{S}})/p(M_{\rm{G}})\right]$ can be dropped from the expression.

In this version of the Viterbi algorithm, rather than storing a value proportional to the log-probabilities as in Sec. II.3, here we store a value proportional to the log-odds ratio. Here we take the log-odds ratio defined in Eq. 21 and add the log-prior odds $p(\bm{\nu}\mid M_{S})/(p(\bm{\nu}\mid M_{N})+p(\bm{\nu}\mid M_{L}))$ which is the log-prior or any particular track. By assuming that the track transitions for the line and noise model are equally probable for any jump, we set the denominator of the prior-odds is a constant $b$. This then means Eq. 9 is modified to,

[TABLE]

where $\hat{V}$ refers to a log-odds ratio. The maximised statistic now has three tuneable parameters: the width ($w_{S}$) in Eq. 18, on the prior for a signal SNR squared, $p_{\rm S}(\lambda)$, the width ($w_{L}$) of the prior in the case of a line, $p_{\rm L}(\lambda)$, and the ratio of the prior on the line and noise models, $p({M_{\rm L}})/p({M_{\rm G}})$. These parameters are optimised for each search, where we initially estimate the SNR of a signal we hope to be sensitive to in each time slice, then use this as a guide for the width of the signal prior. This is then repeated for an expected line SNR and this is used for the width of the line prior. The ratio of line and noise models runs in the range 0 to 1, we set this limit as we do not expect an instrumental line to be as likely as Gaussian noise in any particular frequency bin.

This line-aware statistic can be applied in a more powerful way when we use multiple detectors and is similar to the approach in Keitel et al. (2014). The multiple-detector algorithm described in Sec. II.3 returns the most probable track of a common signal assumed to be in Gaussian noise. As a consequence the algorithm will return large values of the log-likelihood even if there are inconsistent values of SFT power between the detectors, either from non-Gaussian noise or because the signal is not equally strong in the two detectors. However a signal with unequal power in the two detectors is more likely to be a non-Gaussian instrumental line than an astrophysical signal. The line-aware statistic described in this section is designed to make the search more robust to such instrumental artefacts within realistic non-Gaussian data whilst maintaining sensitivity to astrophysical signals.

For most of the analysis examples presented here we use data which is the incoherent sum of 30-minute normalised SFTs over a day (described in more detail in Sec. II.5). As a result the effects of the detector antenna patterns and of differential Doppler shifts are significantly reduced, and any signal should have a broadly similar summed log-likelihood in the same frequency bin in each detector. The statistic can then be modified such that we expect a similar log-likelihood in each detector.

In a similar way to the single-detector case, we can write out the evidence for each of the three models as follows. If a signal is present we therefore expect the SFT powers in both detectors to follow Eq. 16. Assuming for the moment that the noise variance is the same in both, we can determine the evidence for model $M_{\text{S}}$ by marginalising over $\lambda$,

[TABLE]

We set the prior on $\lambda$ the same as in the single detector case in Eq. 18. In this case, if an instrumental line is present in one of the detectors we expect to see signal-like power in that detector and noise-like power in the other. The evidence for this ‘line’ model ($M_{\text{L}}$) is therefore

[TABLE]

The third option is the simple case of approximately Gaussian noise in both of the detectors,

[TABLE]

We can now find the posterior odds ratio for the presence of a signal over noise or a line by following the same steps as in Eq. 20. Once again we write this as a log-odds ratio,

[TABLE]

The factor $\log\left[p(M_{\rm{S}})/p(M_{\rm{G}})\right]$ can again be dropped from the expression.

For the multi-detector case we then modify Eq. 11 to,

[TABLE]

where $\hat{V}$ refers to a log-odds ratio. This is then optimised over the same three parameters as the single detector case.

III Testing the algorithm

The sensitivity of the algorithm was tested by searching for artificial signals from isolated pulsars added to three types of noise-like data: continuous Gaussian noise, Gaussian noise but with periods of missing data, and real detector data (the S6 mock data challenge (MDC) Walsh et al. (2016)). The S6 MDC refers to a standardised set of simulated signals which are injected into real data, this set is also what is used for the injections into the two Gaussian noise cases. We describe each of the tests in more detail in Sec. III.1,III.2 and III.3, but several common pre-processing steps are performed before running these datasets through the Viterbi algorithm:

We read SFTs generated from 1800 s stretches of data in 2 Hz bands between 100 and 200 Hz. The SFTs length is chosen to ensure that any signal is likely to be contained within the width of a single frequency bin during the length of one day, rather than being split across the bin edges (see below). 2. 2.

We estimate the noise PSD for each SFT by calculating a running median over frequency using LALSuite code XLALSFTtoRnmed LIGO Scientific Collaboration (2018), this includes a bias factor to convert this to the mean and has a width of 100 bins. We then normalise the SFT by dividing it by its running median, giving the noise-like parts of the spectrum a mean power of approximately one. 3. 3.

The SFTs are then summed over one day, as described in Sec. II.5. The signal parameters are chosen so that within the frequencies of the search, the signal will not fall in more than two frequency bins over this period.

The differential Doppler shift of a signal seen at two detector sites due to the Earth’s rotation $\Delta f^{(1,2)}_{\rm{rot}}$ is simply

[TABLE]

where ${\bm{v}}^{(1,2)}$ is the velocity of detector $1,2$ in an inertial reference frame, $f_{0}$ is the instantaneous signal frequency in the frame, $\hat{\bm{s}}$ is the unit vector in the direction of the source and $c$ is the speed of light. The maximum difference in frequency seen by the two LIGO detectors is

[TABLE]

so the frequency measured from a source in the equatorial plane with $f_{0}=200$ Hz will differ by up to $1.3\times 10^{-4}$ Hz in the two detectors. This is $\sim 4$ times smaller than the frequency bin width of 1800 s SFTs ($5.6\times 10^{-4}$ Hz), so signals at frequencies lower than this are likely to appear in the same frequency bin in the two detectors. Therefore, whilst at higher frequencies we still allow the signal to be in different frequency bins between the detectors, in the following searches, we do not allow this. 4. 4.

The data is then split into 0.1 Hz sub-bands which are overlapping by 0.05 Hz. These were chosen to ensure that signals are contained within a sub-band over the year. On these timescales the important contributions to the frequency evolution are the spin-down rate of the pulsar and the Doppler shift due to the earth orbit. To investigate the doppler shift, we can look at a signal at 200 Hz, using Eq. 28 we can calculate the maximum shift in frequency due to the earths orbit as,

[TABLE]

where $T_{\text{o}}$ and $R_{\text{o}}$ are the earth orbit time and radius. This gives a maximum doppler shift of $0.019\;\rm{Hz}$, this is a $\sim 1/5$ of the width of a sub-band, therefore, is more likely to be totally contained within a sub-band than crossing over the edge. To account for the cases where the signal frequency crosses over the edge of a sub-band, the sub-bands overlap by 0.05 Hz so that the majority of the signals should be completely contained within at least one of the sub-bands. To investigate the spin-down of the pulsar, we look at the length of data, $T=4.05\times 10^{7}$ s and we choose a sub-band width of 0.1 Hz. For a signal to drift over the width of a whole sub-band we would need f-dot of,

[TABLE]

The majority of the injections that follow satisfy this condition, signals which are greater than this, and therefore drift over multiple bands, are vetoed from the search. 5. 5.

The two detector Viterbi algorithm is then run using the line aware statistic (see Sec. II.6). There are 4 parameters which we optimise in this search. The transition probabilities, where we have one parameter $\tau$ which is the ratio of the probability of going straight to the probability of going either up or down. Due to the averaging procedure, the signals received at each detector are forced to follow a common track which is equal to the ‘imaginary’ detectors track. The other three parameters, $w_{S},w_{L}\;\rm{and}\;p(M_{L})/p(M_{N})$, are described in Sec. II.6. 6. 6.

The algorithm then returns the most probable track though the data, and the value $\propto$ the log-odds in the final time step, i.e., the maximum final value, $\max_{j}(V_{N,j})$, in Eq. 22, which is then our detection statistic.

As an example of what the algorithm returns, Fig. 2 shows the tracks in the two detectors, H1 and L1. This also shows the log-odds ratio of ending in any frequency bin, i.e., all the elements in Eq. 22. In this figure, each time segment of the odds ratios have been normalised such that the sum of the odds ratios is 1.

In the following tests there are two main quantities which we use to determine the sensitivity. These are sensitivity depth $\mathcal{D}$ and the optimal SNR $\rho$. The sensitivity depth, $\mathcal{D}$, is defined in Behnke et al. (2015) as,

[TABLE]

where $S_{h}(f)$ is the single-sided noise PSD and $h_{0}$ is the GW amplitude. The optimal SNR is defined as,

[TABLE]

where $X$ indexes the detectors and $\tilde{h}(f)$ is the Fourier transform of the time series of the signal $h(t)$. This expression is defined in Prix (2007) for a double-sided PSD and we have defined it for the more common single-sided case.

III.1 S6 injections into gapless Gaussian noise

The first test involves injecting signals into Gaussian noise. The power spectrum of a Gaussian noise time-series follows a $\chi^{2}$ distribution with two degrees of freedom, therefore, as we search through the power spectrum, we generate spectrograms which follow a $\chi^{2}$ distribution. These spectrograms are 0.1 Hz wide and are set at 0.05 Hz intervals between 100 Hz and 200 Hz. The bins are $1./1800$ Hz wide and 1800s long, where the total length of data is the same as S6, i.e. $\sim 1.3$ years. We then generate the signals, where the pulsars parameters are fixed to the same values as the injections in the S6 MDC in this band, these values are outlined in Walsh et al. (2016).

The values of $f_{0}$ for the injections were not always centred in a sub-band, therefore a number of sub-bands contained only part of the injected signal. These sub-bands were ignored as they contaminated the signal statistics and only the sub-band which contained the whole signal was accepted. This reduced the number of sub-bands from 2000 to 1762 with the removal of 238 sub-bands containing only part of a signal. This set also includes signals that drift across multiple sub-bands due to their high spin-down rate. Only two signals were removed due to their spin-down values, which were $>5\times 10^{-9}$ Hz/s, these were the two hardware injections in the 100-200 Hz band.

For each injection the SOAP algorithm returns the detection statistic described in Sec. II.6 and III. We calculate a false alarm rate, which is the fraction of bands that have no injection that do exceed a given threshold. This is set to 1% and is used as a detection threshold. We then take all of the bands and if they pass the threshold we set them as detected, i.e. 1, and if they do not they are set as not detected, i.e. 0. This then leaves us with a set of binomial data, where the efficiency curves later in the paper are sigmoids which have been fitted to this. The sigmoid follows,

[TABLE]

The fit is done by sampling the posterior, i.e.,

[TABLE]

where $p(x_{0},k)$ is the prior and we set to a flat prior and $p(x\mid x_{0},k)$ is the likelihood function which is defined by,

[TABLE]

To plot the efficiency curves and lower and upper error bounds, we sample Eq. 35 using Markov-Chain Monte Carlo (MCMC) and then take the mean and the 5th and 95th percentiles respectively for each point in SNR or depth and plot these. Fig. 3(a) and 3(c) then shows the efficiency curves for the analyses plotted against the signals optimal SNR and depth respectively. The parameters of the search and their optimised values are shown in Tab. 1. Where we set the prior on the line model to 0 as this part is irrelevant to this search due to the lack of lines in the data.

From this we can determine that in Gaussian noise without gaps, the Viterbi algorithm can detect to an SNR of $\sim 60$ and a depth of $\sim 33$ Hz*-1/2* with 95% efficiency at a 1% false alarm.

Fig. 3(b) and 3(d), show the RMS of the difference between the injected signal track and the track found by Viterbi for SNR and sensitivity depth respectively. This shows that at SNR of 60, where we are detecting signals with a 95% efficiency, the signals have a mean RMS of $\sim 2$ frequency bins.

III.2 S6 injections into Gaussian noise with gaps

In the second test, we attempt to more closely mirror the S6 MDC Walsh et al. (2016) in two stages. The first uses the same injection method as Sec. III.1 however, removes the SFTs where there are gaps in S6. The second uses the same injection method again including gaps, however, uses a different value for the noise floor for each SFT, this is calculated for each band and SFT from S6 data.

Both detectors in S6 had a duty cycle of $\sim$50% Aasi et al. (2015c), which means that there are sections of time where there is no data in either one or both detectors. In the sections where one detector is observing but the other is not, the multi detector statistic will not behave correctly as it only has access to data from a single detector. In these sections we switch from using the multi-detector statistic to the single-detector statistic using the same parameters, these are both defined in defined in Sec. II.6.

The process of removing sub-bands and generating efficiency curves is the same as in Sec. III.1.

We set a 1% false alarm rate and generate an efficiency curve for SNR and depth in Fig. 3(a) and Fig. 3(c) respectively. From these efficiency plots we can see to an SNR of $\sim 72$ or a depth of $\sim 13$ Hz*-1/2* at a 95% confidence with a false alarm of 1%.

The parameters of the search which were optimised and their optimised values are shown in Tab. 2.

In Fig. 3(b) and  3(d) show the RMS of the difference between the injected signal track and the track found by Viterbi for SNR and sensitivity depth respectively. This shows that at SNR of 72, where we are detecting signals with a 95% efficiency, the signals have a mean RMS of $\sim 10$ frequency bins.

III.3 Tests on the S6 MDC

For a more direct comparison to other CW searches and to see how the algorithm performs with real data, we test the two detector SOAP algorithm using the S6 MDC. We focus this search on the 100-200 Hz band, there are two main reasons for this, one being that this is LIGOs most sensitive band and the other is that for much higher frequencies the signal will drift over larger frequency ranges, therefore, our SFT length will have to be changed. Here the 1800 s SFTs are split as in Sec. III, whereafter normalisation, the data is split into 0.1 Hz wide sub-bands overlapping by 0.05 Hz.

The two detector SOAP algorithm using the line-aware statistic in Sec. II.6 is then run on each sub-band under the assumption that the detectors have the same sensitivity. For this search we have four parameters which we optimise, the ranges and optimised values are shown in Tab. 3.

As in Sec. III.2, only the sub-bands which contained the entire frequency evolution of the signal were selected. Out of the 2000 sub-bands, 238 were removed due the sub-band only containing part of the signals frequency evolution. The main difference between the analysis for Gaussian noise and real data is that the real data is contaminated with instrumental lines. This means that whilst the techniques described in Sec. II.6 reduce the number of contaminated bands with a high statistic value, there are still instrumental lines which are coincident between the detectors and which could not be removed with these techniques. Within the data there are large number of lines at integer Hertz, which are seen in coincidence between the two detectors, these are thought to originate from digital electronics Coughlin et al. (2010). Therefore the frequency bins $\pm 1$ bin of each integer frequency in Hertz were removed and filled with the expectation value of the noise. To remove instrumental effects at other frequencies, the sub-bands which gave values of our statistic above a chosen threshold were investigated by eye. In this case 344 sub-bands were investigated, and any which were contaminated were vetoed. From these 344 sub-bands, 193 were removed from the analysis. The predominant feature in the bands which were removed were broad spectral features which lasted the whole run. Therefore, out of the 2000 sub-bands which are searched over, a total number of 431 sub-bands were removed.

The process to calculate the efficiency curves is the same as in Sec. III.2 and III.1.

Fig. 3(c) and Fig. 3(a) show the efficiency curves for SNR and depth respectively. These show that we can detect and SNR of $\sim 74$ and a sensitivity depth of $\sim 13$ Hz*-1/2* with an efficiency of 95% at a false alarm of 1%. These results can then be compared to other searches in the S6 MDC comparison paper Walsh et al. (2016). Whilst we only search in the 100 - 200 Hz range, the closest comparison in Walsh et al. (2016) is the test in the 40 - 500 Hz range, such as in Fig. 4 in Walsh et al. (2016). Here our algorithm sits roughly in the middle of all other searches in terms of sensitivity.

III.4 Computational cost

One of the main strengths of this search is the drastically reduced computational cost when compared to other current CW searches. The scaling of the computing cost can be estimated for a single detector by looking at the number of calculations that need to be made. The number of calculations for a single detector search, $N^{(1)}_{\rm{calcs}}$ is,

[TABLE]

where $n_{1}$ is the size of the transition matrix, $N$ is the number of SFTs, $M$ is the number of frequency bins and $m$ is the amount of memory described in Sec. II.4. Where the computing cost scales linearly with the number of frequency bins and SFTs. In the following test we ignore ‘memory’ and look at the time taken for the single detector search where the time taken to read and save data is ignored. Here the data is the same size as the S6 MDC for a single detector search and the search is over a 0.1 Hz band, where we set $n_{1}=3$. This test, and the following test, was run locally on a MacBook Air with a 1.3 GHz Intel Core i5 processor .We can then write the time taken ,$T$, as,

[TABLE]

where $N_{\rm{bands}}$ is the number of different frequency bands. For the multiple, $Q$, detector case, we can then generalise Eq. 37 and write the number of calculations $N^{(Q)}_{\rm{calcs}}$ as,

[TABLE]

where $n_{1}$ is the first dimension of the transition matrix, $Q$ is the number of detectors and $n_{q+1}$ is the size of the transition matrix element which refers to detector $q$. For our tests we set $n_{1}=n_{q+1}=3$ and use 2 detectors i.e. $Q=2$ which each have the same size data as the previous test. The actual time taken to run however, depends on the version of the algorithm which is run. For example, including the line aware statistic slows the search slightly. For the two detector where two SFT powers are summed,

[TABLE]

The same search now including the line aware statistic, which is implemented using a lookup table, changes this to,

[TABLE]

IV Discussion

In this paper we describe an application of the Viterbi algorithm, called SOAP, to search for continuous sources of gravitational waves. This paper outlines the method and derives the statistics behind the method in a consistent Bayesian formalism. It then presents the results from the first set of tests of sensitivity for the SOAP algorithm on three separate datasets.

We tested SOAP on a set of fake isolated pulsar signals in the 100 – 200 Hz range, based on 1800s SFTs summed over one day. The three datasets that inclded these signals comprised continuous Gaussian noise, Gaussian noise but with temporal gaps corresponding to LIGO dead times in the S6 data run, and real data, i.e. the S6 MDC. Although a major attraction of SOAP is its sensitivity to a wide range of signal types, in the tests above it was optimised to detect isolated pulsar signals below 100 Hz with low spin-down to offer a comparison with other CW searches. From these tests, by setting a 95% efficiency and a false alarm of 1%, we found that in the case of continuous Gaussian data we could detect a signal with an optimal SNR of $\sim 60$ and a depth of $\sim 33$ Hz*-1/2* with an RMS of the difference between the injected and Viterbi track being $\sim 2$ frequency bins. When gaps were introduced into the data to simulate S6 we could detect a signal with an SNR $\sim 72$ and a depth of $\sim 10$ Hz*-1/2*, with an RMS of $\sim 10$ bins. The drop in sensitivity here is simply because there is $\sim 50\%$ less data compared to the previous case. Finally, in the S6 MDC we could detect a signal with an SNR $\sim 74$ and a depth of $\sim 13$ Hz*-1/2*. These real data contain non-Gaussian artefacts such as instrumental lines and this causes a further drop in sensitivity. Whilst not a full comparison to other searches in the S6 MDC Walsh et al. (2016), as we only tested on a subset of the bands, this search has a sensitivity which is comparable to some other CW searches, however offers a massive increase in speed.

We chose the specific frequency band to search over as the data which we used, i.e. the summed data, becomes less effective at frequencies much higher than 200 Hz, and using the parameters of our simulations, signals can spread over many frequency bins in a day, reducing sensitivity further, however this can be mitigated by using shorter SFTs or performing their summation over 12 (rather than 24) hours.

The methods described in this paper present a basic approach for gravitational-wave signal searches using SOAP. However there are several further developments that could increase its sensitivity. Some of these are outlined below:

One of the main features which reduces the sensitivity of the search is non-gaussianities within the data, namely instrumental lines. Although we have a statistic which penalises these features, in some cases it will also penalise a strong signal. For example, when the amplitude of the noise floor is high for one detector or the duty factor is lower, the signal will appear more like an instrumental line to this statistic. We hope to improve the search statistic in the future by searching for consistent amplitudes as opposed to consistent SNR, i.e the statistic will take the amplitude of the noise floor and the duty factor into account.

One variation of this method which has been described in this paper is ‘memory’, which is where the tracks jump in frequency is determined by the previous $n$ jumps. This has yet to be fully tested, however, we expect that this will increase our sensitivity to signals where have a better idea of their frequency evolution. This however, comes at a cost in computational time which we can estimate given Eq. 39 in Sec. III.4.

Further additions to the search include using the Fourier transform of the SFT power along the Viterbi track as a detection statistic. If the Viterbi track follows that from an astrophysical signal, then we should see the effects of the antenna pattern in this Fourier transform as a peak at half a sidereal day. If the track follows something which is not astrophysical then this should not be seen this peak in this Fourier transform. This only applies to the search directly on the SFTs not the summed data, as the antenna pattern variations will have been averaged out in the summing.

As well as searching for astrophysical signals, SOAP can also be used to search for and identify instrumental lines. Here we use single detector data, or multiple channels from a single detector, to identify quasi-monochromatic features on the data for further study.

Whilst this paper presents initial tests on sensitivity, further tests will be needed for a full comparison to other CW search methods. This search, however, aims to look for signals which may not follow the standard frequency evolution and is intended to return potentially interesting candidates for a more sensitive followup.

V Acknowledgements

We would like to acknowledge Dr Matt Pitkin for his help during this project, as well as the continuous wave working group of LIGO-Virgo Collaboration. This research is supported by the Science and Technology Facilities Council., G.W. and C.M. are supported by the Science and Technology Research Council (grant No. ST/L000946/1) . The authors are grateful for computational resources provided by the LIGO Laboratory and Cardiff University and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459 and STFC grant ST/I006285/1.