Consequently one can detect such a signal by correlating the data with the waveform of the signal and maximizing the correlation with respect to the parameters of the waveform. Using this method one can pick up a weak signal from the noise by building a large signal-to-noise ratio over a wide bandwidth of the detector [79] . This observation has led to a rapid development of the theory of gravitational-wave data analysis. It became clear that the detectability of sources is determined by optimal signal-to-noise ratio, Equation ( 22 ), which is the power spectrum of the signal divided by the power spectrum of the noise integrated over the bandwidth of the detector.

An important landmark was a workshop entitled Gravitational Wave Data Analysis held in Dyffryn House and Gardens, St. Nicholas near Cardiff, in July 1987 [74] . The meeting acquainted physicists interested in analyzing gravitational-wave data with the basics of the statistical theory of signal detection and its application to detection of gravitational-wave sources. As a result of subsequent studies the Fisher information matrix was introduced to the theory of the analysis of gravitational-wave data [33, 49] . The diagonal elements of the Fisher matrix give lower bounds on the variances of the estimators of the parameters of the signal and can be used to assess the quality of astrophysical information that can be obtained from detections of gravitational-wave signals [27, 48, 15] . It was also realized that application of matched-filtering to some sources, notably to continuous sources originating from neutron stars, will require extraordinary large computing resources. This gave a further stimulus to the development of optimal and efficient algorithms and data analysis methods [75] .

A very important development was the work by Cutler et al. [26] where it was realized that for the case of coalescing binaries matched filtering was sensitive to very small post-Newtonian effects of the waveform. Thus these effects can be detected. This leads to a much better verification of Einstein's theory of relativity and provides a wealth of astrophysical information that would make a laser interferometric gravitational-wave detector a true astronomical observatory complementary to those utilizing the electromagnetic spectrum. As further developments of the theory methods were introduced to calculate the quality of suboptimal filters [7] , to calculate the number of templates to do a search using matched-filtering [66] , to determine the accuracy of templates required [20] , and to calculate the false alarm probability and thresholds [42] . An important point is the reduction of the number of parameters that one needs to search for in order to detect a signal. Namely estimators of a certain type of parameters, called extrinsic parameters, can be found in a closed analytic form and consequently eliminated from the search. Thus a computationally intensive search needs only be performed over a reduced set of intrinsic parameters [49, 42, 51] .

Techniques reviewed in this paper have been used in the data analysis of prototypes of gravitational-wave detectors [65, 63, 6] and in the data analysis of presently working gravitational-wave detectors [77, 12, 3, 2, 1] .

We use units such that the velocity of light $c=1$  .

2 Response of a Detector to a Gravitational-Wave Signal

There are two main methods to detect gravitational waves which have been implemented in the currently working instruments. One method is to measure changes induced by gravitational waves on the distances between freely moving test masses using coherent trains of electromagnetic waves.

The other method is to measure the deformation of large masses at their resonance frequencies induced by gravitational waves. The first idea is realized in laser interferometric detectors and Doppler tracking experiments [39, 56] whereas the second idea is implemented in resonant mass detectors [10] .

Let us consider the response to a plane gravitational wave of a freely falling configuration of masses. It is enough to consider a configuration of three masses shown in Figure  1 to obtain the response for all currently working and planned detectors. Two masses model a Doppler tracking experiment where one mass is the Earth and the other one is a distant spacecraft. Three masses model a ground-based laser interferometer where the masses are suspended from seismically isolated supports or a space-borne interferometer where the three masses are shielded in satellites driven by drag-free control systems.

Figure 1 : Schematic configuration of freely falling masses as a detector of gravitational waves. The masses are labelled 1, 2, and 3. The optical and radio paths are denoted by $L_i$  , where the index $i$  corresponds to the opposite mass. The unit vectors ${\hat{\mathbf{n}}}_i$  point between pairs of masses, with the orientation indicated.

Let ${\nu }_0$  be the frequency of the coherent beam used in the detector (laser light in the case of an interferometer and radio waves in the case of Doppler tracking). Let us assume for simplicity that the distance between the masses is the same and equal to $L$  . Let ${\hat{\mathbf{n}}}_i$  , $i=1,2,3$  , be the unit vectors along the lines joining the test masses, let $O$  be the point lying in the plane of the three masses and equidistant from the masses, and let ${\mathbf{p}}_i$  , $i=1,2,3$  , be the vectors of length $l=L/\sqrt{3}$  joining $O$  and the masses. Let $y_{21}$  be the relative change $\Delta \nu /{\nu }_0$  of frequency induced by a transverse, traceless, plane gravitational wave on the coherent beam travelling from the mass $2$  to the mass $1$  , and let $y_{31}$  be a relative change of frequency induced on the beam travelling from the mass 3 to the mass 1. The equations for $y_{21}$  and $y_{31}$  are given by (see [31, 8] and also [71] for a coordinate free derivation)

$$\begin{array}{ccc}{y}_{21}\left(t\right)& =& \left(1-\hat{\mathbf{k}}\cdot {\hat{\mathbf{n}}}_3\right)\left({\Psi }_3(t+\hat{\mathbf{k}}\cdot {\mathbf{p}}_2-L)-{\Psi }_3(t+\hat{\mathbf{k}}\cdot {\mathbf{p}}_1)\right),\end{array}$$

(1)

$$\begin{array}{ccc}{y}_{31}\left(t\right)& =& \left(1+\hat{\mathbf{k}}\cdot {\hat{\mathbf{n}}}_2\right)\left({\Psi }_2(t+\hat{\mathbf{k}}\cdot {\mathbf{p}}_3-L)-{\Psi }_2(t+\hat{\mathbf{k}}\cdot {\mathbf{p}}_1)\right),\end{array}$$

(2)

where

$$\begin{array}{c}{\Psi }_j\left(t\right):=\frac{{\Phi }_j\left(t\right)}{1-(\hat{\mathbf{k}}\cdot {\hat{\mathbf{n}}}_j{)}^2},{\Phi }_j\left(t\right):=\frac{1}{2}{\hat{\mathbf{n}}}_j^T\cdot \tilde{\mathsf{H}}\left(t\right)\cdot {\hat{\mathbf{n}}}_j,j=1,2,3.\end{array}$$

(3)

In Equation ( 3 ) $\tilde{\mathsf{H}}$  is the three-dimensional matrix of the spatial metric perturbation produced by the wave in the proper reference frame of the detector, $\hat{\mathbf{k}}$  is the unit vector in this reference frame directed from the center $O$  to the source of the gravitational wave and T denotes matrix transposition.

In the source frame the three-dimensional matrix $\mathsf{H}$  of the spatial metric perturbation produced by the gravitational wave is given by

$$\begin{array}{c}\mathsf{H}\left(t\right)=\left(\begin{array}{ccc}{h}_{+}\left(t\right)& h_{\times }\left(t\right)& 0\\ h_{\times }\left(t\right)& -h_{+}\left(t\right)& 0\\ 0& 0& 0\end{array}\right),\end{array}$$

(4)

where $h_{+}$  and $h_{\times }$  are the two polarizations of the wave. In general the detector is moving with respect to the source of a gravitational wave and the signal registered by the detector will be modulated by this motion. To obtain the matrix $\tilde{\mathsf{H}}$  it is convenient to introduce a coordinate system with respect to which the source is fixed. The origin of this coordinate system is usually chosen to be the solar system barycenter (SSB). Then we obtain $\tilde{\mathsf{H}}$  by the transformation

$$\begin{array}{c}\tilde{\mathsf{H}}\left(t\right)={\mathsf{O}}_2^{-1}\left(t\right){\mathsf{O}}_1\left(t\right)\mathsf{H}\left(t\right){\mathsf{O}}_1^{-1}\left(t\right){\mathsf{O}}_2\left(t\right).\end{array}$$

(5)

Here ${\mathsf{O}}_1$  is the transformation matrix from the source frame to SSB and ${\mathsf{O}}_2$  is the transformation matrix from the detector frame to SSB (see [35, 19, 42, 51] for details).

The difference of the phase fluctuations $\Delta \phi \left(t\right)$  measured, say, by a photo detector, is related to the corresponding relative frequency fluctuations $\Delta \nu$  by

$$\begin{array}{c}\frac{\Delta \nu }{{\nu }_0}=\frac{1}{2\pi {\nu }_0}\frac{d\Delta \phi \left(t\right)}{dt}.\end{array}$$

(6)

For a standard Michelson, equal-arm interferometric configuration $\Delta \nu$  is given in terms of one-way frequency changes $y_{21}$  and $y_{31}$  [see Equations ( 1 ) and ( 2 )] by the expression [80] 

$$\begin{array}{c}\frac{\Delta \nu }{{\nu }_0}=\left(y_{31}\left(t\right)+y_{13}(t-L)\right)-\left(y_{21}\left(t\right)+y_{12}(t-L)\right).\end{array}$$

(7)

In the long-wavelength approximation Equation ( 7 ) reduces to

$$\begin{array}{c}\frac{\Delta \nu }{{\nu }_0}=2L\left(\frac{1}{2}{\hat{\mathbf{n}}}_i^T\cdot \dot{\tilde{\mathsf{H}}}\left(t\right)\cdot {\hat{\mathbf{n}}}_i-\frac{1}{2}{\mathbf{n}}_j^T\cdot \dot{\tilde{\mathsf{H}}}\left(t\right)\cdot {\mathbf{n}}_j\right).\end{array}$$

(8)

Consequently the phase change is given by

$$\begin{array}{c}\Delta \phi \left(t\right)=4\pi {\nu }_{0}Lh\left(t\right),\end{array}$$

(9)

where the function

$$\begin{array}{c}h\left(t\right)=\frac{1}{2}{\hat{\mathbf{n}}}_i^T\cdot \tilde{\mathsf{H}}\left(t\right)\cdot {\hat{\mathbf{n}}}_i-\frac{1}{2}{\hat{\mathbf{n}}}_j^T\cdot \tilde{\mathsf{H}}\left(t\right)\cdot {\hat{\mathbf{n}}}_j\end{array}$$

(10)

is the response of the interferometer to a gravitational-wave signal in the long-wavelength approximation. In this approximation the response of a laser interferometer is usually derived from the equation of geodesic deviation (where the response is defined as the relative change of the length of the two arms, i.e., $h\left(t\right):=\Delta L\left(t\right)/L$  ). There are important cases where the long-wavelength approximation is not valid. These include the space-borne LISA detector for gravitational-wave frequencies larger than a few mHz and satellite Doppler tracking measurements.

In the case of a bar detector the long-wavelength approximation is very accurate and the detector's response $h_B\left(t\right)=\Delta L\left(t\right)/L$  , where $\Delta L$  is the change of the length $L$  of the bar, is given by

$$\begin{array}{c}{h}_B\left(t\right)={\hat{\mathbf{n}}}^T\cdot \tilde{\mathsf{H}}\left(t\right)\cdot \hat{\mathbf{n}},\end{array}$$

(11)

where $\hat{\mathbf{n}}$  is the unit vector along the symmetry axis of the bar.

In most cases of interest the response of the detector to a gravitational-wave signal can be written as a linear combination of four constant amplitudes $a^{\left(k\right)}$  ,

$$\begin{array}{c}h(t;a^{\left(k\right)},{\xi }^{\mu })={\sum }_{k=1}^4{a}^{\left(k\right)}{h}^{\left(k\right)}(t;{\xi }^{\mu })={\mathbf{a}}^T\cdot \mathbf{h}(t;{\xi }^{\mu }),\end{array}$$

(12)

where the four functions $h^{\left(k\right)}$  depend on a set of parameters ${\xi }^{\mu }$  but are independent of the parameters $a^{\left(k\right)}$  . The parameters $a^{\left(k\right)}$  are called extrinsic parameters whereas the parameters ${\xi }^{\mu }$  are called intrinsic. In the long-wavelength approximation the functions $h^{\left(k\right)}$  are given by

$$\begin{array}{c}\begin{array}{ccc}{h}^{\left(1\right)}(t;{\xi }^{\mu })& =& u(t;{\xi }^{\mu })cos\phi (t;{\xi }^{\mu }),\\ h^{\left(2\right)}(t;{\xi }^{\mu })& =& v(t;{\xi }^{\mu })cos\phi (t;{\xi }^{\mu }),\\ h^{\left(3\right)}(t;{\xi }^{\mu })& =& u(t;{\xi }^{\mu })sin\phi (t;{\xi }^{\mu }),\\ h^{\left(4\right)}(t;{\xi }^{\mu })& =& v(t;{\xi }^{\mu })sin\phi (t;{\xi }^{\mu }),\end{array}\end{array}$$

(13)

where $\phi (t;{\xi }^{\mu })$  is the phase modulation of the signal and $u(t;{\xi }^{\mu })$  , $v(t;{\xi }^{\mu })$  are slowly varying amplitude modulations.

Equation ( 12 ) is a model of the response of the space-based detector LISA to gravitational waves from a binary system [51] , whereas Equation ( 13 ) is a model of the response of a ground-based detector to a continuous source of gravitational waves like a rotating neutron star [42] . The gravitational-wave signal from spinning neutron stars may consist of several components of the form ( 12 ). For short observation times over which the amplitude modulation functions are nearly constant, the response can be approximated by

$$\begin{array}{c}h(t;A_0,{\phi }_0,{\xi }^{\mu })=A_{0}g(t;{\xi }^{\mu })cos\left(\phi (t;{\xi }^{\mu })-{\phi }_0\right),\end{array}$$

(14)

where $A_0$  and ${\phi }_0$  are constant amplitude and initial phase, respectively, and $g(t;{\xi }^{\mu })$  is a slowly varying function of time. Equation ( 14 ) is a good model for a response of a detector to the gravitational-wave signal from a coalescing binary system [79, 18] . We would like to stress that not all deterministic gravitational-wave signals may be cast into the general form ( 12 ).

3 Statistical Theory of Signal Detection

The gravitational-wave signal will be buried in the noise of the detector and the data from the detector will be a random process. Consequently the problem of extracting the signal from the noise is a statistical one. The basic idea behind the signal detection is that the presence of the signal changes the statistical characteristics of the data $x$  , in particular its probability distribution.

When the signal is absent the data have probability density function (pdf ) $p_0\left(x\right)$  , and when the signal is present the pdf is $p_1\left(x\right)$  .

A full exposition of the statistical theory of signal detection that is outlined here can be found in the monographs [87, 47, 83, 82, 57, 37, 68] . A general introduction to stochastic processes is given in [85] .

Advanced treatment of the subject can be found in [55, 86] .

The problem of detecting the signal in noise can be posed as a statistical hypothesis testing problem. The null hypothesis $H_0$  is that the signal is absent from the data and the alternative hypothesis $H_1$  is that the signal is present. A hypothesis test (or decision rule) $\delta$  is a partition of the observation set into two sets, $\mathcal{ℛ}$  and its complement ${\mathcal{ℛ}}^{\prime }$  . If data are in $\mathcal{ℛ}$  we accept the null hypothesis, otherwise we reject it. There are two kinds of errors that we can make. A type I error is choosing hypothesis $H_1$  when $H_0$  is true and a type II error is choosing $H_0$  when $H_1$  is true. In signal detection theory the probability of a type I error is called the false alarm probability, whereas the probability of a type II error is called the false dismissal probability. $1-\text{(false dismissal probability)}$  is the probability of detection of the signal. In hypothesis testing the probability of a type I error is called the significance of the test, whereas $1-\text{(probability of type II error)}$  is called the power of the test.

The problem is to find a test that is in some way optimal. There are several approaches to find such a test. The subject is covered in detail in many books on statistics, for example see references [45, 34, 53] .

3.1 Bayesian approach

In the Bayesian approach we assign costs to our decisions; in particular we introduce positive numbers $C_{ij}$  , $i,j=0,1$  , where $C_{ij}$  is the cost incurred by choosing hypothesis $H_i$  when hypothesis $H_j$  is true. We define the conditional risk $R$  of a decision rule $\delta$  for each hypothesis as

$$\begin{array}{c}{R}_j\left(\delta \right)=C_{0j}{P}_j(\mathcal{ℛ})+C_{1j}{P}_j\left({\mathcal{ℛ}}^{\prime }\right),j=0,1,\end{array}$$

(15)

where $P_j$  is the probability distribution of the data when hypothesis $H_j$  is true. Next we assign probabilities ${\pi }_0$  and ${\pi }_1=1-{\pi }_0$  to the occurrences of hypothesis $H_0$  and $H_1$  , respectively. These probabilities are called a priori probabilities or priors. We define the Bayes risk as the overall average cost incurred by the decision rule $\delta$  :

$$\begin{array}{c}r\left(\delta \right)={\pi }_0{R}_0\left(\delta \right)+{\pi }_1{R}_1\left(\delta \right).\end{array}$$

(16)

Finally we define the Bayes rule as the rule that minimizes the Bayes risk $r\left(\delta \right)$  .

3.2 Minimax approach

Very often in practice we do not have the control over or access to the mechanism generating the state of nature and we are not able to assign priors to various hypotheses. In such a case one criterion is to seek a decision rule that minimizes, over all $\delta$  , the maximum of the conditional risks, $R_0\left(\delta \right)$  and $R_1\left(\delta \right)$  . A decision rule that fulfills that criterion is called minimax rule.

3.3 Neyman–Pearson approach

In many problems of practical interest the imposition of a specific cost structure on the decisions made is not possible or desirable. The Neyman–Pearson approach involves a trade-off between the two types of errors that one can make in choosing a particular hypothesis. The Neyman–Pearson design criterion is to maximize the power of the test (probability of detection) subject to a chosen significance of the test (false alarm probability).

3.4 Likelihood ratio test

It is remarkable that all three very different approaches – Bayesian, minimax, and Neyman–Pearson – lead to the same test called the likelihood ratio test [28] . The likelihood ratio $\Lambda$  is the ratio of the pdf when the signal is present to the pdf when it is absent:

$$\begin{array}{c}\Lambda \left(x\right):=\frac{p_1\left(x\right)}{p_0\left(x\right)}.\end{array}$$

(17)

We accept the hypothesis $H_1$  if $\Lambda >k$  , where $k$  is the threshold that is calculated from the costs $C_{ij}$  , priors ${\pi }_i$  , or the significance of the test depending on what approach is being used.

3.4.1 Gaussian case – The matched filter

Let $h$  be the gravitational-wave signal and let $n$  be the detector noise. For convenience we assume that the signal $h$  is a continuous function of time $t$  and that the noise $n$  is a continuous random process. Results for the discrete time data that we have in practice can then be obtained by a suitable sampling of the continuous-in-time expressions. Assuming that the noise is additive the data $x$  can be written as

$$\begin{array}{c}x\left(t\right)=n\left(t\right)+h\left(t\right).\end{array}$$

(18)

In addition, if the noise is a zero-mean, stationary, and Gaussian random process, the log likelihood function is given by

$$\begin{array}{c}log\Lambda =\left(x\right|h)-\frac{1}{2}(h\left|h\right),\end{array}$$

(19)

where the scalar product $(\cdot |\cdot )$  is defined by

$$\begin{array}{c}\left(x\right|y):=4\Re {\int }_0^{\infty }\frac{\tilde{x}\left(f\right){\tilde{y}}^{*}\left(f\right)}{\tilde{S}\left(f\right)}df.\end{array}$$

(20)

In Equation ( 20 ) $\Re$  denotes the real part of a complex expression, the tilde denotes the Fourier transform, the asterisk is complex conjugation, and $\tilde{S}$  is the one-sided spectral density of the noise in the detector, which is defined through equation

$$\begin{array}{c}E\left[\tilde{n}\left(f\right){\tilde{n}}^{*}\left(f^{\prime }\right)\right]=\frac{1}{2}\delta (f-f^{\prime })\tilde{S}\left(f\right),\end{array}$$

(21)

where E denotes the expectation value.

From the expression ( 19 ) we see immediately that the likelihood ratio test consists of correlating the data $x$  with the signal $h$  that is present in the noise and comparing the correlation to a threshold.

Such a correlation is called the matched filter. The matched filter is a linear operation on the data.

An important quantity is the optimal signal-to-noise ratio $\rho$  defined by

$$\begin{array}{c}{\rho }^2:=\left(h\right|h)=4\Re {\int }_0^{\infty }\frac{\left|\tilde{h}\right(f){|}^2}{\tilde{S}\left(f\right)}df.\end{array}$$

(22)

We see in the following that $\rho$  determines the probability of detection of the signal. The higher the signal-to-noise ratio the higher the probability of detection.

An interesting property of the matched filter is that it maximizes the signal-to-noise ratio over all linear filters [28] . This property is independent of the probability distribution of the noise.

4 Parameter Estimation

Very often we know the waveform of the signal that we are searching for in the data in terms of a finite number of unknown parameters. We would like to find optimal procedures of estimating these parameters. An estimator of a parameter $\theta$  is a function $\hat{\theta }\left(x\right)$  that assigns to each data the “best” guess of the true value of $\theta$  . Note that because $\hat{\theta }\left(x\right)$  depends on the random data it is a random variable. Ideally we would like our estimator to be (i) unbiased, i.e., its expectation value to be equal to the true value of the parameter, and (ii) of minimum variance. Such estimators are rare and in general difficult to find. As in the signal detection there are several approaches to the parameter estimation problem. The subject is exposed in detail in reference [54] . See also [88] for a concise account.

4.1 Bayesian estimation

We assign a cost function $C({\theta }^{\prime },\theta )$  of estimating the true value of $\theta$  as ${\theta }^{\prime }$  . We then associate with an estimator $\hat{\theta }$  a conditional risk or cost averaged over all realizations of data $x$  for each value of the parameter $\theta$  :

$$\begin{array}{c}{R}_{\theta }\left(\hat{\theta }\right)=E_{\theta }\left[C\right(\hat{\theta },\theta \left)\right]={\int }_{X}C\left(\hat{\theta }\left(x\right),\theta \right)p(x,\theta )dx,\end{array}$$

(23)

where $X$  is the set of observations and $p(x,\theta )$  is the joint probability distribution of data $x$  and parameter $\theta$  . We further assume that there is a certain a priori probability distribution $\pi \left(\theta \right)$  of the parameter $\theta$  . We then define the Bayes estimator as the estimator that minimizes the average risk defined as

$$\begin{array}{c}r\left(\hat{\theta }\right)=E\left[R_{\theta }\right(\hat{\theta }\left)\right]={\int }_X{\int }_{\Theta }C\left(\hat{\theta }\left(x\right),\theta \right)p(x,\theta )\pi \left(\theta \right)d\theta dx,\end{array}$$

(24)

where E is the expectation value with respect to an a priori distribution $\pi$  , and $\Theta$  is the set of observations of the parameter $\theta$  . It is not difficult to show that for a commonly used cost function

$$\begin{array}{c}C({\theta }^{\prime },\theta )=({\theta }^{\prime }-\theta {)}^2,\end{array}$$

(25)

the Bayesian estimator is the conditional mean of the parameter $\theta$  given data $x$  , i.e.,

$$\begin{array}{c}\hat{\theta }\left(x\right)=E\left[\theta \right|x]={\int }_{\Theta }\hat{\theta }(x\left)p\right(\theta \left|x\right)d\theta ,\end{array}$$

(26)

where $p\left(\theta \right|x)$  is the conditional probability density of parameter $\theta$  given the data $x$  .

4.2 Maximum a posteriori probability estimation

Suppose that in a given estimation problem we are not able to assign a particular cost function $C({\theta }^{\prime },\theta )$  . Then a natural choice is a uniform cost function equal to $0$  over a certain interval $I_{\theta }$  of the parameter $\theta$  . From Bayes theorem [17] we have

$$\begin{array}{c}p\left(\theta \right|x)=\frac{p(x,\theta )\pi \left(\theta \right)}{p\left(x\right)},\end{array}$$

(27)

where $p\left(x\right)$  is the probability distribution of data $x$  . Then from Equation ( 24 ) one can deduce that for each data $x$  the Bayes estimate is any value of $\theta$  that maximizes the conditional probability $p\left(\theta \right|x)$  . The density $p\left(\theta \right|x)$  is also called the a posteriori probability density of parameter $\theta$  and the estimator that maximizes $p\left(\theta \right|x)$  is called the maximum a posteriori (MAP) estimator. It is denoted by ${\hat{\theta }}_{MAP}$  . We find that the MAP estimators are solutions of the following equation

$$\begin{array}{c}\frac{\partial logp(x,\theta )}{\partial \theta }=-\frac{\partial log\pi \left(\theta \right)}{\partial \theta },\end{array}$$

(28)

which is called the MAP equation.

4.3 Maximum likelihood estimation

Often we do not know the a priori probability density of a given parameter and we simply assign to it a uniform probability. In such a case maximization of the a posteriori probability is equivalent to maximization of the probability density $p(x,\theta )$  treated as a function of $\theta$  . We call the function $l(\theta ,x):=p(x,\theta )$  the likelihood function and the value of the parameter $\theta$  that maximizes $l(\theta ,x)$  the maximum likelihood (ML) estimator. Instead of the function $l$  we can use the function $\Lambda (\theta ,x)=l(\theta ,x)/p\left(x\right)$  (assuming that $p\left(x\right)>0$  ). $\Lambda$  is then equivalent to the likelihood ratio [see Equation ( 17 )] when the parameters of the signal are known. Then the ML estimators are obtained by solving the equation

$$\begin{array}{c}\frac{\partial log\Lambda (\theta ,x)}{\partial \theta }=0,\end{array}$$

(29)

which is called the ML equation.

4.3.1 Gaussian case

For the general gravitational-wave signal defined in Equation ( 12 ) the log likelihood function is given by

$$\begin{array}{c}log\Lambda ={\mathbf{a}}^T\cdot \mathbf{N}-\frac{1}{2}{\mathbf{a}}^T\cdot \mathsf{M}\cdot \mathbf{a},\end{array}$$

(30)

where the components of the column vector $\mathbf{N}$  and the matrix $\mathsf{M}$  are given by

$$\begin{array}{c}{N}^{\left(k\right)}:=\left(x\right|h^{\left(k\right)}),{\mathsf{M}}^{\left(k\right)\left(l\right)}:=(h^{\left(k\right)}\left|h^{\left(l\right)}\right),\end{array}$$

(31)

with $x\left(t\right)=n\left(t\right)+h\left(t\right)$  , and where $n\left(t\right)$  is a zero-mean Gaussian random process. The ML equations for the extrinsic parameters $\mathbf{a}$  can be solved explicitly and their ML estimators $\hat{\mathbf{a}}$  are given by

$$\begin{array}{c}\hat{\mathbf{a}}={\mathsf{M}}^{-1}\cdot \mathbf{N}.\end{array}$$

(32)

Substituting $\hat{\mathbf{a}}$  into $log\Lambda$  we obtain a function

$$\begin{array}{c}\mathcal{ℱ}=\frac{1}{2}{\mathbf{N}}^T\cdot {\mathsf{M}}^{-1}\cdot \mathbf{N},\end{array}$$

(33)

that we call the $\mathcal{ℱ}$  -statistic. The $\mathcal{ℱ}$  -statistic depends (nonlinearly) only on the intrinsic parameters ${\xi }^{\mu }$  .

Thus the procedure to detect the signal and estimate its parameters consists of two parts. The first part is to find the (local) maxima of the $\mathcal{ℱ}$  -statistic in the intrinsic parameter space. The ML estimators of the intrinsic parameters are those for which the $\mathcal{ℱ}$  -statistic attains a maximum. The second part is to calculate the estimators of the extrinsic parameters from the analytic formula ( 32 ), where the matrix $\mathsf{M}$  and the correlations $\mathbf{N}$  are calculated for the intrinsic parameters equal to their ML estimators obtained from the first part of the analysis. We call this procedure the maximum likelihood detection. See Section  4.8 for a discussion of the algorithms to find the (local) maxima of the $\mathcal{ℱ}$  -statistic.

4.4 Fisher information

It is important to know how good our estimators are. We would like our estimator to have as small variance as possible. There is a useful lower bound on variances of the parameter estimators called Cramèr–Rao bound. Let us first introduce the Fisher information matrix $\Gamma$  with the components defined by

$$\begin{array}{c}{\Gamma }_{ij}:=E\left[\frac{\partial ln\Lambda }{\partial {\theta }_i}\frac{\partial ln\Lambda }{\partial {\theta }_j}\right]=-E\left[\frac{{\partial }^{2}ln\Lambda }{\partial {\theta }_i\partial {\theta }_j}\right].\end{array}$$

(34)

The Cramèr–Rao bound states that for unbiased estimators the covariance matrix of the estimators $C\ge {\Gamma }^{-1}$  . (The inequality $A\ge B$  for matrices means that the matrix $A-B$  is nonnegative definite.) A very important property of the ML estimators is that asymptotically (i.e., for a signal-to-noise ratio tending to infinity) they are (i) unbiased, and (ii) they have a Gaussian distribution with covariance matrix equal to the inverse of the Fisher information matrix.

4.4.1 Gaussian case

In the case of Gaussian noise the components of the Fisher matrix are given by

$$\begin{array}{c}{\Gamma }_{ij}=\left(\frac{\partial h}{\partial {\theta }_i}|\frac{\partial h}{\partial {\theta }_j}\right).\end{array}$$

(35)

For the case of the general gravitational-wave signal defined in Equation ( 12 ) the set of the signal parameters $\mathit{\theta }$  splits naturally into extrinsic and intrinsic parameters: $\mathit{\theta }=(a^{\left(k\right)},{\xi }^{\mu })$  . Then the Fisher matrix can be written in terms of block matrices for these two sets of parameters as

$$\begin{array}{c}\Gamma =\left(\begin{array}{cc}\mathsf{M}& \mathsf{F}\cdot \mathbf{a}\\ {\mathbf{a}}^T\cdot {\mathsf{F}}^T& {\mathbf{a}}^T\cdot \mathsf{S}\cdot \mathbf{a}\end{array}\right),\end{array}$$

(36)

where the top left block corresponds to the extrinsic parameters, the bottom right block corresponds to the intrinsic parameters, the superscript T denotes here transposition over the extrinsic parameter indices, and the dot stands for the matrix multiplication with respect to these parameters. Matrix $\mathsf{M}$  is given by Equation ( 31 ), and the matrices $\mathsf{F}$  and $\mathsf{S}$  are defined as follows:

$$\begin{array}{c}{F}_{\mu }^{\left(k\right)\left(l\right)}:=\left(h^{\left(k\right)}|\frac{\partial h^{\left(l\right)}}{\partial {\xi }^{\mu }}\right),S_{\mu \nu }^{\left(k\right)\left(l\right)}:=\left(\frac{\partial h^{\left(k\right)}}{\partial {\xi }^{\mu }}|\frac{\partial h^{\left(l\right)}}{\partial {\xi }^{\nu }}\right).\end{array}$$

(37)

The covariance matrix $\mathsf{C}$  , which approximates the expected covariances of the ML parameter estimators, is defined as ${\Gamma }^{-1}$  . Using the standard formula for the inverse of a block matrix [58] we have

$$\begin{array}{c}\mathsf{C}=\left(\begin{array}{ccc}{\mathsf{M}}^{-1}+{\mathsf{M}}^{-1}\cdot (\mathsf{F}\cdot \mathbf{a})\cdot {\overline{\Gamma }}^{-1}\cdot (\mathsf{F}\cdot \mathbf{a}{)}^T\cdot {\mathsf{M}}^{-1}& & -{\mathsf{M}}^{-1}\cdot (\mathsf{F}\cdot \mathbf{a})\cdot {\overline{\Gamma }}^{-1}\\ -{\overline{\Gamma }}^{-1}\cdot (\mathsf{F}\cdot \mathbf{a}{)}^T\cdot {\mathsf{M}}^{-1}& & {\overline{\Gamma }}^{-1}\end{array}\right),\end{array}$$

(38)

where

$$\begin{array}{c}\overline{\Gamma }:={\mathbf{a}}^T\cdot (\mathsf{S}-{\mathsf{F}}^T\cdot {\mathsf{M}}^{-1}\cdot \mathsf{F})\cdot \mathbf{a}.\end{array}$$

(39)

We call ${\overline{\Gamma }}^{\mu \nu }$  (the Schur complement of $\mathsf{M}$  ) the projected Fisher matrix (onto the space of intrinsic parameters). Because the projected Fisher matrix is the inverse of the intrinsic-parameter submatrix of the covariance matrix $\mathsf{C}$  , it expresses the information available about the intrinsic parameters that takes into account the correlations with the extrinsic parameters. Note that ${\overline{\Gamma }}^{\mu \nu }$  is still a function of the putative extrinsic parameters.

We next define the normalized projected Fisher matrix

$$\begin{array}{c}{\overline{\Gamma }}_n:=\frac{\overline{\Gamma }}{{\rho }^2}=\frac{{\mathbf{a}}^T\cdot (\mathsf{S}-{\mathsf{F}}^T\cdot {\mathsf{M}}^{-1}\cdot \mathsf{F})\cdot \mathbf{a}}{{\mathbf{a}}^T\cdot \mathsf{M}\cdot \mathbf{a}},\end{array}$$

(40)

where $\rho =\sqrt{{\mathbf{a}}^T\cdot \mathsf{M}\cdot \mathbf{a}}$  is the signal-to-noise ratio. From the Rayleigh principle [58] follows that the minimum value of the component ${\overline{\Gamma }}_n^{\mu \nu }$  is given by the smallest eigenvalue (taken with respect to the extrinsic parameters) of the matrix ${\left((\mathsf{S}-{\mathsf{F}}^T\cdot {\mathsf{M}}^{-1}\cdot \mathsf{F})\cdot {\mathsf{M}}^{-1}\right)}^{\mu \nu }$  . Similarly, the maximum value of the component ${\overline{\Gamma }}_n^{\mu \nu }$  is given by the largest eigenvalue of that matrix. Because the trace of a matrix is equal to the sum of its eigenvalues, the matrix

$$\begin{array}{c}\tilde{\Gamma }:=\frac{1}{4}tr\left[\left(\mathsf{S}-{\mathsf{F}}^T\cdot {\mathsf{M}}^{-1}\cdot \mathsf{F}\right)\cdot {\mathsf{M}}^{-1}\right],\end{array}$$

(41)

where the trace is taken over the extrinsic-parameter indices, expresses the information available about the intrinsic parameters, averaged over the possible values of the extrinsic parameters. Note that the factor 1/4 is specific to the case of four extrinsic parameters. We call ${\tilde{\Gamma }}^{\mu \nu }$  the reduced Fisher matrix. This matrix is a function of the intrinsic parameters alone. We see that the reduced Fisher matrix plays a key role in the signal processing theory that we review here. It is used in the calculation of the threshold for statistically significant detection and in the formula for the number of templates needed to do a given search.

For the case of the signal

$$\begin{array}{c}h(t;A_0,{\phi }_0,{\xi }^{\mu })=A_{0}g(t;{\xi }^{\mu })cos\left(\phi (t;{\xi }^{\mu })-{\phi }_0\right),\end{array}$$

(42)

the normalized projected Fisher matrix ${\overline{\Gamma }}_n$  is independent of the extrinsic parameters $A_0$  and ${\phi }_0$  , and it is equal to the reduced matrix $\tilde{\Gamma }$  [66] . The components of $\tilde{\Gamma }$  are given by

$$\begin{array}{c}{\tilde{\Gamma }}^{\mu \nu }={\Gamma }_0^{\mu \nu }-\frac{{\Gamma }_0^{{\phi }_0\mu }{\Gamma }_0^{{\phi }_0\nu }}{{\Gamma }_0^{{\phi }_0{\phi }_0}},\end{array}$$

(43)

where ${\Gamma }_0^{ij}$  is the Fisher matrix for the signal $g(t;{\xi }^{\mu })cos\left(\phi (t;{\xi }^{\mu })-{\phi }_0\right)$  .

4.5 False alarm and detection probabilities – Gaussian case

4.5.1 Statistical properties of the $\mathcal{ℱ}$  -statistic

We first present the false alarm and detection pdfs when the intrinsic parameters of the signal are known. In this case the statistic $\mathcal{ℱ}$  is a quadratic form of the random variables that are correlations of the data. As we assume that the noise in the data is Gaussian and the correlations are linear functions of the data, $\mathcal{ℱ}$  is a quadratic form of the Gaussian random variables. Consequently $\mathcal{ℱ}$  -statistic has a distribution related to the ${\chi }^2$  distribution. One can show (see Section III B in [41] ) that for the signal given by Equation ( 12 ), $2\mathcal{ℱ}$  has a ${\chi }^2$  distribution with 4 degrees of freedom when the signal is absent and noncentral ${\chi }^2$  distribution with 4 degrees of freedom and non-centrality parameter equal to signal-to-noise ratio $\left(h\right|h)$  when the signal is present.

As a result the pdfs $p_0$  and $p_1$  of $\mathcal{ℱ}$  when the intrinsic parameters are known and when respectively the signal is absent and present are given by

$$\begin{array}{ccc}{p}_0(\mathcal{ℱ})& =& \frac{{\mathcal{ℱ}}^{n/2-1}}{(n/2-1)!}exp(-\mathcal{ℱ}),\end{array}$$

(44)

$$\begin{array}{ccc}{p}_1(\rho ,\mathcal{ℱ})& =& \frac{(2\mathcal{ℱ}{)}^{(n/2-1)/2}}{{\rho }^{n/2-1}}{I}_{n/2-1}\left(\rho \sqrt{2\mathcal{ℱ}}\right)exp\left(-\mathcal{ℱ}-\frac{1}{2}{\rho }^2\right),\end{array}$$

(45)

where $n$  is the number of degrees of freedom of ${\chi }^2$  distributions and $I_{n/2-1}$  is the modified Bessel function of the first kind and order $n/2-1$  . The false alarm probability $P_F$  is the probability that $\mathcal{ℱ}$  exceeds a certain threshold ${\mathcal{ℱ}}_0$  when there is no signal. In our case we have

$$\begin{array}{c}{P}_F\left({\mathcal{ℱ}}_0\right):={\int }_{{\mathcal{ℱ}}_0}^{\infty }{p}_0(\mathcal{ℱ})d\mathcal{ℱ}=exp(-{\mathcal{ℱ}}_0){\sum }_{k=0}^{n/2-1}\frac{{\mathcal{ℱ}}_0^k}{k!}.\end{array}$$

(46)

The probability of detection $P_D$  is the probability that $\mathcal{ℱ}$  exceeds the threshold ${\mathcal{ℱ}}_0$  when the signal-to-noise ratio is equal to $\rho$  :

$$\begin{array}{c}{P}_D(\rho ,{\mathcal{ℱ}}_0):={\int }_{{\mathcal{ℱ}}_0}^{\infty }{p}_1(\rho ,\mathcal{ℱ})d\mathcal{ℱ}.\end{array}$$

(47)

The integral in the above formula can be expressed in terms of the generalized Marcum $Q$  -function [81, 37] , $Q(\alpha ,\beta )=P_D(\alpha ,{\beta }^2/2)$  . We see that when the noise in the detector is Gaussian and the intrinsic parameters are known, the probability of detection of the signal depends on a single quantity: the optimal signal-to-noise ratio $\rho$  .

4.5.2 False alarm probability

Next we return to the case when the intrinsic parameters $\mathit{\xi }$  are not known. Then the statistic $\mathcal{ℱ}\left(\mathit{\xi }\right)$  given by Equation ( 33 ) is a certain generalized multiparameter random process called the random field (see Adler's monograph [4] for a comprehensive discussion of random fields). If the vector $\mathit{\xi }$  has one component the random field is simply a random process. For random fields we can define the autocovariance function $\mathcal{C}$  just in the same way as we define such a function for a random process:

$$\begin{array}{c}\mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime }):=E_0[\mathcal{ℱ}(\mathit{\xi })\mathcal{ℱ}(\mathit{\xi }{}^{\prime }\left)\right]-E_0[\mathcal{ℱ}(\mathit{\xi }\left)\right]E_0[\mathcal{ℱ}(\mathit{\xi }{}^{\prime }\left)\right],\end{array}$$

(48)

where $\mathit{\xi }$  and $\mathit{\xi }{}^{\prime }$  are two values of the intrinsic parameter set, and $E_0$  is the expectation value when the signal is absent. One can show that for the signal ( 12 ) the autocovariance function $\mathcal{C}$  is given by

$$\begin{array}{c}\mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime })=\frac{1}{4}tr\left({\mathsf{Q}}^T\cdot {\mathsf{M}}^{-1}\cdot \mathsf{Q}\cdot {{\mathsf{M}}^{\prime }}^{-1}\right),\end{array}$$

(49)

where

$$\begin{array}{c}{\mathsf{Q}}^{\left(k\right)\left(l\right)}:=\left(h^{\left(k\right)}(t;\mathit{\xi })\left|h^{\left(l\right)}\right(t;\mathit{\xi }{}^{\prime })\right),{{\mathsf{M}}^{\prime }}^{\left(k\right)\left(l\right)}:=\left(h^{\left(k\right)}(t;\mathit{\xi }{}^{\prime })\left|h^{\left(l\right)}\right(t;\mathit{\xi }{}^{\prime })\right).\end{array}$$

(50)

We have $\mathcal{C}(\mathit{\xi },\mathit{\xi })=1$  .

One can estimate the false alarm probability in the following way [42] . The autocovariance function $\mathcal{C}$  tends to zero as the displacement $\Delta \mathit{\xi }=\mathit{\xi }{}^{\prime }-\mathit{\xi }$  increases (it is maximal for $\Delta \mathit{\xi }=0$  ). Thus we can divide the parameter space into elementary cells such that in each cell the autocovariance function $\mathcal{C}$  is appreciably different from zero. The realizations of the random field within a cell will be correlated (dependent), whereas realizations of the random field within each cell and outside the cell are almost uncorrelated (independent). Thus the number of cells covering the parameter space gives an estimate of the number of independent realizations of the random field. The correlation hypersurface is a closed surface defined by the requirement that at the boundary of the hypersurface the correlation $\mathcal{C}$  equals half of its maximum value. The elementary cell is defined by the equation

$$\begin{array}{c}\mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime })=\frac{1}{2}\end{array}$$

(51)

for $\mathit{\xi }$  at cell center and $\mathit{\xi }{}^{\prime }$  on cell boundary. To estimate the number of cells we perform the Taylor expansion of the autocorrelation function up to the second-order terms:

$$\begin{array}{c}\mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime })\sim =1+\frac{\partial \mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime })}{\partial {\xi }_i^{\prime }}{|}_{\mathit{\xi }{}^{\prime }=\mathit{\xi }}\Delta {\xi }_i+\frac{1}{2}\frac{{\partial }^2\mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime })}{\partial {\xi }_i^{\prime }\partial {\xi }_j^{\prime }}{|}_{\mathit{\xi }{}^{\prime }=\mathit{\xi }}\Delta {\xi }_i\Delta {\xi }_j.\end{array}$$

(52)

As $\mathcal{C}$  attains its maximum value when $\mathit{\xi }-\mathit{\xi }{}^{\prime }=0$  , we have

$$\begin{array}{c}\frac{\partial \mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime })}{\partial {\xi }_i^{\prime }}{|}_{\mathit{\xi }{}^{\prime }=\mathit{\xi }}=0.\end{array}$$

(53)

Let us introduce the symmetric matrix

$$\begin{array}{c}{G}_{ij}:=-\frac{1}{2}\frac{{\partial }^2\mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime })}{\partial {\xi }_i^{\prime }\partial {\xi }_j^{\prime }}{|}_{{\mathit{\xi }}^{\mathsf{\prime }}=\mathit{\xi }}.\end{array}$$

(54)

Then the approximate equation for the elementary cell is given by

$$\begin{array}{c}{G}_{ij}\Delta {\xi }_i\Delta {\xi }_j=\frac{1}{2}.\end{array}$$

(55)

It is interesting to find a relation between the matrix $G$  and the Fisher matrix. One can show (see [51] , Appendix B) that the matrix $G$  is precisely equal to the reduced Fisher matrix $\tilde{\Gamma }$  given by Equation ( 41 ). Let $K$  be the number of the intrinsic parameters. If the components of the matrix $G$  are constant (independent of the values of the parameters of the signal) the above equation is an equation for a hyperellipse. The $K$  -dimensional Euclidean volume $V_{cell}$  of the elementary cell defined by Equation ( 55 ) equals

$$\begin{array}{c}{V}_{cell}=\frac{(\pi /2{)}^{K/2}}{\Gamma (K/2+1)\sqrt{detG}},\end{array}$$

(56)

where $\Gamma$  denotes the Gamma function. We estimate the number $N_c$  of elementary cells by dividing the total Euclidean volume $V$  of the $K$  -dimensional parameter space by the volume $V_{cell}$  of the elementary cell, i.e. we have

$$\begin{array}{c}{N}_c=\frac{V}{V_{cell}}.\end{array}$$

(57)

The components of the matrix $G$  are constant for the signal $h(t;A_0,{\phi }_0,{\xi }^{\mu })=A_{0}cos\left(\phi (t;{\xi }^{\mu })-{\phi }_0\right)$  when the phase $\phi (t;{\xi }^{\mu })$  is a linear function of the intrinsic parameters ${\xi }^{\mu }$  .

To estimate the number of cells in the case when the components of the matrix $G$  are not constant, i.e. when they depend on the values of the parameters, we write Equation ( 57 ) as

$$\begin{array}{c}{N}_c=\frac{\Gamma (K/2+1)}{(\pi /2{)}^{K/2}}{\int }_V\sqrt{detG}dV.\end{array}$$

(58)

This procedure can be thought of as interpreting the matrix $G$  as the metric on the parameter space. This interpretation appeared for the first time in the context of gravitational-wave data analysis in the work by Owen [66] , where an analogous integral formula was proposed for the number of templates needed to perform a search for gravitational-wave signals from coalescing binaries.

The concept of number of cells was introduced in [42] and it is a generalization of the idea of an effective number of samples introduced in [30] for the case of a coalescing binary signal.

We approximate the probability distribution of $\mathcal{ℱ}\left(\mathit{\xi }\right)$  in each cell by the probability $p_0(\mathcal{ℱ})$  when the parameters are known [in our case by probability given by Equation ( 44 )]. The values of the statistic $\mathcal{ℱ}$  in each cell can be considered as independent random variables. The probability that $\mathcal{ℱ}$  does not exceed the threshold ${\mathcal{ℱ}}_0$  in a given cell is $1-P_F\left({\mathcal{ℱ}}_0\right)$  , where $P_F\left({\mathcal{ℱ}}_0\right)$  is given by Equation ( 46 ). Consequently the probability that $\mathcal{ℱ}$  does not exceed the threshold ${\mathcal{ℱ}}_0$  in all the $N_c$  cells is $[1-P_F({\mathcal{ℱ}}_0){]}^{N_c}$  . The probability $P_F^T$  that $\mathcal{ℱ}$  exceeds ${\mathcal{ℱ}}_0$  in one or more cell is thus given by

$$\begin{array}{c}{P}_F^T\left({\mathcal{ℱ}}_0\right)=1-[1-P_F({\mathcal{ℱ}}_0){]}^{N_c}.\end{array}$$

(59)

This by definition is the false alarm probability when the phase parameters are unknown. The number of false alarms $N_F$  is given by

$$\begin{array}{c}{N}_F=N_c{P}_F^T\left({\mathcal{ℱ}}_0\right).\end{array}$$

(60)

A different approach to the calculation of the number of false alarms using the Euler characteristic of level crossings of a random field is described in [41] .

It was shown (see [25] ) that for any finite ${\mathcal{ℱ}}_0$  and $N_c$  , Equation ( 59 ) provides an upper bound for the false alarm probability. Also in [25] a tighter upper bound for the false alarm probability was derived by modifying a formula obtained by Mohanty [59] . The formula amounts essentially to introducing a suitable coefficient multiplying the number of cells $N_c$  .

4.5.3 Detection probability

When the signal is present a precise calculation of the pdf of $\mathcal{ℱ}$  is very difficult because the presence of the signal makes the data random process $x\left(t\right)$  non-stationary. As a first approximation we can estimate the probability of detection of the signal when the parameters are unknown by the probability of detection when the parameters of the signal are known [given by Equation ( 47 )].

This approximation assumes that when the signal is present the true values of the phase parameters fall within the cell where $\mathcal{ℱ}$  has a maximum. This approximation will be the better the higher the signal-to-noise ratio $\rho$  is.

4.6 Number of templates

To search for gravitational-wave signals we evaluate the $\mathcal{ℱ}$  -statistic on a grid in parameter space. The grid has to be sufficiently fine such that the loss of signals is minimized. In order to estimate the number of points of the grid, or in other words the number of templates that we need to search for a signal, the natural quantity to study is the expectation value of the $\mathcal{ℱ}$  -statistic when the signal is present. We have

$$\begin{array}{c}E[\mathcal{ℱ}]=\frac{1}{2}\left(4+{\mathbf{a}}^T\cdot {\mathsf{Q}}^T\cdot {{\mathsf{M}}^{\prime }}^{-1}\cdot \mathsf{Q}\cdot \mathbf{a}\right).\end{array}$$

(61)

The components of the matrix $\mathsf{Q}$  are given in Equation ( 50 ). Let us rewrite the expectation value ( 61 ) in the following form,

$$\begin{array}{c}E[\mathcal{ℱ}]=\frac{1}{2}\left(4+{\rho }^2\frac{{\mathbf{a}}^T\cdot {\mathsf{Q}}^T\cdot {{\mathsf{M}}^{\prime }}^{-1}\cdot \mathsf{Q}\cdot \mathbf{a}}{{\mathbf{a}}^T\cdot \mathsf{M}\cdot \mathbf{a}}\right),\end{array}$$

(62)

where $\rho$  is the signal-to-noise ratio. Let us also define the normalized correlation function

$$\begin{array}{c}{\mathcal{C}}_n:=\frac{{\mathbf{a}}^T\cdot {\mathsf{Q}}^T\cdot {{\mathsf{M}}^{\prime }}^{-1}\cdot \mathsf{Q}\cdot \mathbf{a}}{{\mathbf{a}}^T\cdot \mathsf{M}\cdot \mathbf{a}}.\end{array}$$

(63)

From the Rayleigh principle [58] it follows that the minimum of the normalized correlation function is equal to the smallest eigenvalue of the normalized matrix ${\mathsf{Q}}^T\cdot {{\mathsf{M}}^{\prime }}^{-1}\cdot \mathsf{Q}\cdot {\mathsf{M}}^{-1}$  , whereas the maximum is given by its largest eigenvalue. We define the reduced correlation function as

$$\begin{array}{c}\mathcal{C}(\mathit{\xi },\mathit{\xi }{}^{\prime }):=\frac{1}{4}tr\left({\mathsf{Q}}^T\cdot {\mathsf{M}}^{-1}\cdot \mathsf{Q}\cdot {{\mathsf{M}}^{\prime }}^{-1}\right).\end{array}$$

(64)

As the trace of a matrix equals the sum of its eigenvalues, the reduced correlation function $\mathcal{C}$  is equal to the average of the eigenvalues of the normalized correlation function ${\mathcal{C}}_n$  . In this sense we can think of the reduced correlation function as an “average” of the normalized correlation function. The advantage of the reduced correlation function is that it depends only on the intrinsic parameters $\mathit{\xi }$  , and thus it is suitable for studying the number of grid points on which the $\mathcal{ℱ}$  -statistic needs to be evaluated. We also note that the normalized correlation function $\mathcal{C}$  precisely coincides with the autocovariance function $\mathcal{C}$  of the $\mathcal{ℱ}$  -statistic given by Equation ( 49 ).

Like in the calculation of the number of cells in order to estimate the number of templates we perform a Taylor expansion of $\mathcal{C}$  up to second order terms around the true values of the parameters, and we obtain an equation analogous to Equation ( 55 ),

$$\begin{array}{c}{G}_{ij}\Delta {\xi }_i\Delta {\xi }_j=1-C_0,\end{array}$$

(65)

where $G$  is given by Equation ( 54 ). By arguments identical to those in deriving the formula for the number of cells we arrive at the following formula for the number of templates:

$$\begin{array}{c}{N}_t=\frac{1}{(1-C_0{)}^{K/2}}\frac{\Gamma (K/2+1)}{{\pi }^{K/2}}{\int }_V\sqrt{detG}dV.\end{array}$$

(66)

When $C_0=1/2$  the above formula coincides with the formula for the number $N_c$  of cells, Equation ( 58 ). Here we would like to place the templates sufficiently closely so that the loss of signals is minimized. Thus $1-C_0$  needs to be chosen sufficiently small. The formula ( 66 ) for the number of templates assumes that the templates are placed in the centers of hyperspheres and that the hyperspheres fill the parameter space without holes. In order to have a tiling of the parameter space without holes we can place the templates in the centers of hypercubes which are inscribed in the hyperspheres. Then the formula for the number of templates reads

$$\begin{array}{c}{N}_t=\frac{1}{(1-C_0{)}^{K/2}}\frac{K^{K/2}}{2^K}{\int }_V\sqrt{detG}dV.\end{array}$$

(67)

For the case of the signal given by Equation ( 14 ) our formula for number of templates is equivalent to the original formula derived by Owen [66] . Owen [66] has also introduced a geometric approach to the problem of template placement involving the identification of the Fisher matrix with a metric on the parameter space. An early study of the template placement for the case of coalescing binaries can be found in [72, 29, 16] . Applications of the geometric approach of Owen to the case of spinning neutron stars and supernova bursts are given in [20, 9] .

The problem of how to cover the parameter space with the smallest possible number of templates, such that no point in the parameter space lies further away from a grid point than a certain distance, is known in mathematical literature as the covering problem [24] . The maximum distance of any point to the next grid point is called the covering radius $R$  . An important class of coverings are lattice coverings. We define a lattice in $K$  -dimensional Euclidean space ${\mathbb{R}}^K$  to be the set of points including $0$  such that if $u$  and $v$  are lattice points, then also $u+v$  and $u-v$  are lattice points. The basic building block of a lattice is called the fundamental region. A lattice covering is a covering of ${\mathbb{R}}^K$  by spheres of covering radius $R$  , where the centers of the spheres form a lattice. The most important quantity of a covering is its thickness $\Theta$  defined as

$$\begin{array}{c}\Theta :=\frac{\text{volume of one}K\text{-dimensional sphere}}{\text{volume of the fundamental region}}.\end{array}$$

(68)

In the case of a two-dimensional Euclidean space the best covering is the hexagonal covering and its thickness $\simeq 1.21$  . For dimensions higher than 2 the best covering is not known. We know however the best lattice covering for dimensions $K\le 23$  . These are so-called $A_K^{*}$  lattices which have a thickness ${\Theta }_{A_K^{*}}$  equal to

$$\begin{array}{c}{\Theta }_{A_K^{*}}=V_K\sqrt{K+1}{\left(\frac{K(K+2)}{12(K+1)}\right)}^{K/2},\end{array}$$

(69)

where $V_K$  is the volume of the $K$  -dimensional sphere of unit radius.

For the case of spinning neutron stars a 3-dimensional grid was constructed consisting of prisms with hexagonal bases [13] . This grid has a thickness around 1.84 which is much better than the cubic grid which has thickness of approximately 2.72. It is worse than the best lattice covering which has the thickness around 1.46. The advantage of an $A_K^{*}$  lattice over the hypercubic lattice grows exponentially with the number of dimensions.

4.7 Suboptimal filtering

To extract signals from the noise one very often uses filters that are not optimal. We may have to choose an approximate, suboptimal filter because we do not know the exact form of the signal (this is almost always the case in practice) or in order to reduce the computational cost and to simplify the analysis. The most natural and simplest way to proceed is to use as our statistic the $\mathcal{ℱ}$  -statistic where the filters $h_k^{\prime }(t;\mathit{\zeta })$  are the approximate ones instead of the optimal ones matched to the signal. In general the functions $h_k^{\prime }(t;\mathit{\zeta })$  will be different from the functions $h_k(t;\mathit{\xi })$  used in optimal filtering, and also the set of parameters $\mathit{\zeta }$  will be different from the set of parameters $\mathit{\xi }$  in optimal filters. We call this procedure the suboptimal filtering and we denote the suboptimal statistic by ${\mathcal{ℱ}}_s$  .

We need a measure of how well a given suboptimal filter performs. To find such a measure we calculate the expectation value of the suboptimal statistic. We get

$$\begin{array}{c}E\left[{\mathcal{ℱ}}_s\right]=\frac{1}{2}\left(4+{\mathbf{a}}^T\cdot {\mathsf{Q}}_s^T\cdot {\mathsf{M}}_{s^{\prime }}^{-1}\cdot {\mathsf{Q}}_s\cdot \mathbf{a}\right),\end{array}$$

(70)

where

$$\begin{array}{c}\begin{array}{ccc}{\mathsf{M}}_{s^{\prime }}^{\left(k\right)\left(l\right)}& :=& \left({h^{\prime }}^{\left(k\right)}(t;\mathit{\zeta })\left|{h^{\prime }}^{\left(l\right)}\right(t;\mathit{\zeta })\right),\\ {\mathsf{Q}}_s^{\left(k\right)\left(l\right)}& :=& \left(h^{\left(k\right)}(t;\mathit{\xi })\left|{h^{\prime }}^{\left(l\right)}\right(t;\mathit{\zeta })\right).\end{array}\end{array}$$

(71)

Let us rewrite the expectation value $E\left[{\mathcal{ℱ}}_s\right]$  in the following form,

$$\begin{array}{c}E\left[{\mathcal{ℱ}}_s\right]=\frac{1}{2}\left(4+{\rho }^2\frac{{\mathbf{a}}^T\cdot {\mathsf{Q}}_s^T\cdot {\mathsf{M}}_{s^{\prime }}^{-1}\cdot {\mathsf{Q}}_s\cdot \mathbf{a}}{{\mathbf{a}}^T\cdot \mathsf{M}\cdot \mathbf{a}}\right),\end{array}$$

(72)

where $\rho$  is the optimal signal-to-noise ratio. The expectation value $E\left[{\mathcal{ℱ}}_s\right]$  reaches its maximum equal to $2+{\rho }^2/2$  when the filter is perfectly matched to the signal. A natural measure of the performance of a suboptimal filter is the quantity FF defined by

$$\begin{array}{c}FF:={max}_{(\mathbf{a},\mathit{\zeta })}\sqrt{\frac{{\mathbf{a}}^T\cdot {\mathsf{Q}}_s^T\cdot {\mathsf{M}}_{s^{\prime }}^{-1}\cdot {\mathsf{Q}}_s\cdot \mathbf{a}}{{\mathbf{a}}^T\cdot \mathsf{M}\cdot \mathbf{a}}}.\end{array}$$

(73)

We call the quantity FF the generalized fitting factor.

In the case of a signal given by

$$\begin{array}{c}s(t;A_0,\mathit{\xi })=A_{0}h(t;\mathit{\xi }),\end{array}$$

(74)

the generalized fitting factor defined above reduces to the fitting factor introduced by Apostolatos [7] :

$$\begin{array}{c}FF={max}_{\mathit{\zeta }}\frac{\left(h(t;\mathit{\xi })\left|h^{\prime }\right(t;\mathit{\zeta })\right)}{\sqrt{\left(h(t;\mathit{\xi })\left|h\right(t;\mathit{\xi })\right)}\sqrt{\left(h^{\prime }(t;\mathit{\zeta })\left|h^{\prime }\right(t;\mathit{\zeta })\right)}}.\end{array}$$

(75)

The fitting factor is the ratio of the maximal signal-to-noise ratio that can be achieved with suboptimal filtering to the signal-to-noise ratio obtained when we use a perfectly matched, optimal filter. We note that for the signal given by Equation ( 74 ), FF is independent of the value of the amplitude $A_0$  . For the general signal with 4 constant amplitudes it follows from the Rayleigh principle that the fitting factor is the maximum of the largest eigenvalue of the matrix ${\mathsf{Q}}^T\cdot {{\mathsf{M}}^{\prime }}^{-1}\cdot \mathsf{Q}\cdot {\mathsf{M}}^{-1}$  over the intrinsic parameters of the signal.

For the case of a signal of the form

$$\begin{array}{c}s(t;A_0,{\phi }_0,\mathit{\xi })=A_{0}cos\left(\phi (t;\mathit{\xi })+{\phi }_0\right),\end{array}$$

(76)

where ${\phi }_0$  is a constant phase, the maximum over ${\phi }_0$  in Equation ( 75 ) can be obtained analytically.

Moreover, assuming that over the bandwidth of the signal the spectral density of the noise is constant and that over the observation time $cos\left(\phi (t;\mathit{\xi })\right)$  oscillates rapidly, the fitting factor is approximately given by

$$\begin{array}{c}FF\sim ={max}_{\mathit{\zeta }}{\left[{\left({\int }_0^{T_0}cos\left(\phi (t;\mathit{\xi })-{\phi }^{\prime }(t;\mathit{\zeta })\right)dt\right)}^2+{\left({\int }_0^{T_0}sin\left(\phi (t;\mathit{\xi })-{\phi }^{\prime }(t;\mathit{\zeta })\right)dt\right)}^2\right]}^{1/2}.\end{array}$$

(77)

In designing suboptimal filters one faces the issue of how small a fitting factor one can accept.

A popular rule of thumb is accepting $FF=0.97$  . Assuming that the amplitude of the signal and consequently the signal-to-noise ratio decreases inversely proportional to the distance from the source this corresponds to 10% loss of the signals that would be detected by a matched filter.

Proposals for good suboptimal (search) templates for the case of coalescing binaries are given in [22, 78] and for the case spinning neutron stars in [41, 13] .

4.8 Algorithms to calculate the $\mathcal{ℱ}$  -statistic

4.8.1 The two-step procedure

In order to detect signals we search for threshold crossings of the $\mathcal{ℱ}$  -statistic over the intrinsic parameter space. Once we have a threshold crossing we need to find the precise location of the maximum of $\mathcal{ℱ}$  in order to estimate accurately the parameters of the signal. A satisfactory procedure is the two-step procedure. The first step is a coarse search where we evaluate $\mathcal{ℱ}$  on a coarse grid in parameter space and locate threshold crossings. The second step, called fine search, is a refinement around the region of parameter space where the maximum identified by the coarse search is located.

There are two methods to perform the fine search. One is to refine the grid around the threshold crossing found by the coarse search [61, 59, 78, 76] , and the other is to use an optimization routine to find the maximum of $\mathcal{ℱ}$  [41, 51] . As initial value to the optimization routine we input the values of the parameters found by the coarse search. There are many maximization algorithms available.

One useful method is the Nelder–Mead algorithm [52] which does not require computation of the derivatives of the function being maximized.

4.8.2 Evaluation of the $\mathcal{ℱ}$  -statistic

Usually the grid in parameter space is very large and it is important to calculate the optimum statistic as efficiently as possible. In special cases the $\mathcal{ℱ}$  -statistic given by Equation ( 33 ) can be further simplified. For example, in the case of coalescing binaries $\mathcal{ℱ}$  can be expressed in terms of convolutions that depend on the difference between the time-of-arrival (TOA) of the signal and the TOA parameter of the filter. Such convolutions can be efficiently computed using Fast Fourier Transforms (FFTs). For continuous sources, like gravitational waves from rotating neutron stars observed by ground-based detectors [41] or gravitational waves form stellar mass binaries observed by space-borne detectors [51] , the detection statistic $\mathcal{ℱ}$  involves integrals of the general form

$$\begin{array}{c}{\int }_0^{T_0}x\left(t\right)m(t;\omega ,{\tilde{\xi }}^{\mu })exp\left(i\omega {\phi }_{mod}(t;{\tilde{\xi }}^{\mu })\right)exp\left(i\omega t\right)dt,\end{array}$$

(78)

where ${\tilde{\xi }}^{\mu }$  are the intrinsic parameters excluding the frequency parameter $\omega$  , $m$  is the amplitude modulation function, and $\omega {\phi }_{mod}$  the phase modulation function. The amplitude modulation function is slowly varying comparing to the exponential terms in the integral ( 78 ). We see that the integral ( 78 ) can be interpreted as a Fourier transform (and computed efficiently with an FFT), if ${\phi }_{mod}=0$  and if $m$  does not depend on the frequency $\omega$  . In the long-wavelength approximation the amplitude function $m$  does not depend on the frequency. In this case Equation ( 78 ) can be converted to a Fourier transform by introducing a new time variable $t_b$  [75] ,

$$\begin{array}{c}{t}_b:=t+{\phi }_{mod}(t;{\tilde{\xi }}^{\mu }).\end{array}$$

(79)

Thus in order to compute the integral ( 78 ), for each set of the intrinsic parameters ${\tilde{\xi }}^{\mu }$  we multiply the data by the amplitude modulation function $m$  , resample according to Equation ( 79 ), and perform the FFT. In the case of LISA detector data when the amplitude modulation $m$  depends on frequency we can divide the data into several band-passed data sets, choosing the bandwidth for each set sufficiently small so that the change of $mexp\left(i\omega {\phi }_{mod}\right)$  is small over the band. In the integral ( 78 ) we can then use as the value of the frequency in the amplitude and phase modulation function the maximum frequency of the band of the signal (see [51] for details).

4.8.3 Comparison with the Cramèr–Rao bound

In order to test the performance of the maximization method of the $\mathcal{ℱ}$  statistic it is useful to perform Monte Carlo simulations of the parameter estimation and compare the variances of the estimators with the variances calculated from the Fisher matrix. Such simulations were performed for various gravitational-wave signals [46, 16, 41] . In these simulations we observe that above a certain signal-to-noise ratio, that we call the threshold signal-to-noise ratio, the results of the Monte Carlo simulations agree very well with the calculations of the rms errors from the inverse of the Fisher matrix. However, below the threshold signal-to-noise ratio they differ by a large factor.

This threshold effect is well-known in signal processing [82] . There exist more refined theoretical bounds on the rms errors that explain this effect, and they were also studied in the context of the gravitational-wave signal from a coalescing binary [64] . Here we present a simple model that explains the deviations from the covariance matrix and reproduces well the results of the Monte Carlo simulations. The model makes use of the concept of the elementary cell of the parameter space that we introduced in Section  4.5.2 . The calculation given below is a generalization of the calculation of the rms error for the case of a monochromatic signal given by Rife and Boorstyn [70] .

When the values of parameters of the template that correspond to the maximum of the functional $\mathcal{ℱ}$  fall within the cell in the parameter space where the signal is present, the rms error is satisfactorily approximated by the inverse of the Fisher matrix. However, sometimes as a result of noise the global maximum is in the cell where there is no signal. We then say that an outlier has occurred. In the simplest case we can assume that the probability density of the values of the outliers is uniform over the search interval of a parameter, and then the rms error is given by

$$\begin{array}{c}{\sigma }_{out}^2=\frac{{\Delta }^2}{12},\end{array}$$

(80)

where $\Delta$  is the length of the search interval for a given parameter. The probability that an outlier occurs will be the higher the lower the signal-to-noise ratio is. Let $q$  be the probability that an outlier occurs. Then the total variance ${\sigma }^2$  of the estimator of a parameter is the weighted sum of the two errors

$$\begin{array}{c}{\sigma }^2={\sigma }_{out}^{2}q+{\sigma }_{CR}^2(1-q),\end{array}$$

(81)

where ${\sigma }_{CR}$  is the rms errors calculated from the covariance matrix for a given parameter. One can show [41] that the probability $q$  can be approximated by the following formula:

$$\begin{array}{c}q=1-{\int }_0^{\infty }{p}_1(\rho ,\mathcal{ℱ}){\left({\int }_0^{\mathcal{ℱ}}{p}_0\left(y\right)dy\right)}^{N_c-1}d\mathcal{ℱ},\end{array}$$

(82)

where $p_0$  and $p_1$  are the probability density functions of false alarm and detection given by Equations ( 44 ) and ( 45 ), respectively, and where $N_c$  is the number of cells in the parameter space.

Equation ( 82 ) is in good but not perfect agreement with the rms errors obtained from the Monte Carlo simulations (see [41] ). There are clearly also other reasons for deviations from the Cramèr–Rao bound. One important effect (see [64] ) is that the functional $\mathcal{ℱ}$  has many local subsidiary maxima close to the global one. Thus for a low signal-to-noise the noise may promote the subsidiary maximum to a global one.

4.9 Upper limits

Detection of a signal is signified by a large value of the $\mathcal{ℱ}$  -statistic that is unlikely to arise from the noise-only distribution. If instead the value of $\mathcal{ℱ}$  is consistent with pure noise with high probability we can place an upper limit on the strength of the signal. One way of doing this is to take the loudest event obtained in the search and solve the equation

$$\begin{array}{c}{P}_D({\rho }_{UL},{\mathcal{ℱ}}_L)=\beta \end{array}$$

(83)

for signal-to-noise ratio ${\rho }_{UL}$  , where $P_D$  is the detection probability given by Equation ( 47 ), ${\mathcal{ℱ}}_L$  is the value of the $\mathcal{ℱ}$  -statistic corresponding to the loudest event, and $\beta$  is a chosen confidence [12, 1] .

Then ${\rho }_{UL}$  is the desired upper limit with confidence $\beta$  .

When gravitational-wave data do not conform to a Gaussian probability density assumed in Equation ( 47 ), a more accurate upper limit can be obtained by injecting the signals into the detector's data and thereby estimating the probability of detection $P_D$  [3] .

4.10 Network of detectors

Several gravitational-wave detectors can observe gravitational waves from the same source. For example a network of bar detectors can observe a gravitational-wave burst from the same supernova explosion, or a network of laser interferometers can detect the inspiral of the same compact binary system. The space-borne LISA detector can be considered as a network of three detectors that can make three independent measurements of the same gravitational-wave signal. Simultaneous observations are also possible among different types of detectors, for example a search for supernova bursts can be performed simultaneously by bar and laser detectors [14] .

We consider the general case of a network of detectors. Let $\mathbf{h}$  be the signal vector and let n be the noise vector of the network of detectors, i.e., the vector component $h_k$  is the response of the gravitational-wave signal in the $k$  th detector with noise $n_k$  . Let us also assume that each $n_k$  has zero mean. Assuming that the noise in all detectors is additive the data vector $\mathbf{x}$  can be written as

$$\begin{array}{c}\mathbf{x}\left(t\right)=\mathbf{n}\left(t\right)+\mathbf{h}\left(t\right).\end{array}$$

(84)

In addition if the noise is a stationary, Gaussian, and continuous random process the log likelihood function is given by

$$\begin{array}{c}log\Lambda =(\mathbf{x}|\mathbf{h})-\frac{1}{2}(\mathbf{h}|\mathbf{h}).\end{array}$$

(85)

In Equation ( 85 ) the scalar product $(\cdot |\cdot )$  is defined by

$$\begin{array}{c}(\mathbf{x}|\mathbf{y}):=4\Re {\int }_0^{\infty }{\tilde{\mathbf{x}}}^T{\tilde{\mathbf{S}}}^{-1}{\tilde{\mathbf{y}}}^{*}df,\end{array}$$

(86)

where $\tilde{\mathbf{S}}$  is the one-sided cross spectral density matrix of the noises of the detector network which is defined by (here E denotes the expectation value)

$$\begin{array}{c}E\left[\tilde{\mathbf{n}}\left(f\right){\tilde{\mathbf{n}}}^{*T}\left(f^{\prime }\right)\right]=\frac{1}{2}\delta (f-f^{\prime })\tilde{\mathbf{S}}\left(f\right).\end{array}$$

(87)

The analysis is greatly simplified if the cross spectrum matrix $\mathbf{S}$  is diagonal. This means that the noises in various detectors are uncorrelated. This is the case when the detectors of the network are in widely separated locations like for example the two LIGO detectors. However, this assumption is not always satisfied. An important case is the LISA detector where the noises of the three independent responses are correlated. Nevertheless for the case of LISA one can find a set of three combinations for which the noises are uncorrelated [69, 62] . When the cross spectrum matrix is diagonal the optimum $\mathcal{ℱ}$  -statistic is just the sum of $\mathcal{ℱ}$  -statistics in each detector.

A derivation of the likelihood function for an arbitrary network of detectors can be found in [32] , and applications of optimal filtering for the special cases of observations of coalescing binaries by networks of ground-based detectors are given in [40, 27, 67] and for the case of stellar mass binaries observed by LISA space-borne detector in [51] . A least square fit solution for the estimation of the sky location of a source of gravitational waves by a network of detectors for the case of a broad band burst was obtained in [36] .

There is also another important method for analyzing the data from a network of detectors – the search for coincidences of events among detectors. This analysis is particularly important when we search for supernova bursts the waveforms of which are not very well known. Such signals can be easily mimicked by non-Gaussian behavior of the detector noise. The idea is to filter the data optimally in each of the detector and obtain candidate events. Then one compares parameters of candidate events, like for example times of arrivals of the bursts, among the detectors in the network. This method is widely used in the search for supernovae by networks of bar detectors [11] .

4.11 Non-stationary, non-Gaussian, and non-linear data

Equations ( 32 ) and ( 33 ) provide maximum likelihood estimators only when the noise in which the signal is buried is Gaussian. There are general theorems in statistics indicating that the Gaussian noise is ubiquitous. One is the central limit theorem which states that the mean of any set of variates with any distribution having a finite mean and variance tends to the normal distribution.

The other comes from the information theory and says that the probability distribution of a random variable with a given mean and variance which has the maximum entropy (minimum information) is the Gaussian distribution. Nevertheless, analysis of the data from gravitational-wave detectors shows that the noise in the detector may be non-Gaussian (see, e.g., Figure 6 in [10] ). The noise in the detector may also be a non-linear and a non-stationary random process.

The maximum likelihood method does not require that the noise in the detector is Gaussian or stationary. However, in order to derive the optimum statistic and calculate the Fisher matrix we need to know the statistical properties of the data. The probability distribution of the data may be complicated, and the derivation of the optimum statistic, the calculation of the Fisher matrix components and the false alarm probabilities may be impractical. There is however one important result that we have already mentioned. The matched-filter which is optimal for the Gaussian case is also a linear filter that gives maximum signal-to-noise ratio no matter what is the distribution of the data. Monte Carlo simulations performed by Finn [32] for the case of a network of detectors indicate that the performance of matched-filtering (i.e., the maximum likelihood method for Gaussian noise) is satisfactory for the case of non-Gaussian and stationary noise.

In the remaining part of this section we review some statistical tests and methods to detect non-Gaussianity, non-stationarity, and non-linearity in the data. A classical test for a sequence of data to be Gaussian is the Kolmogorov–Smirnov test [23] . It calculates the maximum distance between the cumulative distribution of the data and that of a normal distribution, and assesses the significance of the distance. A similar test is the Lillifors test [23] , but it adjusts for the fact that the parameters of the normal distribution are estimated from the data rather than specified in advance. Another test is the Jarque–Bera test [43] which determines whether sample skewness and kurtosis are unusually different from their Gaussian values.

Let $x_k$  and $u_l$  be two discrete in time random processes ( $-\infty <k,l<\infty$  ) and let $u_l$  be independent and identically distributed (i.i.d.). We call the process $x_k$  linear if it can be represented by

$$\begin{array}{c}{x}_k={\sum }_{l=0}^N{a}_l{u}_{k-l},\end{array}$$

(88)

where $a_l$  are constant coefficients. If $u_l$  is Gaussian (non-Gaussian), we say that $x_l$  is linear Gaussian (non-Gaussian). In order to test for linearity and Gaussianity we examine the third-order cumulants of the data. The third-order cumulant $C_{kl}$  of a zero mean stationary process is defined by

$$\begin{array}{c}{C}_{kl}:=E\left[x_m{x}_{m+k}{x}_{m+l}\right].\end{array}$$

(89)

The bispectrum $S_2(f_1,f_2)$  is the two-dimensional Fourier transform of $C_{kl}$  . The bicoherence is defined as

$$\begin{array}{c}B(f_1,f_2):=\frac{S_2(f_1,f_2)}{S(f_1+f_2)S\left(f_1\right)S\left(f_2\right)},\end{array}$$

(90)

where $S\left(f\right)$  is the spectral density of the process $x_k$  . If the process is Gaussian then its bispectrum and consequently its bicoherence is zero. One can easily show that if the process is linear then its bicoherence is constant. Thus if the bispectrum is not zero, then the process is non-Gaussian; if the bicoherence is not constant then the process is also non-linear. Consequently we have the following hypothesis testing problems:

• 1. $H_1$  : The bispectrum of $x_k$  is nonzero.
• 2. $H_0$  : The bispectrum of $x_k$  is zero.

If Hypothesis  1 holds, we can test for linearity, that is, we have a second hypothesis testing problem:

• 3. $H_1^{\prime }$  : The bicoherence of $x_k$  is not constant.
• 4. $H_1^{\prime \prime }$  : The bicoherence of $x_k$  is a constant.

If Hypothesis  4 holds, the process is linear.

Using the above tests we can detect non-Gaussianity and, if the process is non-Gaussian, non-linearity of the process. The distribution of the test statistic $B(f_1,f_2)$  , Equation ( 90 ), can be calculated in terms of ${\chi }^2$  distributions. For more details see [38] .

It is not difficult to examine non-stationarity of the data. One can divide the data into short segments and for each segment calculate the mean, standard deviation and estimate the spectrum.

One can then investigate the variation of these quantities from one segment of the data to the other.

This simple analysis can be useful in identifying and eliminating bad data. Another quantity to examine is the autocorrelation function of the data. For a stationary process the autocorrelation function should decay to zero. A test to detect certain non-stationarities used for analysis of econometric time series is the Dickey–Fuller test [21] . It models the data by an autoregressive process and it tests whether values of the parameters of the process deviate from those allowed by a stationary model. A robust test for detection non-stationarity in data from gravitational-wave detectors has been developed by Mohanty [60] . The test involves applying Student's t-test to Fourier coefficients of segments of the data.

5 Acknowledgments

One of us (AK) acknowledges support from the National Research Council under the Resident Research Associateship program at the Jet Propulsion Laboratory, California Institute of Technology and from Max-Planck-Institut für Gravitationsphysik. We would like to thank Dr. Michele Vallisneri for preparing the figure. This research was in part funded by the Polish KBN Grant No. 1 P03B 029 27.

References

1. Abbott, B. et al. (LIGO Scientific Collaboration), “Analysis of LIGO data for gravitational waves from binary neutron stars”, Phys. Rev. D, 69, 122001–1–16, (2004). Related online version (cited on 8 January 2005): . ☻ open access ✓
2. Abbott, B. et al. (LIGO Scientific Collaboration), “First upper limits from LIGO on gravitational wave bursts”, Phys. Rev. D, 69, 102001–1–21, (2004). Related online version (cited on 8 January 2005): . ☻ open access ✓
3. Abbott, B. et al. (LIGO Scientific Collaboration), “Setting upper limits on the strength of periodic gravitational waves from PSR J1939+2134 using the first science data from the GEO 600 and LIGO detectors”, Phys. Rev. D, 69, 082004–1–16, (2004). Related online version (cited on 8 January 2005): . ☻ open access ✓
4. Adler, R.J., The Geometry of Random Fields, (Wiley, Chichester, U.K.; New York, U.S.A., 1981).
5. Allen, B., “The stochastic gravity-wave background: Sources and detection”, in Marck, J.-A., and Lasota, J.-P., eds., Relativistic gravitation and gravitational radiation, Proceedings of the Les Houches School of Physics, held in Les Houches, Haute Savoie, 26 September 6 October, 1995, Cambridge Contemporary Astrophysics, (Cambridge University Press, Cambridge, U.K., 1997). Related online version (cited on 8 January 2005): . ☻ open access ✓
6. Allen, B., Blackburn, J.K., Brady, P.R., Creighton, J.D., Creighton, T., Droz, S., Gillespie, A.D., Hughes, S.A., Kawamura, S., Lyons, T.T., Mason, J.E., Owen, B.J., Raab, F.J., Regehr, M.W., Sathyaprakash, B.S., Savage Jr, R.L., Whitcomb, S., and Wiseman, A.G., “Observational Limit on Gravitational Waves from Binary Neutron Stars in the Galaxy”, Phys. Rev. Lett., 83, 1498–1501, (1999). Related online version (cited on 8 January 2005): . ☻ open access ✓
7. Apostolatos, T.A., “Search templates for gravitational waves from precessing, inspiraling binaries”, Phys. Rev. D, 52, 605–620, (1995).
8. Armstrong, J.W., Estabrook, F.B., and Tinto, M., “Time-delay interferometry for space-based gravitational wave searches”, Astrophys. J., 527, 814–826, (1999).
9. Arnaud, N., Barsuglia, M., Bizouard, M., Brisson, V., Cavalier, F., Davier, M., Hello, P., Kreckelbergh, S., and Porter, E.K., “Coincidence and coherent data analysis methods for gravitational wave bursts in a network of interferometric detectors”, Phys. Rev. D, 68, 102001–1–18, (2003). Related online version (cited on 8 January 2005): . ☻ open access ✓
10. Astone, P., “Long-term operation of the Rome “Explorer” cryogenic gravitational wave detector”, Phys. Rev. D, 47, 362–375, (1993).
11. Astone, P., Babusci, D., Baggio, L., Bassan, M., Blair, D.G., Bonaldi, M., Bonifazi, P., Busby, D., Carelli, P., Cerdonio, M., Coccia, E., Conti, L., Cosmelli, C., D'Antonio, S., Fafone, V., Falferi, P., Fortini, P., Frasca, S., Giordano, G., Hamilton, W.O., Heng, I.S., Ivanov, E.N., Johnson, W.W., Marini, A., Mauceli, E., McHugh, M.P., Mezzena, R., Minenkov, Y., Modena, I., Modestino, G., Moleti, A., Ortolan, A., Pallottino, G.V., Pizzella, G., Prodi, G.A., Quintieri, L., Rocchi, A., Rocco, E., Ronga, F., Salemi, F., Santostasi, G., Taffarello, L., Terenzi, R., Tobar, M.E., Torrioli, G., Vedovato, G., Vinante, A., Visco, M., Vitale, S., and Zendri, J.P., “Methods and results of the IGEC search for burst gravitational waves in the years 1997–2000”, Phys. Rev. D, 68, 022001–1–33, (2003).
12. Astone, P., Babusci, D., Bassan, M., Borkowski, K.M., Coccia, E., D'Antonio, S., Fafone, V., Giordano, G., Jaranowski, P., Królak, A., Marini, A., Minenkov, Y., Modena, I., Modestino, G., Moleti, A., Pallottino, G.V., Pietka, M., Pizzella, G., Quintieri, L., Rocchi, A., Ronga, F., Terenzi, R., and Visco, M., “All-sky upper limit for gravitational radiation from spinning neutron stars”, Class. Quantum Grav., 20, S665–S676, (2003). Paper from the 7th Gravitational Wave Data Analysis Workshop, Kyoto, Japan, 17–19 December 2002.
13. Astone, P., Borkowski, K.M., Jaranowski, P., and Królak, A., “Data Analysis of gravitational-wave signals from spinning neutron stars. IV. An all-sky search”, Phys. Rev. D, 65, 042003–1–18, (2002).
14. Astone, P., Lobo, J.A., and Schutz, B.F., “Coincidence experiments between interferometric and resonant bar detectors of gravitational waves”, Class. Quantum Grav., 11, 2093–2112, (1994).
15. Balasubramanian, R., and Dhurandhar, S.V., “Estimation of parameters of gravitational waves from coalescing binaries”, Phys. Rev. D, 57, 3408–3422, (1998).
16. Balasubramanian, R., Sathyaprakash, B.S., and Dhurandhar, S.V., “Gravitational waves from coalescing binaries: detection strategies and Monte Carlo estimation of parameters”, Phys. Rev. D, 53, 3033–3055, (1996). Related online version (cited on 8 January 2005): . Erratum: Phys. Rev. D 54 (1996) 1860. ☻ open access ✓
17. Bayes, T., “An essay towards solving a problem in doctrine of chances”, Philos. Trans. R. Soc. London, 53, 293–315, (1763).
18. Blanchet, L., “Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries”, Living Rev. Relativity, 5, lrr-2002-3, (2002). URL (cited on 8 January 2005): . ☻ open access ✓
19. Bonazzola, S., and Gourgoulhon, E., “Gravitational waves from pulsars: Emission by the magnetic field induced distortion”, Astron. Astrophys., 312, 675–690, (1996). Related online version (cited on 8 January 2005): . ☻ open access ✓
20. Brady, P.R., Creighton, T., Cutler, C., and Schutz, B.F., “Searching for periodic sources with LIGO”, Phys. Rev. D, 57, 2101–2116, (1998). Related online version (cited on 8 January 2005): . ☻ open access ✓
21. Brooks, C., Introductory econometrics for finance, (Cambridge University Press, Cambridge, U.K.; New York, U.S.A., 2002).
22. Buonanno, A., Chen, Y., and Vallisneri, M., “Detection template families for gravitational waves from the final stages of binary-black-hole inspirals: Nonspinning case”, Phys. Rev. D, 67, 024016–1–50, (2003). Related online version (cited on 8 January 2005): . ☻ open access ✓
23. Conover, W.J., Practical Nonparametric Statistic, (Wiley, New York, U.S.A., 1980), 2nd edition.
24. Conway, J.H., and Sloane, N.J.A., Sphere Packings, Lattices and Groups, vol. 290 of Grundlehren der mathematischen Wissenschaften, (Springer, New York, U.S.A., 1999), 3rd edition.
25. Croce, R.P., Demma, T., Longo, M., Marano, S., Matta, V., Pierro, V., and Pinto, I.M., “Correlator bank detection of gravitational wave chirps—False-alarm probability, template density, and thresholds: Behind and beyond the minimal-match issue”, Phys. Rev. D, 70, 122001–1–19, (2004). Related online version (cited on 8 January 2005): . ☻ open access ✓
26. Cutler, C., Apostolatos, T.A., Bildsten, L., Finn, L.S., Flanagan, É.É., Kennefick, D., Markovic, D.M., Ori, A., Poisson, E., Sussman, G.J., and Thorne, K.S., “The Last Three Minutes: Issues in Gravitational Wave Measurements of Coalescing Compact Binaries”, Phys. Rev. Lett., 70, 2984–2987, (1993). Related online version (cited on 8 January 2005): . ☻ open access ✓
27. Cutler, C., and Flanagan, É.É., “Gravitational Waves from Merging Compact Binaries: How Accurately Can One Extract the Binary's Parameters from the Inspiral Waveform?”, Phys. Rev. D, 49, 2658–2697, (1994). Related online version (cited on 8 January 2005): . ☻ open access ✓
28. Davis, M.H.A., “A Review of Statistical Theory of Signal Detection”, in Schutz, B.F., ed., Gravitational Wave Data Analysis, Proceedings of the NATO Advanced Research Workshop, held at Dyffryn House, St. Nichols, Cardiff, Wales, 6–9 July 1987, vol. 253 of NATO ASI Series C, 73–94, (Kluwer, Dordrecht, Netherlands; Boston, U.S.A., 1989).
29. Dhurandhar, S.V., and Sathyaprakash, B.S., “Choice of filters for the detection of gravitational waves from coalescing binaries. II. Detection in colored noise”, Phys. Rev. D, 49, 1707–1722, (1994).
30. Dhurandhar, S.V., and Schutz, B.F., “Filtering coalescing binary signals: Issues concerning narrow banding, thresholds, and optimal sampling”, Phys. Rev. D, 50, 2390–2405, (1994).
31. Estabrook, F.B., and Wahlquist, H.D., “Response of Doppler spacecraft tracking to gravitational radiation”, Gen. Relativ. Gravit., 439, 439–447, (1975).
32. Finn, L.S., “Aperture synthesis for gravitational-wave data analysis: Deterministic Sources”, Phys. Rev. D, 63, 102001–1–18, (2001). Related online version (cited on 8 January 2005): . ☻ open access ✓
33. Finn, L.S., and Chernoff, D.F., “Observing binary inspiral in gravitational radiation: One interferometer”, Phys. Rev. D, 47, 2198–2219, (1993). Related online version (cited on 8 January 2005): . ☻ open access ✓
34. Fisz, M., Probability Theory and Mathematical Statistics, (Wiley, New York, U.S.A., 1963).
35. Giampieri, G., “On the antenna pattern of an orbiting interferometer”, Mon. Not. R. Astron. Soc., 289, 185–195, (1997).
36. Gürsel, Y., and Tinto, M., “Nearly optimal solution to the inverse problem for gravitational-wave bursts”, Phys. Rev. D, 40, 3884–3938, (1989).
37. Helström, C.W., Statistical Theory of Signal Detection, vol. 9 of International Series of Monographs in Electronics and Instrumentation, (Pergamon Press, Oxford, U.K., New York, U.S.A., 1968), 2nd edition.
38. Hinich, M.J., “Testing for Gaussianity and linearity of a stationary time series”, J. Time Series Anal., 3, 169–176, (1982).
39. Hough, J., and Rowan, S., “Gravitational Wave Detection by Interferometry (Ground and Space)”, Living Rev. Relativity, 3, lrr-2000-3, (2000). URL (cited on 8 January 2005): . ☻ open access ✓
40. Jaranowski, P., and Królak, A., “Optimal solution to the inverse problem for the gravitational wave signal of a coalescing binary”, Phys. Rev. D, 49, 1723–1739, (1994).
41. Jaranowski, P., and Królak, A., “Data Analysis of gravitational-wave signals from spinning neutron stars. III. Detection statistics and computational requirements”, Phys. Rev. D, 61, 062001–1–32, (2000).
42. Jaranowski, P., Królak, A., and Schutz, B.F., “Data Analysis of gravitational-wave signals from spinning neutron stars: The signal and its detection”, Phys. Rev. D, 58, 063001–1–24, (1998).
43. Judge, G.G., Hill, R.C., Griffiths, W.E., Lutkepohl, H., and Lee, T.-C., The Theory and Practice of Econometrics, (Wiley, New York, U.S.A., 1980).
44. Kafka, P., “Optimal Detection of Signals through Linear Devices with Thermal Noise Sources and Application to the Munich-Frascati Weber-Type Gravitational Wave Detectors”, in De Sabbata, V., and Weber, J., eds., Topics in Theoretical and Experimental Gravitation Physics, Proceedings of the International School of Cosmology and Gravitation held in Erice, Trapani, Sicily, March 13–25, 1975, vol. 27 of NATO ASI Series B, 161, (Plenum Press, New York, U.S.A., 1977).
45. Kendall, M., and Stuart, A., The Advanced Theory of Statistics. Vol. 2: Inference and Relationship, number 2, (C. Griffin, London, 1979).
46. Kokkotas, K.D., Królak, A., and Tsegas, G., “Statistical analysis of the estimators of the parameters of the gravitational-wave signal from a coalescing binary”, Class. Quantum Grav., 11, 1901–1918, (1994).
47. Kotelnikov, V.A., The theory of optimum noise immunity, (McGraw-Hill, New York, U.S.A., 1959).
48. Królak, A., Kokkotas, D., and Schäfer, G., “On estimation of the post-Newtonian parameters in the gravitational-wave emission of a coalescing binary”, Phys. Rev. D, 52, 2089–2111, (1995). Related online version (cited on 8 January 2005): . ☻ open access ✓
49. Królak, A., Lobo, J.A., and Meers, B.J., “Estimation of the parameters of the gravitational-wave signal of a coalescing binary system”, Phys. Rev. D, 48, 3451–3462, (1993).
50. Królak, A., and Schutz, B.F., “Coalescing binaries – probe to the Universe”, Gen. Relativ. Gravit., 19, 1163–1171, (1987).
51. Królak, A., Tinto, M., and Vallisneri, M., “Optimal filtering of the LISA data”, Phys. Rev. D, 70, 022003–1–24, (2004). Related online version (cited on 8 January 2005): . ☻ open access ✓
52. Lagarias, J.C., Reeds, J.A., Wright, M.H., and Wright, P.E., “Convergence properties of the Nelder–Mead simplex method in low dimensions”, SIAM J. Optimiz., 9, 112–147, (1998).
53. Lehmann, E.L., Testing Statistical Hypothesis, (Wiley, New York, U.S.A., 1959).
54. Lehmann, E.L., Theory of Point Estimation, (Wiley, New York, U.S.A., 1983).
55. Liptser, R.S., and Shiryaev, A.N., Statistics of Random Processes, 2 vols., Applications of Mathematics, (Springer, New York, U.S.A., 1977).
56. LISA: Pre-phase A report, December 1998, MPQ 223, (Max-Planck-Institut für Quantenoptik, Garching, Germany, 1998).
57. McDonough, R.N., and Whalen, A.D., Detection of signals in noise, (Academic Press, San Diego, U.S.A., 1995), 2nd edition.
58. Meyer, C., Matrix Analysis and Applied Linear Algebra, (SIAM, Philadelphia, U.S.A., 2000).
59. Mohanty, S.D., “Hierarchical search strategy for the detection of gravitational waves from coalescing binaries: Extension to post-Newtonian waveforms”, Phys. Rev. D, 57, 630–658, (1998). Related online version (cited on 8 January 2005): . ☻ open access ✓
60. Mohanty, S.D., “A robust test for detecting non-stationarity in data from gravitational wave detectors”, Phys. Rev. D, 61, 122002–1–12, (2000). Related online version (cited on 8 January 2005): . ☻ open access ✓
61. Mohanty, S.D., and Dhurandhar, S.V., “Hierarchical search strategy for the detection of gravitational waves from coalescing binaries”, Phys. Rev. D, 54, 7108–7128, (1996).
62. Nayak, K.R., Pai, A., Dhurandhar, S.V., and Vinet, J.-Y., “Improving the sensitivity of LISA”, Class. Quantum Grav., 20, 1217–1232, (2003).
63. Nicholson, D., Dickson, C.A., Watkins, W.J., Schutz, B.F., Shuttleworth, J., Jones, G.S., Robertson, D.I., MacKenzie, N.L., Strain, K.A., Meers, B.J., Newton, G.P., Ward, H., Cantley, C.A., Robertson, N.A., Hough, J., Danzmann, K., Niebauer, T.M., Rüdiger, A., Schilling, R., Schnupp, L., and Winkler, W., “Results of the first coincident observations by two laser-interferometric gravitational wave detectors”, Phys. Lett. A, 218, 175–180, (1996).
64. Nicholson, D., and Vecchio, A., “Bayesian bounds on parameter estimation accuracy for compact coalescing binary gravitational wave signals”, Phys. Rev. D, 57, 4588–4599, (1998). Related online version (cited on 8 January 2005): . ☻ open access ✓
65. Niebauer, T.M., Rüdiger, A., Schilling, R., Schnupp, L., Winkler, W., and Danzmann, K., “Pulsar search using data compression with the Garching gravitational wave detector”, Phys. Rev. D, 47, 3106–3123, (1993).
66. Owen, B.J., “Search templates for gravitational waves from inspiraling binaries: Choice of template spacing”, Phys. Rev. D, 53, 6749–6761, (1996). Related online version (cited on 8 January 2005): . ☻ open access ✓
67. Pai, A., Dhurandhar, S.V., and Bose, S., “Data-analysis strategy for detecting gravitational-wave signals from inspiraling compact binaries with a network of laser-interferometric detectors”, Phys. Rev. D, 64, 042004–1–30, (2001). Related online version (cited on 8 January 2005): . ☻ open access ✓
68. Poor, H.V., An Introduction to Signal Detection and Estimation, (Springer, New York, U.S.A., 1994), 2nd edition.
69. Prince, T.A., Tinto, M., Larson, S.L., and Armstrong, J.W., “The LISA optimal sensitivity”, Phys. Rev. D, 66, 122002–1–7, (2002).
70. Rife, D.C., and Boorstyn, R.R., “Single tone parameter estimation from discrete-time observations”, IEEE Trans. Inform. Theory, 20, 591–598, (1974).
71. Rubbo, L.J., Cornish, N.J., and Poujade, O., “Forward modeling of space-borne gravitational wave detectors”, Phys. Rev. D, 69, 082003–1–14, (2003). Related online version (cited on 8 January 2005): . ☻ open access ✓
72. Sathyaprakash, B.S., and Dhurandhar, S.V., “Choice of filters for the detection of gravitational waves from coalescing binaries”, Phys. Rev. D, 44, 3819–3834, (1991).
73. Schutz, B.F., “Determining the nature of the Hubble constant”, Nature, 323, 310–311, (1986).
74. Schutz, B.F., ed., Gravitational Wave Data Analysis, Proceedings of the NATO Advanced Research Workshop held at Dyffryn House, St. Nichols, Cardiff, Wales, 6–9 July 1987, vol. 253 of NATO ASI Series C, (Kluwer, Dordrecht, Netherlands; Boston, U.S.A., 1989).
75. Schutz, B.F., “Data processing, analysis and storage for interferometric antennas”, in Blair, D.G., ed., The Detection of Gravitational Waves, 406–452, (Cambridge University Press, Cambridge, U.K.; New York, U.S.A., 1991).
76. Sengupta, S.A., Dhurandhar, S.V., and Lazzarini, A., “Faster implementation of the hierarchical search algorithm for detection of gravitational waves from inspiraling compact binaries”, Phys. Rev. D, 67, 082004–1–14, (2003).
77. Tagoshi, H., Kanda, N., Tanaka, T., Tatsumi, D., Telada, S., Ando, M., Arai, K., Araya, A., Asada, H., Barton, M.A., Fujimoto, M.-K., Fukushima, M., Futamase, T., Heinzel, G., Horikoshi, G., Ishizuka, H., Kamikubota, N., Kawabe, K., Kawamura, S., Kawashima, N., Kojima, Y., Kozai, Y., Kuroda, K., Matsuda, N., Matsumura, S., Miki, S., Mio, N., Miyakawa, O., Miyama, S.M., Miyoki, S., Mizuno, E., Moriwaki, S., Musha, M., Nagano, S., Nakagawa, K., Nakamura, T., Nakao, K., Numata, K., Ogawa, Y., Ohashi, M., Ohishi, N., Okutomi, A., Oohara, K., Otsuka, S., Saito, Y., Sasaki, M., Sato, S., Sekiya, A., Shibata, M., Shirakata, K., Somiya, K., Suzuki, T., Takahashi, R., Takamori, A., Taniguchi, S., Tochikubo, K., Tomaru, T., Tsubono, K., Tsuda, N., Uchiyama, T., Ueda, A., Ueda, K., Waseda, K., Watanabe, Y., Yakura, H., Yamamoto, K., and Yamazaki, T. (The TAMA Collaboration), “First search for gravitational waves from inspiraling compact binaries using TAMA300 data”, Phys. Rev. D, 63, 062001–1–5, (2001). Related online version (cited on 8 January 2005): . ☻ open access ✓
78. Tanaka, T., and Tagoshi, H., “Use of new coordinates for the template space in hierarchical search for gravitational waves from inspiraling binaries”, Phys. Rev. D, 62, 082001–1–8, (2000). Related online version (cited on 8 January 2005): . ☻ open access ✓
79. Thorne, K.S., “Gravitational radiation”, in Hawking, S.W., and Israel, W., eds., Three Hundred Years of Gravitation, 330–458, (Cambridge University Press, Cambridge, U.K.; New York, U.S.A., 1987).
80. Tinto, M., and Armstrong, J.W., “Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation”, Phys. Rev. D, 59, 102003–1–11, (1999).
81. Table of Q Functions, RAND Research Memorandum, M-339, (U.S. Air Force, Rand Corporation, Santa Monica, U.S.A., 1950).
82. Van Trees, H.L., Detection, Estimation and Modulation Theory. Part 1: Detection, Estimation, and Linear Modulation Theory, number 1, (Wiley, New York, U.S.A., 1968).
83. Wainstein, L.A., and Zubakov, V.D., Extraction of signals from noise, (Prentice-Hall, Englewood Cliffs, U.S.A., 1962).
84. Weber, J., “Evidence for Discovery of Gravitational Radiation”, Phys. Rev. Lett., 22, 1320–1324, (1969).
85. Wong, E., Introduction to Random Processes, (Springer, New York, U.S.A., 1983).
86. Wong, E., and Hajek, B., Stochastic Processes in Engineering Systems, (Springer, New York, U.S.A., 1985).
87. Woodward, P.M., Probability and information theory with applications to radar, (Pergamon Press, London, U.K., 1953).
88. Zieliński, R., “Theory of parameter estimation”, in Królak, A., ed., Mathematics of Gravitation. Part II: Gravitational Wave Detection, Proceedings of the Workshop on Mathematical Aspects of Theories of Gravitation, held in Warsaw, February 29–March 30, 1996, vol. 41(II) of Banach Center Publications, 209–220, (Institute of Mathematics, Polish Academy of Sciences, Warsaw, Poland, 1997).

Note: The reference version of this article is published by Living Reviews in Relativity