next up previous contents
Next: How does the test Up: GRASP Routines: Gravitational Radiation Previous: Vetoing techniques (time domain   Contents


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

$\displaystyle h(t)$ $\textstyle =$ $\displaystyle C(t) + n(t) = \alpha T(t-t_0) + n(t)$  
  $\textstyle \iff$    
$\displaystyle \tilde h(f)$ $\textstyle =$ $\displaystyle \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$:

\begin{displaymath}
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}.
\end{displaymath} (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$:
\begin{displaymath}
4 \int_0^{f_{\rm Ny}} df \; { \vert \tilde T(f)\vert^2 \over
S_h(\vert f\vert)} = 1.
\end{displaymath} (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:
$\displaystyle \Delta f_1$ $\textstyle =$ $\displaystyle \{ f \; \vert \; 0< f < f_1 \}$  
$\displaystyle \Delta f_2$ $\textstyle =$ $\displaystyle \{ f \; \vert \; f_1 < f < f_2 \}$  
$\displaystyle \cdots$      
$\displaystyle \Delta f_p$ $\textstyle =$ $\displaystyle \{ 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:
\begin{displaymath}
4 \int_{\Delta f_i} df \; { \vert \tilde T(f)\vert^2 \over S_h(\vert f\vert)} = {1 \over p}
\end{displaymath} (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$.
\begin{figure}\begin{center}\begin{picture}(300,20)(0,0)
\put(0,10){\line(1,0){3...
...(-10,0){$f=0$}
\put(290,0){$f=f_{\rm Ny}$}
\end{picture}\end{center}\end{figure}

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:

\begin{displaymath}
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.
\end{displaymath} (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:
\begin{displaymath}
\sigma = \langle S_i^2 \rangle - \langle S_i \rangle^2 = {1 \over p}.
\end{displaymath} (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

\begin{displaymath}
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}}
\end{displaymath} (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
\begin{displaymath}
\Delta S_i \equiv S_i-S/p.
\end{displaymath} (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
$\displaystyle \int dy_1 \cdots dy_p$ $\textstyle \bar P (y_1,\cdots,y_p)$ $\displaystyle F(y_1,\cdots,y_p)$  
  $\textstyle =$    
$\displaystyle \int dx_1 \cdots dx_p$ $\textstyle P (x_1,\cdots,x_p)$ $\displaystyle 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
\begin{displaymath}
\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 ).
\end{displaymath} (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
$\displaystyle W$ $\textstyle =$ $\displaystyle x_1 + \cdots + x_p$  
$\displaystyle z_1$ $\textstyle =$ $\displaystyle x_1 - W/p$  
  $\textstyle \cdots$    
$\displaystyle z_{p-1}$ $\textstyle =$ $\displaystyle x_{p-1} - W/p$ (6.24.116)

which can be inverted to yield
$\displaystyle x_1$ $\textstyle =$ $\displaystyle z_1 + W/p$  
  $\textstyle \cdots$    
$\displaystyle x_{p-1}$ $\textstyle =$ $\displaystyle z_{p-1} + W/p$ (6.24.117)
$\displaystyle x_p$ $\textstyle =$ $\displaystyle W/p - z_1 - \cdots - z_{p-1}$  

The Jacobian of this coordinate transformation is:
\begin{displaymath}
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]
\end{displaymath} (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

$\displaystyle \bar P (\Delta S_1,\cdots,\Delta S_p)$ $\textstyle =$ $\displaystyle \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}$  
$\displaystyle \times \delta(z_1-\Delta S_1)$ $\textstyle \cdots$ $\displaystyle \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
\begin{displaymath}
(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
\end{displaymath} (6.24.120)

and thus the integral yields
$\displaystyle \bar P (\Delta S_1,\cdots,\Delta S_p)$ $\textstyle =$ $\displaystyle \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)$  
  $\textstyle =$ $\displaystyle (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

\begin{displaymath}
\langle \Delta S_i^2 \rangle = {1 \over p}\left( 1 - {1 \over p} \right)
\end{displaymath} (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

\begin{displaymath}
r^2 = \sum_{i=1}^p (\Delta S_i)^2.
\end{displaymath} (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

\begin{displaymath}
\langle r^2 \rangle = p \langle \Delta S_i^2 \rangle = 1 - {1 \over p}
\end{displaymath} (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
$\displaystyle P(r>R)$ $\textstyle =$ $\displaystyle \int_{r^2 > R^2} \bar P$  
  $\textstyle =$ $\displaystyle (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}$:
$\displaystyle P(r>R)$ $\textstyle =$ $\displaystyle (2 \pi \sigma)^{(1-p)/2} \int_{r^2 > R^2}
{\rm e}^{-r^2/2\sigma} du_1 \cdots du_{p-1}$  
  $\textstyle =$ $\displaystyle \Omega_{p-2} \int_R^\infty r^{p-2} {\rm e}^{-r^2/2\sigma} dr$  
  $\textstyle =$ $\displaystyle {1 \over \Gamma({p-1 \over 2}) } \int_{R^2/2 \sigma}^\infty x^{(p-3)/2} {\rm e}^{-x} dx$  
  $\textstyle =$ $\displaystyle 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:

$\displaystyle h(t)$ $\textstyle =$ $\displaystyle C(t) + n(t) = \alpha T_0(t-t_0) + \beta T_{90}(t-t_0) +
n(t)$  
  $\textstyle \iff$    
$\displaystyle \tilde h(f)$ $\textstyle =$ $\displaystyle \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
\begin{displaymath}
T = \left( {\alpha \over \sqrt{\alpha^2 + \beta^2}} \right) ...
... \left( {\beta \over
\sqrt{\alpha^2 + \beta^2}} \right) T_{90}
\end{displaymath} (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

\begin{displaymath}
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}.
\end{displaymath} (6.24.129)

The templates for the individual phases are normalized as before:
\begin{displaymath}
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.
\end{displaymath} (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:
\begin{displaymath}
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.
\end{displaymath} (6.24.131)

The mean value of the signal in each frequency band is
\begin{displaymath}
\langle S_i \rangle = \langle S \rangle /p = (\alpha + i \beta)/p,
\end{displaymath} (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:
\begin{displaymath}
\langle \vert S_i \vert ^2 \rangle - \vert \langle S_i \rangle \vert^2 = {2 \over p}.
\end{displaymath} (6.24.133)

The signal values are now characterized by the probability distribution
\begin{displaymath}
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}}.
\end{displaymath} (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:
\begin{displaymath}
\Delta S_i \equiv S_i-S/p.
\end{displaymath} (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
\begin{displaymath}
\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)
\end{displaymath} (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:
\begin{displaymath}
\langle \vert\Delta S_i\vert^2 \rangle = {2 \over p}\left( 1 - {1 \over p} \right).
\end{displaymath} (6.24.137)

The $r^2$ statistic is now defined by
\begin{displaymath}
r^2 = \sum_{i=1}^p \left\vert \Delta S_i \right\vert ^2.
\end{displaymath} (6.24.138)

and has an expectation value which as twice as large as in the single-phase case:
\begin{displaymath}
\langle r^2 \rangle = p \langle \left\vert \Delta S_i \right\vert^2 \rangle = 2 - {2 \over p}.
\end{displaymath} (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
$\displaystyle P(r>R)$ $\textstyle =$ $\displaystyle (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$  
  $\textstyle =$ $\displaystyle (2 \pi \sigma)^{(1-p)} \int_{r^2 > R^2}
{\rm e}^{-r^2/2\sigma} du_1 \cdots du_{p-1}$  
  $\textstyle =$ $\displaystyle 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:
\begin{displaymath}
\langle \chi^2 \rangle = 2 p - 2
\end{displaymath} (6.24.141)

the relationship between the $r^2$ and $\chi^2$ statistic may be obtained by comparing equations ([*]) and ([*]), giving
\begin{displaymath}
\chi^2 = p r^2.
\end{displaymath} (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}$.


next up previous contents
Next: How does the test Up: GRASP Routines: Gravitational Radiation Previous: Vetoing techniques (time domain   Contents
Bruce Allen 2000-11-19