Vetoing techniques ($r^2$ time/frequency test)

The second technique vetoing or discrimination test operates in the frequency domain, and is described here. It is a very stringent test, which determines if the hypothetical chirp which has been found in the data stream is consistent with a true binary inspiral chirp summed with Gaussian interferometer noise. If this is true, it should be possible to subtract the (best fit) chirp from the signal, and be left with a signal stream that is consistent with Gaussian IFO noise. One of the nice features of this technique is that it can be statistically characterized in a rigorous way. We follow the same general course as in Section on Wiener filtering, first considering the case of a ``single phase" phase signal, then considering the case of an ``unknown phase" signal.

In the single-phase case, suppose that one of our optimal chirp filters $\tilde Q$ is triggered with a large SNR at time $t_0$. We suppose that the signal which was responsible for this trigger may be written in either the time or the frequency domain as

$$h(t) = C(t) + n(t) = \alpha T(t-t_0) + n(t)$$

$$\iff$$

$$\tilde h(f) = \alpha \tilde T(f) {\rm e}^{2 \pi i f t_0} + \tilde n.$$ (6.24.105)

We assume that we have identified what is believed to be the ``correct" template $T$, by the procedure already described of maximizing the SNR over arrival time and template, and have used this to estimate $\alpha$. We assume that $t_0$ has been determined exactly (a good approximation since it can be estimated to high accuracy). Our goal is to construct a statistic which will indicate if our estimate of $\alpha$ and identification of $T$ are credible.

We will denote the signal value at time offset $t_0$ by the real number $S$:

$$S = 2 \int_{-f_{\rm Ny}}^{f_{\rm Ny}} df \; { \tilde h(f) \t... ...^* (f) \over S_h(\vert f\vert)} \; {\rm e}^{- 2 \pi i f t_0}.$$ (6.24.106)

(Here, $f_{\rm Ny}$ denotes the Nyqist frequency, one-half of the sampling rate.) The expected value of $S$ is $\langle S \rangle = \alpha$, although of course since we are discussing a single instance, it's actual value will in general be different. The chirp template $T$ is normalized so that the expected value $\langle N^2 \rangle=1$:

$$4 \int_0^{f_{\rm Ny}} df \; { \vert \tilde T(f)\vert^2 \over S_h(\vert f\vert)} = 1.$$ (6.24.107)

We are going to investigate if this signal is ``really" due to a chirp by investigating the way in which $S$ gets its contribution from different ranges of frequencies. To do this, break up the integration region in this integral into a set of $p$ disjoint subintervals $\Delta f_1,\cdots,\Delta f_p$ whose union is the entire range of frequencies from DC to Nyquist. Here $p$ is a small integer (for example, $p=8$). This splitup can be performed using the GRASP function splitup(). The frequency intervals:

$$\Delta f_1 = \{ f \; \vert \; 0< f < f_1 \}$$

$$\Delta f_2 = \{ f \; \vert \; f_1 < f < f_2 \}$$

$$\cdots$$

$$\Delta f_p = \{ f \; \vert \; f_{p-1}< f < f_{\rm Ny} \},$$ (6.24.108)

are defined by the condition that the expected signal contributions in each frequency band from a chirp are equal:

$$4 \int_{\Delta f_i} df \; { \vert \tilde T(f)\vert^2 \over S_h(\vert f\vert)} = {1 \over p}$$ (6.24.109)

A typical set of frequency intervals in shown in Figure .

Figure: A typical set of frequency intervals $\Delta f_i$ for the case $p=4$.

Because the filter is optimal, this also means that the expected noise contributions in each band from the chirp is the same. The frequency subintervals $\Delta f_i$ are narrow in regions of frequency space where the interferometer is quiet, and are broad in regions where the IFO is noisy.

Now, define a set of $p$ signal values, one for each frequency interval:

$$S_i = 2 \int_{-\Delta f_i \cup \Delta f_i} df \; { \tilde h(... ... \; {\rm e}^{- 2 \pi i f t_0} \quad {\rm for\ } i=1,\cdots,p.$$ (6.24.110)

We have included both the positive and negative frequency subintervals to ensure that the $S_i$ are real. If the detector output is Gaussian noise plus a true chirp, the $S_i$ are $p$ normal random variables, with a mean value of $\langle S_i \rangle =\langle S \rangle /p$ and a variance determined by the expected value of the noise-squared:

$$\sigma = \langle S_i^2 \rangle - \langle S_i \rangle^2 = {1 \over p}.$$ (6.24.111)

>From these signal values we can construct a useful time/frequency statistic.

To characterize the statistic, we will need the probability distribution of the $S_i$. Because each of these values is a sum over different (non-overlapping) frequency bins, they are independent random Gaussian variables with unknown mean values. Their a-priori probability distribution is

$$P(S_1,\cdots,S_p) = \prod_{i=1}^p (2 \pi \sigma)^{-1/2} {\rm e}^{- { \left(S_i - \alpha /p \right )^2 / 2 \sigma}}$$ (6.24.112)

The statistic that we will construct addresses the question, ``are the measured values of $S_i$ consistent with the assumption that the measured signal is Gaussian detector noise plus $\alpha T$?" One small difficulty is that the value of $\alpha$ that appears in () is not known to us: we only have an estimate of its value. To overcome this, we first construct a set of values denoted

$$\Delta S_i \equiv S_i-S/p.$$ (6.24.113)

These are not independent normal random variables: they are correlated since $\sum_{i=1}^p \Delta S_i$ vanishes. To proceed, we need to calculate the probability distribution of the $\Delta S_i$, which we denote by $\bar P(\Delta S_1,\cdots,\Delta S_p)$. This quantity is defined by the relation that the integral of any function of $p$ variables $F(y_1,\cdots,y_p)$ with respect to the measure defined by this probability distribution satisfies

$$\int dy_1 \cdots dy_p \bar P (y_1,\cdots,y_p) F(y_1,\cdots,y_p)$$

$$=$$

$$\int dx_1 \cdots dx_p P (x_1,\cdots,x_p) F(x_1 - {1 \over p} \sum_{i=1}^p x_i,\cdots,x_p - {1 \over p} \sum_{i=1}^p x_i ).$$ (6.24.114)

[Note that in this expression and the following ones, all integrals are from $-\infty$ to $\infty$.] This may be used to find a closed form for $\bar P$: let $F(y_1,\cdots,y_p) = \delta(y_1 - \Delta S_1) \cdots \delta(y_p - \Delta S_p)$. This gives

$$\bar P (\Delta S_1,\cdots,\Delta S_p) = \int \prod_{i=1}^pdx... ...gma} \delta(x_i - \Delta S_i - {1 \over p} \sum_{j=1}^p x_j ).$$ (6.24.115)

To evaluate the integral, change variables from $(x_1,\cdots,x_p)$ to $(z_1,\cdots,z_{p-1},W)$ defined by

$$W = x_1 + \cdots + x_p$$

$$z_1 = x_1 - W/p$$

$$\cdots$$

$$z_{p-1} = x_{p-1} - W/p$$ (6.24.116)

which can be inverted to yield

$$x_1 = z_1 + W/p$$

$$\cdots$$

$$x_{p-1} = z_{p-1} + W/p$$ (6.24.117)

$$x_p = W/p - z_1 - \cdots - z_{p-1}$$

The Jacobian of this coordinate transformation is:

$$J= \det \left[ { \partial(x_1,\cdots,x_p ) \over \partial(z_... ...\cdots & 1 & 1/p \cr -1 &-1 & \cdots & -1 & 1/p \cr } \right]$$ (6.24.118)

Using the linearity in rows of the determinant, it is straightforward to show that $J=1$.

The integral may now be written as

$$\bar P (\Delta S_1,\cdots,\Delta S_p) = \int dz_1 \cdots dz_{p-1} dW (2 \pi \sigma)^{-p/2} {\rm e}^{-[(x_1 - \alpha/p)^2 + \cdots + (x_p - \alpha/p)^2]/ 2 \sigma}$$

$$\times \delta(z_1-\Delta S_1) \cdots \delta(z_{p-1}-\Delta S_{p-1}) \delta(z_1+\cdots + z_{p-1}+\Delta S_p).$$ (6.24.119)

A few moments of algebra shows that the exponent may be expressed in terms of the new integration variables as

$$(x_1 - \alpha/p)^2 + \cdots + (x_p - \alpha/p)^2 = z_1^2 + \cdots + z_{p-1}^2 + (W-\alpha)^2/p + (z_1 + \cdots + z_{p-1})^2$$ (6.24.120)

and thus the integral yields

$$\bar P (\Delta S_1,\cdots,\Delta S_p) = \int dW (2 \pi \sigma)^{-p/2} {\rm e}^{-[\Delta S_1^2 + \cdots + \Delta S_p^2 + (W-\alpha)^2/p]/2\sigma} \delta(\Delta S_1 + \cdots +\Delta S_p)$$

$$= (2 \pi \sigma)^{-p/2} (2 \pi \sigma p)^{1/2} {\rm e}^{-[\Delta S_1^2 + \cdots + \Delta S_p^2]/2 \sigma} \delta(\Delta S_1 + \cdots +\Delta S_p)$$ (6.24.121)

This probability distribution arises because we do not know the true mean value of $S$ which is $\alpha=\langle S \rangle$ but can only estimate it using the actual measured value of $S$. Similar problems arise whenever the mean of a distribution is not know but must be estimated (problem 14-7 of [24]). This probability distribution is ``as close as you can get to Gaussian" subject to the constraint that the sum of the $\Delta S_i$ must vanish. It is significant that this probability density function is completely independent of $\alpha$, which means that the properties of the $\Delta S_i$ do not depend upon whether a signal is present or not.

The individual $\Delta S_i$ have identical mean and variance, which may be easily calculated from the probability distribution function (). For example the mean is zero: $\langle \Delta S_i \rangle = 0$. To calculate the variance, let $F(y_1,\cdots,y_p) = y_1^2$ in (). One finds

$$\langle \Delta S_i^2 \rangle = {1 \over p}\left( 1 - {1 \over p} \right)$$ (6.24.122)

Now that we have calculated the probability distribution of the $\Delta S_i$ it is straightforward to construct and characterize a $\chi^2$-like statistic which we will call $r^2$.

Define the statistic

$$r^2 = \sum_{i=1}^p (\Delta S_i)^2.$$ (6.24.123)

>From () it is obvious that for Gaussian noise plus a chirp the statistical properties of $r^2$ are independent of $\alpha$: it has the same statistical properties if a chirp signal is present or not. For this reason, the value of $r^2$ provides a powerful method of testing the hypothesis that the detector's output is Gaussian noise plus a chirp. If the detector's output is of this form, then the value of $r^2$ is unlikely to be much larger than its expected value (this statement is quantified below). On the other hand, if the filter was triggered by a spurious noise event that does not have the correct time/frequency distribution, then $r^2$ will typically have a value that is very different than the value that it has for Gaussian noise alone (or equivalently, for Gaussian noise plus a chirp).

The expected value of $r^2$ is trivial to calculate

$$\langle r^2 \rangle = p \langle \Delta S_i^2 \rangle = 1 - {1 \over p}$$ (6.24.124)

One can also easily compute the probability distribution of $r$ using (). The probability that $r>R$ in the presence of a true chirp signal is the integral of () over the region $r>R$. In the $p$-dimensional space, the integrand vanishes except on a $p-1$-plane, where it is spherically-symmetric. To evaluate the integral, introduce a new set of orthonormal coordinates $(u_1,\cdots,u_p)$ obtained from any orthogonal transformation on $(\Delta S_1,\cdots,\Delta S_p)$ for which the new $p$'th coordinate is orthogonal to the hyperplane $\Delta S_1 + \cdots +\Delta S_p = 0$. Hence $u_p = p^{-1/2} \left[ \Delta S_1 + \cdots +\Delta S_p \right]$. Our statistic $r^2$ is also the squared radius $r^2 = u_1^2 + \cdots + u_p^2$ in these coordinates. Hence

$$P(r>R) = \int_{r^2 > R^2} \bar P$$

$$= (2 \pi \sigma)^{-p/2} (2 \pi \sigma p)^{1/2} \int_{r^2 > R^2} {\rm e}^{-r^2/2\sigma} \delta(\sqrt{p} u_p) du_1 \cdots du_p.$$ (6.24.125)

It's now easy to do the integral over the coordinate $u_p$, and having done this, we are left with a spherically-symmetric integral over $R^{p-1}$:

$$P(r>R) = (2 \pi \sigma)^{(1-p)/2} \int_{r^2 > R^2} {\rm e}^{-r^2/2\sigma} du_1 \cdots du_{p-1}$$

$$= \Omega_{p-2} \int_R^\infty r^{p-2} {\rm e}^{-r^2/2\sigma} dr$$

$$= {1 \over \Gamma({p-1 \over 2}) } \int_{R^2/2 \sigma}^\infty x^{(p-3)/2} {\rm e}^{-x} dx$$

$$= Q( {p-1 \over 2} ,{ R^2 \over 2 \sigma}),$$ (6.24.126)

where $\Omega_{n} = {2 \pi^{(n+1)/2} \over \Gamma({n+1 \over 2})}$ is the $n-$volume of a unit-radius $n-$sphere $S^n$. The incomplete gamma function $Q$ is the same function that describes the likelihood function in the traditional $\chi^2$ test [the Numerical Recipes function gammq(a,x)]. Figure show a graph of $P(r>R)$ for some different values of the parameter $p$. The appearance of $p-1$ in these expressions reflects the fact that although $r^2$ is a sum of the squares of $p$ Gaussian variables, these variables are subject to a single constraint (their sum vanishes) and thus the number of degrees of freedom is $p-1$.

In practice (based on CIT 40-meter data) breaking up the frequency range into $p=8$ intervals provides a very reliable veto for rejecting events that trigger an optimal filter, but which are not themselves chirps. The value of $Q(3.5,10.0) = 0.0056\cdots$ so if $r^2>2.5$ then one can conclude that the likelihood that a given trigger is actually due to a chirp is less than $0.6\%$; rejecting or vetoing such events will only reduce the ``true event" rate by $0.6\%$. However in practice it eliminates almost all other events that trigger an optimal filter; a noisy event that stimulates a binary chirp filter typically has $r^2 \approx 100$ or larger!

The previous analysis for the ``single-phase" case assumes that we have found the correct template $T$ describing the signal. In searching for a binary inspiral chirp however, the signal is a linear combination of the two different possible phases:

$$h(t) = C(t) + n(t) = \alpha T_0(t-t_0) + \beta T_{90}(t-t_0) + n(t)$$

$$\iff$$

$$\tilde h(f) = \left[ \alpha \tilde T_0(f) + \beta \tilde T_{90}(f) \right] {\rm e}^{2 \pi i f t_0} + \tilde n.$$ (6.24.127)

and the amplitudes $\alpha$ and $\beta$ are unknown. The reader might well wonder why we can't simply construct a single properly normalized template as

$$T = \left( {\alpha \over \sqrt{\alpha^2 + \beta^2}} \right) ... ... \left( {\beta \over \sqrt{\alpha^2 + \beta^2}} \right) T_{90}$$ (6.24.128)

and then use the previously-described ``single phase" method. In principle, this would work properly. The problem is that we do not know the correct values of $\alpha$ and $\beta$ . Since $\alpha = {\rm Re} \> \langle S(t_0) \rangle$ and $\beta = {\rm Im} \> \langle S(t_0) \rangle$, we can estimate the values of $\alpha$ and $\beta$ from the real and imaginary parts of the measured signal, however these estimates will not give the true values. For this reason, an $r^2$ statistic and test can be constructed for the ``two-phase" case, but it has twice the number of degrees of freedom as the ``single-phase" case.

The description and characterization of the $r^2$ test for the two phase case is similar to the single-phase case. For the two phase case, the signal is a complex number

$$S = 2 \int_{-f_{\rm Ny}}^{f_{\rm Ny}} df \; { \tilde h(f) \l... ...right] \over S_h(\vert f\vert)} \; {\rm e}^{- 2 \pi i f t_0}.$$ (6.24.129)

The templates for the individual phases are normalized as before:

$$4 \int_0^{f_{\rm Ny}} df \; { \vert \tilde T_0 (f)\vert^2 \o... ...\tilde T_0(f) \tilde T^*_{90}(f) \over S_h(\vert f\vert)} = 0.$$ (6.24.130)

This assume the same adiabatic limit discussed earlier: ${\dot f}/f << f$. In this limit, the frequency intervals $\Delta f_i$ are identical for either template. We define signal values in each frequency band in the same way as before, except now these are complex:

$$S_i = 2 \int_{-\Delta f_i \cup \Delta f_i} df \; { \tilde h(... ... \; {\rm e}^{- 2 \pi i f t_0} \quad {\rm for\ } i=1,\cdots,p.$$ (6.24.131)

The mean value of the signal in each frequency band is

$$\langle S_i \rangle = \langle S \rangle /p = (\alpha + i \beta)/p,$$ (6.24.132)

and the variance of either the real or imaginary part is $\sigma=1/p$ as before, so that the total variance is twice as large as in the single phase case:

$$\langle \vert S_i \vert ^2 \rangle - \vert \langle S_i \rangle \vert^2 = {2 \over p}.$$ (6.24.133)

The signal values are now characterized by the probability distribution

$$P(S_1,\cdots,S_p) = \prod_{i=1}^p (2 \pi \sigma)^{-1} {\rm ... ...t\vert S_i - \alpha /p - i \beta/p \right\vert^2 / 2 \sigma}}.$$ (6.24.134)

Note that the arguments of this function are complex; for this reason the overall normalization factors have changed from the single-phase case. We now construct complex quantities which are the difference between the actual signal measured in a frequency band and the expected value for our templates and phases:

$$\Delta S_i \equiv S_i-S/p.$$ (6.24.135)

The probability distribution of these differences is still defined by () but in that expression, the variables of integration $x_i$ and $y_i$ are integrated over the complex plane (real and imaginary parts from $-\infty$ to $\infty$), and $F$ is any function of $p$ complex variables. As before, we can calculate $\bar P$ by choosing $F$ correctly, in this case as $F(y_1,\cdots,y_p) = \delta^2(y_1 - \Delta S_1) \cdots \delta^2(y_p - \Delta S_p)$, where $\delta^2(z) \equiv \delta({\rm Re} \; z) \delta({\rm Im} \; z)$. The same procedure as before then yields the probability distribution function

$$\bar P (\Delta S_1,\cdots,\Delta S_p) = (2 \pi \sigma)^{-p} ... ...2 \right] /2 \sigma} \delta^2(\Delta S_1 + \cdots +\Delta S_p)$$ (6.24.136)

It is now easy to see that the expectation of the signal differences is still zero $\langle \Delta S_i \rangle = 0$ but the variances are twice as large as in the single-phase case:

$$\langle \vert\Delta S_i\vert^2 \rangle = {2 \over p}\left( 1 - {1 \over p} \right).$$ (6.24.137)

The $r^2$ statistic is now defined by

$$r^2 = \sum_{i=1}^p \left\vert \Delta S_i \right\vert ^2.$$ (6.24.138)

and has an expectation value which as twice as large as in the single-phase case:

$$\langle r^2 \rangle = p \langle \left\vert \Delta S_i \right\vert^2 \rangle = 2 - {2 \over p}.$$ (6.24.139)

The calculation of the distribution function of $r^2$ is similar to the single phase case (but with twice the number of degrees of freedom) and gives the incomplete $\Gamma$-function

$$P(r>R) = (2 \pi \sigma)^{-p} (2 \pi \sigma p) \int_{r^2 > R^2} {\rm e}^{-r^2/2\sigma} \delta^2(\sqrt{p} u_p) du_1 \cdots du_p$$

$$= (2 \pi \sigma)^{(1-p)} \int_{r^2 > R^2} {\rm e}^{-r^2/2\sigma} du_1 \cdots du_{p-1}$$

$$= P(r>R) = Q( { p-1} ,{ R^2 \over 2 \sigma}) = Q( { p-1} ,{ p R^2 \over 2})$$ (6.24.140)

This is precisely the distribution of a $\chi^2$ statistic with $2p-2$ degrees of freedom: each of the $p$ variables $\Delta S_i$ has 2 degrees of freedom, and there are two constraints since the sum of both the real and imaginary parts vanishes. In fact since the expectation value of the $\chi^2$ statistic is just the number of degrees of freedom:

$$\langle \chi^2 \rangle = 2 p - 2$$ (6.24.141)

the relationship between the $r^2$ and $\chi^2$ statistic may be obtained by comparing equations () and (), giving

$$\chi^2 = p r^2.$$ (6.24.142)

Figure: The probability that the $r^2$ statistic exceeds a given threshold $R^2$ is shown for both the single-phase and two-phase test, for $p=8$ and $p=16$ frequency ranges. For example, for the single-phase $p=8$ test, the probability that $r^2 > 2.31$ is 1% for a chirp plus Gaussian noise. For the single-phase test with $p=16$ the probability of exceeding the same threshold is about $10^{-3}$.