FINDCHIRP: An algorithm for detection of gravitational waves from inspiraling compact binaries

Bruce Allen, Warren G. Anderson, Patrick R. Brady, Duncan A. Brown, and Jolien D. E. Creighton Department of Physics, University of Wisconsin–Milwaukee, P.O. Box 413, Milwaukee, Wisconsin 53201, USA Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), Callinstraße 38, D-30167 Hannover, Germany Theoretical Astrophysics 130-33, California Institute of Technology, Pasadena, California 91125, USA LIGO Laboratory, California Institute of Technology, Pasadena, California 91125, USA Department of Physics, Syracuse University, Syracuse, New York 13244, USA (Received 4 August 2011; published 19 June 2012)


I. INTRODUCTION
For the detection of a known signal it is well-known that, in the presence of stationary and Gaussian noise, the use of a matched filter is the optimal detection strategy [1]. Such signals are the anticipated gravitational waveforms from the inspiral of binary systems containing compact objects, i.e., binaries whose components are neutron stars and/or black holes [2,3]. These gravitational-wave sources are among the most promising targets for ground-based gravitational-wave detectors, such as the Laser Interferometer Gravitationalwave Observatory (LIGO) [4] and Virgo [5].
Some practical complications arise for the gravitationalwave detection problem, however, because: (i) the signal is not precisely known-it is parametrized by the binary companions' masses, an initial phase, the time of arrival, and various parameters describing the distance and orientation of the system relative to the detector that can be combined into a single parameter we call the ''effective distance,'' and (ii) the detector noise is not perfectly described as a stationary Gaussian process. Standard techniques for extending the simple matched filter to search over the unknown parameters involve using a quadrature sum of matched-filter outputs for orthogonal-phase waveforms (thereby eliminating the unknown phase), use of Fourier transform to efficiently apply the matched filters for different times of arrival, and use of a bank of templates to cover the parameter space of binary companion masses [6][7][8][9][10][11][12]. Methods for making the matched filter more robust against non-Gaussian noise artifacts, e.g., by examining the relative contributions of frequency-bandlimited matched-filter outputs (vetoing those transients that produce large matched-filter outputs but have a timefrequency decomposition that is inconsistent with the expected waveform), have also been explored [13] and found to be effective.
The FINDCHIRP algorithm is the implementation of matched filtering used by the LIGO Scientific Collaboration (LSC) and Virgo's offline searches for gravitational waves from low-mass (2M < M ¼ m 1 þ m 2 < 35M ), high-mass (25M < M < 100M ), and primordial black hole (0:2M < M < 2M ) coalescing compact binaries [14][15][16][17][18][19][20][21][22][23][24][25]. FINDCHIRP has also been used to search for supermassive black holes in data from the LISA Mock Data Challenge [26][27][28] and for comparisons with numerical relativity waveforms for both high-mass and low-mass binaries [29]. Several aspects of the algorithm have been described in passing in the above references, but here we provide a detailed and comprehensive description of our algorithm as used in the LSC and Virgo's search for coalescing compact binaries.
The FINDCHIRP algorithm is the part of the search that (i) computes the matched-filter response to the interferometer data for each template in a bank of templates, (ii) computes a chi-squared discriminant [13] (if needed) to reject instrumental artifacts that produce large spurious responses of the matched filter but otherwise do not resemble an expected signal, and (iii) selects candidate events or triggers based on the matched-filter and chisquared outputs. This is a fundamental part of the search for coalescing compact binaries, but the search also consists of several other important steps such as data selection and conditioning, template bank generation, rejection of candidate events by vetoes based on auxiliary instrumental channels and multidetector coincidence, and computation of the false alarm rate of candidate triggers.
The entire search pipeline, which is a transformation of raw interferometer data into candidate events, contains all these aspects. The specific details of the pipeline used to search for gravitational waves depends on the target population (e.g., a triggered search for gamma-ray bursts or an all-sky search for low-mass compact binaries) and the particular ''science run'' of the LIGO and/or Virgo detectors. We refer to Refs. [14][15][16]18,19,[21][22][23][24][25][30][31][32] for descriptions of the pipelines that have used the algorithms described here. This paper is not intended to provide documentation for our implementation of the FINDCHIRP algorithm. (This can be found in Refs. [33,34].) Rather, this paper is intended to describe the algorithm itself.

II. NOTATION
Our conventions for the Fourier transform are as follows. For continuous quantities, the forward and inverse Fourier transforms are given bỹ xðtÞe À2ift dt; (2.1a) and xðtÞ ¼ respectively, soxðfÞ is the Fourier transform of xðtÞ.
If these continuous quantities are discretized so that x½j ¼ xðjÁtÞ where 1=Át is the sampling rate and j ¼ 0; . . . ; N À 1 are N sample points, then the discrete approximation to the forward and inverse Fourier transforms arex where Áf ¼ 1=ðNÁtÞ andx½k is an approximation to the value of the continuous Fourier transform at frequency kÁf:x½k %xðkÁfÞ for 0 k bN=2c andx½k % xððk À NÞÁfÞ for bN=2c < k < N (negative frequencies).
Here bac means the greatest integer less than or equal to a. The DC component is k ¼ 0 and, when N is even, k ¼ N=2 corresponds to the Nyquist-frequency. Notice that our convention is to have the Fourier com-ponentsx½k normalized so as to have the same units as the continuous Fourier transformxðfÞ, i.e., the discretized versions of the continuous forward and inverse Fourier transforms carry the normalization constants Át and Áf respectively. Numerical packages instead compute the discrete Fourier transform (DFT) where the minus sign in the exponential refers to the forward DFT and the positive sign refers to the reverse 1 DFT. The DFT is efficiently implemented via the fast Fourier transform (FFT) algorithm. Thus, it is important to write the most computationally demanding terms in the form of Eq. (2.3) so that this computation can be done most efficiently. Throughout this paper we will reserve the indices j to be a time index (which labels a particular time sample), k to be a frequency index (which labels a particular frequency bin), m to be an index over a bank of templates, and n to be an index over analysis segments. Thus, for example, the quantity z m;n ½j will be the jth sample of analysis segment n of the matched-filter output for the mth template, and z m;n ½k will be the kth frequency bin of the Fourier transform of the matched-filter output for the same template and analysis segment.

III. WAVEFORM
For first-generation gravitational-wave detectors (such as Initial LIGO), the gravitational-wave signal from a compact binary inspiral waveform can be accurately modeled by the restricted post-Newtonian waveform [35] below approximately M $ 12M [36] for nonspinning bodies. At higher masses, resummation techniques such as the effective one body (EOB) waveforms better reproduce the waveforms computed by numerically solving the full nonlinear Einstein equations [37,38]. In either case, the two polarizations of the gravitational-wave produced by such a system exhibit a monotonically-increasing frequency and amplitude as the orbital motion radiates away energy and decays. The waveform, often called a chirp waveform, is given for t < t c by where D is the distance from the source, is the angle between the direction to the observer and the orbital angular momentum axis of the binary system, is the total mass of the two companions, the reduced mass is ¼ m 1 m 2 =M, and ¼ =M) is the chirp mass, and ðt À t c ; M; Þ is the orbital phase of the binary (whose evolution also depends on the individual masses of the binary companions) [39,40]. Here, t c and c are the time and phase of the binary coalescence when the waveform is terminated, known as the coalescence time and coalescence phase ð0; M; Þ ¼ 0. Note that in Ref. [39] the integration constant t c and c are implicitly contained in the orbital phase ðtÞ. We have chosen to express them explicitly.
The coalescence phase and time may lack a formal definition (as in the case of numerical waveforms), however for a detection algorithm, the formal definition of the coalescence time is not critical, as long as any offset caused by the choice of coalescence time is constant between the detectors in the network. For restricted post-Newtonian waveforms, we define t c to be the time at which the gravitational-wave frequency becomes infinite within the post-Newtonian formalism. The Initial LIGO science runs (S1-S5) used the second-post-Newtonian waveform, whereas the Enhanced LIGO (S6) analysis uses post-Newtonian waveforms computed to 3.5th order. In this paper we illustrate the FINDCHIRP algorithm using second-order post-Newtonian waveforms. The extension of the algorithm to EOB waveforms (and other waveforms computed in the time-domain) is described in Appendix A and the extension to higher post-Newtonian orders is described in Appendix E. We assume that the binary system's center-of-mass is at rest with respect to the detector frame (otherwise this motion ''redshifts'' the binary's masses). We also assume that the time that the signal is in the detector's sensitive bandwidth is short compared to 24h, so that the detector does not change orientation significantly as the earth rotates.
The gravitational-wave strain induced in a particular detector depends on the detector's antenna response to the two polarizations of the gravitational waveform. The induced strain on the detector is given by where t 0 is the termination time (the time at the detector at which the coalescence occurs, i.e., the detector time when the gravitational-wave frequency, according to restricted post-Newtonian approximation, becomes infinite), t 0 À t c is the propogation time from the source to the detector, and F þ and F Â are the antenna response functions for the incident signal; these functions depend on the location of the source with respect to the reference frame of the detector, which is described by the right ascension and declination of the source, (, ), the arrival time at the detector, t 0 , and on the polarization angle c [41]. The antenna response functions are very nearly constant in time over the duration of the short inspiral signal. Thus the strain on a particular detector can be written as where 0 is the termination phase which is related to the coalescence phase by and is the effective distance of the source. The effective distance of the source is related to the true distance of the source by several geometrical factors that relate the source orientation to the detector orientation. Because the location and orientation of the source are not likely to be known when filtering data from a single detector, it is convenient to combine the geometric factors with the true distance to give a single observable, the effective distance. For an optimally oriented source (one that is directly overhead and is orbiting in the plane of the sky) the effective distance is equal to the true distance; for sub-optimally-oriented sources the effective distance is greater than the true distance. (The location and distance can be estimated using three or more detectors, but we do not consider this here.) Equation (3.3a) gives a waveform that is used as a template for a matched filter. Since FINDCHIRP implements the matched filter via a FFT correlation, it is beneficial to write the Fourier transform of the template and implement it directly rather than taking the FFT of the time-domain waveform of Eq. (3.3a). A frequency-domain version of the waveform can be obtained via the stationary phase approximation [6,42,66]. For f > 0 one has and É has been written to second post-Newtonian order. For f < 0 one hashðfÞ ¼h Ã ðÀfÞ. This template waveform has been expressed in terms of several factors: (1) An overall distance factor involving the effective distance, D eff . For a template waveform, we are free to choose this effective distance as convenient, and in the FINDCHIRP code it is chosen to be 1 Mpc. (2) A constant (in frequency) factor A 1 Mpc ðM; Þ, which has dimensions of ðtimeÞ À1=6 , that depends only on the total and reduced masses, M and , of the particular system. (3) The factor f À7=6 which does not depend on the system parameters. And (4) a phasing factor involving the phase Éðf; M; Þ which is both frequency-dependent and dependent on the system's total and reduced masses. We will see below that an efficient implementation of the matched filter will make use of this factorization of the stationary phase template.
In order to construct a waveform template we need to know how long the binary system will radiate gravitation waves in the sensitivity band of LIGO. A true inspiral chirp waveform would last tens of Myr, but the amount of time that the binary system spends radiating gravitational waves with a frequency above some low-frequency cutoff f low is short: the duration of the chirp or chirp time from a given frequency f low is given to second post-Newtonian order by Eq. (3.3) of Ref. [43] as For the Initial and Enhanced LIGO detectors, FINDCHIRP uses f low ¼ 40 Hz. Higher-mass systems coalesce much more quickly (from a given f low ) than lower mass systems. A search for low-mass systems, such as primordial black holes of mass 0:1M , can require very long waveform templates (of the order of tens of minutes) which can result in a significant computational burden. There is also a high-frequency cutoff for the inspiral waveform. Physically, at some high frequency a binary system will terminate its secular inspiral and the orbit will decay on a dynamical time scale, though identifying the precise frequency is difficult except in the extreme mass ratio limit ! 0. In this limit, that of a test mass orbiting a Schwarzschild black hole, the frequency is known as the innermost stable circular orbit or ISCO. The ISCO gravitational-wave frequency is However, before reaching this frequency, the binary components will be orbiting with sufficiently high orbital velocities that the higher-order corrections to the post-Newtonian waveform will become significant. We regard Eq. (3.6) as an upper limit on the frequency that can be regarded as representing an ''inspiral'' waveform-not as the frequency to which we can trust our inspiral waveform templates. With this caveat, we nevertheless use this as a high-frequency cutoff for the inspiral template waveforms (should this frequency be less than the Nyquist frequency f Nyquist ¼ 1=2Át of the data, where Át is the sample rate). For the lowest mass binary systems (binary neutron stars or subsolar mass black holes) the post-Newtonian template waveforms are reliable within the sensitive band of Initial LIGO so the precise choice of the high-frequency cutoff is not important [36]. For higher-mass systems (containing black holes) the effects of higher-order corrections to the post-Newtonian waveform can be significant. For frequency-domain templates, this can be addressed by the use of ''pseudo post-Newtonian'' terms in the waveform phasing and/or different choices of the high-frequency cutoff [44,45]. When using time-domain EOB waveforms tuned to numerical relativity simulations, the upperfrequency cutoff is set by the frequency of the fundamental l ¼ 2, m ¼ 2 quasinormal ringdown mode [46]. In this case the FINDCHIRP algorithm uses the minimum of this or the Nyquist frequency.

IV. MATCHED FILTER
The matched filter is the optimal filter for detecting a known waveform in stationary Gaussian noise. Suppose that nðtÞ is a stationary Gaussian noise process with one-sided power spectral density S n ðfÞ defined by hñðfÞñ Ã ðf 0 Þi ¼ 1 2 S n ðjfjÞðf À f 0 Þ. Then the matched-filter output of a data stream sðtÞ-which now may contain detector noise alone, sðtÞ ¼ nðtÞ, or a signal in addition to the noise, sðtÞ ¼ nðtÞ þ hðtÞ-with a filter template h template ðtÞ is where the signal h template ðtÞ is implicitly taken to depend on a termination time t 0 as in (3.3a). Notice that the use of a FFT will allow one to search for all possible termination times t 0 efficiently. However, the waveforms described above have additional unknown parameters. These are (i) the amplitude (or effective distance to the source), (ii) the coalescence phase, and (iii) the binary companion masses. The amplitude simply sets a scale for the matchedfilter output, and is unimportant for matched-filter templates (these can be normalized). The ''best match'' unknown phase 0 can be found by maximizing xðt 0 Þ over 0 . xðt 0 Þ can be written as xðt 0 Þ ¼ x re ðt 0 Þ cos2 0 þ x im ðt 0 Þ sin2 0 where x re;im ðt 0 Þ are (respectively) the values of (4.1) with 0 ¼ 0, and with 0 ¼ 0 and Re ! Im. Setting the derivative with regard to 0 to vanish, one finds that at the maximum Hence this maximum value is given by the (more computationally efficient) modulus of the complex filter output where ½h Ã template ðfÞ 0 ¼ ½h Ã template ðfÞ t 0 ¼0; 0 ¼0 . This template is obtained from Eq. (3.4) by deleting the first two terms on the right-hand side of (3.4c).
To search over all the possible binary companion masses it is necessary to construct a bank of matched-filter templates laid out on the ðm 1 ; m 2 Þ plane sufficiently finely that any true system masses will produce a waveform that is close enough to the nearest template. There are wellknown strategies for constructing such a bank [8][9][10][11]. For our purposes, we shall simply introduce an index m ¼ 0; . . . ; N T À 1 labeling the particular waveform template h m ðtÞ in the bank of N T waveform templates.
By convention, the waveform templates are constructed for systems with an effective distance of D eff ¼ 1 Mpc. To construct a signal-to-noise ratio, a normalization constant for each template is computed The quantity m is a measure of the sensitivity of the instrument. For sðtÞ that is purely stationary and Gaussian noise, h<z m ðtÞi ¼ h=z m ðtÞi ¼ 0 and one finds that hz 2 m ðtÞi¼hz 2 m ðtÞi¼ 2 m , while for a detector output that includes a signal at distance D eff ; sðtÞ ¼ nðtÞ þ ðD eff =1 MpcÞ À1 h 1 Mpc;m ðtÞ, hz m ðt 0 Þi¼1Mpc 2 m =D eff . Thus the quantity is the amplitude signal-to-noise ratio of the (quadrature) matched filter. Note that while this is called a signal-to-noise ratio, in the absence of a signal the expectation value is not unity: h 2 m i ¼ 2. It is highly unlikely to obtain m ) 1 for purely stationary and Gaussian noise so a detection strategy usually involves setting a lower threshold on m to identify event candidates. For such candidates, a biased estimate of the effective distance to the candidate system isD eff ¼ ð m = m ÞMpc.
The goal of the FINDCHIRP algorithm is largely to construct the quantity m ðtÞ and to identify the values of the parameters t 0 , 0 , and m that maximize it.

V. DETECTOR OUTPUT AND CALIBRATION
LIGO records several interferometer channels. The gravitational-wave channel (the primary channel for searching for gravitational waves) is formed from the output of a photo-diode at the antisymmetric (or ''dark'') port of the interferometer [47]. This output is used as an error signal for a feedback loop that is needed to keep various optical cavities in the interferometer in resonance or ''in-lock.'' Hence it is often called the error signal eðtÞ. The error signal is not an exact measure of the differential arm displacements of the interferometer so it does not correspond to the gravitational-wave strain. Rather it is part of a linear feedback loop that controls the position of the interferometer mirrors. A gravitational-wave strain-equivalent output, called sðtÞ above, can be obtained from the error signal eðtÞ via a linear filter. This is called calibration. Details on the calibration of the LIGO interferometers can be found in [48,49]. In the frequencydomain, the process of calibration can be thought of as multiplying the error signal by a complex response function, RðfÞs ðfÞ ¼ RðfÞẽðfÞ: The FINDCHIRP algorithm can analyze either calibrated data sðtÞ or it can calibrate the error signal eðtÞ in the frequencydomain. With time-domain calibrated data now available, it is more convenient to use sðtÞ so in this paper we will present the algorithm in terms of strain data rather than the error signal. The detector output is not a continuous signal but rather a time series of samples of sðtÞ taken with a sample rate of 1=Át ¼ 16384 Hz where Át is the sampling interval. Thus, rather than sðtÞ, the input to FINDCHIRP is a discretely sampled set of values s½j ¼ sðt start þ jÁtÞ for some large number of points. The start of the data sample is at time t start . Data from the detector is divided into science segments which are time epochs when the instrument was in-lock and exhibiting normal behavior. However, for practical computational reasons, these science segments are not normally processed as a whole, but are divided into smaller amounts. In this paper we shall call the amount of data processed a data block of duration T block . The data block must be long enough to form a reliable noise power spectral estimate (see below), but not so long as to exhaust a computer's memory or to encompass significant nonstationary changes in the detector noise.
The number of points in a data block is further subdivided into N S overlapping data segments or just segments (not to be confused with the science segments described above) of duration T. The duration of the segment is always an integer multiple of the sample rate Át, so the number of points in a segment N ¼ T=Át is an integer. These segments are used to construct an average noise power spectrum and to perform the matched filtering. The segments are overlapped so that the first segment consists of the points s½j for j ¼ 0; . . . ; N À 1, the second consists of the points j ¼ Á; . . . ; Á þ N À 1 where Á is known as the stride, and so on until the last segment which consists of the points j ¼ ðN S À 1ÞÁ; . . . ; ðN S À 1ÞÁ þ N À 1. Note that All of the LSC's searches choose to overlap the segments by 50% so that the stride is Á ¼ N=2 (and N is always even) and hence there are N S ¼ 2ðT block =TÞ À 1 segments. The values of T block , T, Á, and N S must be commensurate so that these relations hold. Values typically used in the LIGO and/or Virgo analyses for low-mass binaries The FINDCHIRP algorithm (4.2) implements the matched filter by an FFT correlation. Thus the discrete Fourier transforms of the individual data segments, n, is purely real, the components 0 < k bðN À 1Þ=2c are all positive-frequency components corresponding to frequencies kÁf where Áf ¼ 1=ðNÁtÞ, and the components bN=2c < k < N are all negative-frequency components corresponding to frequencies ðk À NÞÁf. If N is even (as it always is for the FINDCHIRP algorithm) then there is also a purely real Nyquist-frequency component in bin k ¼ N=2 corresponding to the Nyquist frequency AENÁf=2 ¼ AE1=ð2ÁtÞ. Recall bac is the greatest integer less than or equal to a. Note that because the data is real, the discrete Fourier transform of it satisfiess Ã n ½k 1 s n ½N À k. Thus, the FINDCHIRP algorithm only stores the frequency components k ¼ 0; . . . ; bN=2c. These can be efficiently computed using a real-to-half-complex forward FFT [50].
If the error signal rather than detector strain is analyzed, we simply replace s with e in Eq. (5.3). Then the detector strain for segment n can be computed by calibrating the error signal where R½k ¼ RðkÁfÞ is the complex response function. As before, sinces n ½k must be the Fourier transform of some real time series, only the frequency components k ¼ 0; . . . ; bN=2c need to be computed. LIGO is sensitive to strains that are smaller than $10 À20 , while the error signal is designed to have typical values much closer to unity. Often the FINDCHIRP algorithm will require quantities that are essentially squares of the measured strain (e.g., the power spectrum described in the next section). To avoid floating-point over-or under-flow problems, the strain can simply be rescaled by a dynamic-range factor . If strain data is input, it is immediately scaled by the factor , so s½j is used rather than s½j. If the error signal is input, the dynamic-range factor is applied to the response function instead, so that R is used rather than R. Choosing a value of $ 10 20 will keep all quantities within representable single-precision IEEE 781 floating-point numbers. It is important to keep track of the factor to make sure it cancels out in all of the results. Essentially this is achieved by multiplying all quantities with ''units'' of strain by the factor within the implementation of the FINDCHIRP algorithm. Thus, in addition to the strain data or response function, the signal template must also be scaled by . On conventional CPUs, storing quantities in single-precision reduces the performance cost of moving quantities to and from memory (which can be dominant). In addition, some CPUs and most graphics processing units (GPUs) perform single-precision arithmetic at least twice as fast as double-precision arithmetic. Therefore it is advantageous to use floating-point (singleprecision) operations for the FINDCHIRP algorithm.
If frequency-domain calibration of the error signal is desired (e.g., if the effects of calibration error are being investigated), the following replacements need to be made in the formulas in this paper:s½k ! R½kẽ½k, S½k ! jR½kj 2 S e ½k, and Q ! jR½kj À2 Q e ½k where S e ½k and Q e ½k are the power spectral density and inverse truncated power spectral density computed from the error signal rather than the calibrated strain data (see Secs. VI and VII).

VI. AVERAGE POWER SPECTRUM
Part of the matched filter involves weighting the data by the inverse of the detector's power spectral density. The detector's power spectrum must be obtained from the detector output. The most common method of power spectral estimation is Welch's method. Welch's method [51] for obtaining the average power spectrum S of the data is is a normalized periodogram for a single segment n which is the modulus-squared of the discrete Fourier transform of windowed data. The data window is given by w½j and W is a normalization constant j¼0 w 2 ½j: (6.2b) FINDCHIRP allows a variety of possible windows, but a Hann window (see, e.g., [52]) is the default choice used by FINDCHIRP. We call this average power spectrum the mean average power spectrum. If uncalibrated data e n ½j is being used rather than s n ½j, then the power spectrum of the detector strain-equivalent noise is related to the power spectrum of the uncalibrated data by S½k ¼ jR½kj 2 S e ½k. The problem with using Welch's method for power spectral estimation is that for detector noise containing significant excursions from ''normal'' behavior (due to instrumental glitches or unexpectedly strong gravitationalwave signals), the mean used in Eq. (6.1) can be significantly biased by the excursion. An alternative that is pursued in the FINDCHIRP algorithm is to replace the mean in Eq. (6.1) by a median, which is a more robust estimator of the average power spectrum where is a required correction factor. When ¼ 1, the expectation value of the median is not equal to the expectation value of the mean in the case of Gaussian noise; hence the factor is introduced to ensure that the same power spectrum results for Gaussian noise. In Ref. [53] and in Appendix B it is shown that if the set fP 0 ½k; P 1 ½k; . . . ; P N S À1 ½kg are independent exponentiallydistributed random variables (as expected for Gaussian noise) then is the correction factor. We call this median estimate of the average power spectrum, corrected by the factor , the median average spectrum. Unfortunately this result is not exactly correct either. Because the segments used to form the individual sample values P n ½k of the power at a given frequency are somewhat overlapping (unless Á ! N), they are not independent random variables, as was assumed in Appendix B. (This is somewhat mitigated by the windowing of the segments of data.) Although the effect is not large, and simply amounts to a slight scaling of what is meant by signal-to-noise ratio, we are led to propose a variant of the median method in which the n ¼ 0; . . . ; N S À 1 overlapping segments are divided into even segments (for which n is even) and the odd segments (for which n is odd). If the stride is Á ! N=2 then no two even segments will depend on the same data so the even segments will be independent; similarly the odd segments will be independent. The average power spectrum can be estimated by taking the weighted mean of the median power spectrum of the ðN S þ 1Þ=2 even segments and the median power spectrum of the ðN S À 1Þ=2 odd segments, each of which are corrected by a factor appropriate for the sample median with, respectively, ðN S AE 1Þ=2 samples. We call this the median-mean average spectrum. Like the median spectrum it is not overly sensitive to a single glitch (or strong gravitationalwave signal).
The FINDCHIRP algorithm can compute the mean average spectrum, the median average spectrum, or the medianmean average spectrum. Traditionally the median average spectrum has been used though we expect that the medianmean average spectrum will be adopted in the future.

VII. DISCRETE MATCHED FILTER
The discretized version of Eq. (4.2) is simply where inh both t 0 and 0 are set to zero. Element j of z n;m ½j corresponds to the matched-filter output (4.2) for s n ½kh Ã 1 Mpc;m ½k 2 S½k k low k bðN À 1Þ=2c 0 bðN À 1Þ=2c < k < N: That is to say, the DC, Nyquist, and negative frequency components are all set to zero, as are all frequencies below some low-frequency cutoff f low ¼ k low Áf (set to some frequency lower than the detector's sensitive band). The low-frequency cutoff limits the time duration of the inspiral template as described below. Our task is to obtain an efficient computation of the factors making upz n;m ½k. Note that there needs to be one reverse FFT performed per segment per template. It is desirable that this (unavoidable) computational cost dominates the evaluation of the matched filter, so we wish to make the computation cost of the calculation ofz n;m ½k for all k and fixed n; m to be less than the computation cost of a FFT. We will consider this in the next section.
One subtlety in the construction of the matched filter is the issue of filter wraparound. The matched filter of Eq. (7.1) can be thought of as digital correlation of a filter h 1 Mpc;m ½j with some suitably over-whitened data stream (the data divided by the noise power spectrum). To simplify the discussion, first assume S½k ¼ 1. Although h 1 Mpc;m ½k is generated in the frequency domain via the stationary phase approximation, we can imagine that it came from a time-domain signal h 1 Mpc;m ½j of duration T chirp;m that is given by Eq. (3.5a) for the low-frequency cutoff f low . By convention of template generation, the time of coalescence corresponds to j ¼ 0. Thus, the entire chirp waveform is nonzero only from t ¼ t 0 À T chirp;m to t ¼ t 0 . Because the discrete Fourier transform presumes that the data is periodic, this is represented by having the chirp begin at point j ¼ N À T chirp;m =Át and end at point j ¼ N À 1. Thus h 1 Mpc;m ½j ¼ 0 for j ¼ 0; . . . ; N À 1 À T chirp;m =Át. The correlation of h 1 Mpc;m ½j with the interferometer data will involve multiplying the T chirp;m =Át points of data before a given time with the T chirp;m =Át points of the chirp. When this is performed by the FFT correlation, this means that the first T chirp;m =Át points of the matched-filter output involve data at times before the start of the segment, which are interpreted as the data values at the end of the segment (since the FFT assumes that the data is periodic). Hence the first T chirp;m =Át points are of the correlation are invalid and must be discarded. That is, of the N points of z n;m ½j in Equation (7.1), only the points j ¼ T chirp;m =Át; . . . ; N À 1 are valid. Recall that the analysis segments of data are overlapped by an amount N À Á: this is to ensure that the matched-filter output is continuous (except at the very beginning of a data block). That is, only points j ¼ T chirp;m =Át; . . . ; N À 1 of z 0;m are valid and only points j ¼ T chirp;m =Át; . . . ; N À 1 of z 1;m are valid, but points j ¼ Á; . . . ; N À 1 of z 0;m ½j correspond to points j ¼ 0; . . . ; N À Á À 1 of z 1;m ½j, and these can be used instead. Therefore FINDCHIRP must ensure that the amount that the data segments overlap, N À Á points, is greater than or equal to the duration, T chirp;m =Át points, of the filter: T chirp;m =Át N À Á.
The quantity that needs to be computed in Eq. (7.1) is more than just a correlation of the data s n ½j with the filter h 1 Mpc;m ½j: it also involves a convolution of the data with the response function (if we are performing frequencydomain calibration) and the inverse of the power spectrum. The interferometer has a relatively short impulse response, so this convolution will only corrupt a short amount of data (though now at the end as well as at the beginning of a analysis segment). However, the inverse of the power spectrum has many very narrow line features that act as sharp notch filters when applied to the data. These filters have an impulse response that is as long as the reciprocal of the resolution of the frequency series, which is set by the amount of data used to compute the periodograms that are used to obtain the average spectrum. Since this is the same duration as the analysis segment duration, the convolution of the data with the inverse power spectrum corrupts the entire matched-filter output.
To resolve this problem we apply a procedure to coarsegrain the inverse power spectrum called inverse spectrum truncation. Our goal is to limit the amount of the matched filter that is corrupted due to the convolution of the data with the inverse power spectrum. To do this we will begin with the time-domain version of the frequency-domain quantity S À1 ½k, truncate it so that it has finite duration, and then find the frequency-domain quantity Q½k corresponding to this truncated time-domain filter. Note that S À1 ½k is real and non-negative, and we want Q½k to share these properties. First, construct the time-domain quantity which can be done via a half-complex-to-real reverse FFT. Since S½k is real and symmetric (S½k ¼ S½N À k), q½j will be real and symmetric (so that q½j ¼ q½N À j). This quantity will be nonzero for all N points, though stronglypeaked near j ¼ 0 and j ¼ N À 1. Now create a truncated quantity q T ½j with a total duration of T spec (T spec =2 at the beginning and T spec =2 at the end) À1 q½j N À T spec =2Át j< N: Since q T is real and symmetric, the discrete Fourier transform of q T will also be real and symmetric, though not necessarily positive. Therefore we construct This quantity is real, positive, and symmetric, as desired.
Multiplying the data by Q½k in the frequency-domain is equivalent to convolving the data with q T ½j in the timedomain twice, which will have the effect of corrupting a duration of T spec of the matched filter z n;m ½j at the beginning and a duration of T spec at the end of the data segment. This is in addition to the duration T chirp;m that is corrupted at the beginning of the data segment due to the correlation with the filter h 1 Mpc;m ½j. Thus the total duration that is corrupted is 2T spec þ T chirp;m , and this must be less than the time that adjacent segments overlap. The net effect of the inverse spectrum truncation is to smear out sharp spectral features and to decrease the resolution of the inverse power spectrum weighting. For simplicity, we normally choose a 50% overlap (so that Á ¼ N=2). Of each data segment the middle half with j ¼ N=4; . . . ; 3N=4 À 1 is assumed to be valid matchedfilter output. Therefore, the inverse truncation duration T spec and the maximum filter duration T chirp;m must satisfy T spec þ T chirp;m T=4 since a time T spec þ T chirp;m is corrupted at the beginning of the data segment.

VIII. WAVEFORM DECOMPOSITION
Our goal is now to construct the quantitỹ z n;m ½k ¼ 4 À2 Q½ks n ½kh Ã 1 Mpc;m ½k (8.1) as efficiently as possible. This quantity must be computed for every segment n, every template m, and every frequency bin in the range k ¼ k low ; . . . ; k high;m À 1 where k low ¼ bf low =Áfc and k high;m is the high-frequency cutoff of the waveform template, which is given by the minimum of the ISCO frequency of Eq. (3.6) and the Nyquistfrequency (recall that the ISCO frequency depends on the binary system's total mass so it is template-dependent). We can factorizez n;m ½k as follows: where A 1 Mpc;m is a template normalization (it needs to be computed once per template but does not depend on k), G m ½k ¼ expðiÉ m ½kÞ with É m ½k is a template phase which must be computed at all values of k for every template (but does not depend on the data segment), and F n ½k is the FINDCHIRP data segment that must be computed for all values of k for each data segment (but does not depend on the template). As before, both t 0 and 0 are set to zero in É. FINDCHIRP first computes and stores the quantities F n ½k for all data segments. Then, for each template m in the bank, the phasing É m ½k is computed once and then applied to all of the data segments (thereby marginalizing the cost of the template generation).
To facilitate the factorization, we rewrite Eq. (3.4a) with The dependence on the data segment is wholly contained in the template-independent quantity F n ½k which is F n ½k ¼ 4 À2 Q½ks n ½kk À7=6 : (8.6) As mentioned earlier, FINDCHIRP computes and stores F n ½k for all segments only once, and then reuses these precomputed spectra in formingz n;m ½k according to Eq. (8.3). The k-dependence on the template is wholly contained in the data-segment-independent quantity G m ½k which is This quantity is known as the FINDCHIRP template. The value of m is also needed in order to normalize z n;m ½j to compute signal-to-noise . It is needs to be computed only once per data block (i.e., only once per power spectrum)-it does not depend on the particular segment within a block or on the particular template in the bank, except in the high-frequency cutoff of the template (if it is less than the Nyquist-frequency). To account for this minimal dependence on the template, the quantity &½k high is precomputed for all values of k high . The division of the matched filter into the data-segmentonly quantity F n ½k and the template-only quantity G m ½k means that FINDCHIRP can efficiently compute the matched filter, or, rather, a quantity that is proportional to it m;n ½j ¼ X Notice that m;n ½j, which is a complex quantity, can be computed using a simple unnormalized reverse FFT of the quantity F n ½kG m ½k. FINDCHIRP computes and stores F n ½k for each of the N S segments in the data block and then, for each template m in the bank, G m ½k is computed and used to filter each of the N S segments. This means that for each data segment and template the computational cost is essentially limited to $N complex multiplications plus one complex reverse FFT. The quantities m;n ½j and z m;n ½j are simply related by a normalization factor z m;n ½j ¼ A 1 Mpc;m m;n ½j: (8.11) Furthermore, the signal-to-noise ratio is related to m;n ½j via 2 m;n ½j ¼ j m;n ½jj 2 =& 2 ½k high;m : (8.12) Rather than applying this normalization to construct the signal-to-noise ratio, FINDCHIRP instead scales the desired signal-to-noise ratio threshold ? to obtain a normalized threshold which can be directly compared to the values j m;n ½jj 2 to determine if there is a candidate event (when j m;n ½jj 2 > 2 ? ). When an event candidate is located, the value of the signal-to-noise ratio can then be recovered for that event along with an estimate of the termination time, where j peak is the point at which j m;n ½jj is peaked; the effective distance of the candidate, D eff ¼ jA 1 Mpc;m j& 2 ½k high;m j m;n ½j peak j Mpc; (8.14b) and the termination phase of the candidate, 2 0 ¼ argðz mn ½j peak Þ ¼ argðÀ m;n ½j peak Þ: (8.14c) The relative sign in the argument arises because z m;n and m;n have opposite sign; this is because A 1Mpc;m , which appears in Eq. (8.11) is negative.

IX. THE CHI-SQUARED VETO
The FINDCHIRP algorithm employs the chi-squared discriminator of Ref. [13] to distinguish between plausible signal candidates and common types of noise artifacts. This method is a type of time-frequency decomposition that ensures that the matched-filter output has the expected accumulation in various frequency bands. (Noise artifacts tend to excite the matched filter at the high-frequency or the low-frequency, but seldom produce the same spectrum as an inspiral.) For data consisting of pure Gaussian noise, the real and imaginary parts of m;n ½j (for a given value of j) are independent Gaussian random variables with zero mean and variance & 2 ½k high;m . If there is a signal present at an effective distance D eff then hRe m;n ½ji¼ ÀA 1Mpc;m & 2 ½k high;m ð1Mpc=D eff Þcos2 0 and hIm m;n ½ji ¼ ÀA 1 Mpc;m & 2 ½k high;m ð1 Mpc=D eff Þ sin2 0 (at the termination time, where 0 is the termination phase). Now consider the contribution to m;n ½j coming from various frequency sub-bands ';m;n ½j ¼ X ðk ' À1Þ k¼k ð'À1Þ F n ½kG m ½ke 2ijk=N (9.1) for ' ¼ 1; . . . ; p. The p sub-bands are defined by the frequency boundaries fk 0 ¼ k low ; k 1 ; . . . ; k p ¼ k high;m g, which are chosen so that a true signal will contribute an equal amount to the total matched filter from each frequency band. This means that the values of k ' must be chosen so that With this choice of bands and for pure Gaussian noise, the real and imaginary parts of ';m;n ½j will be independent Gaussian random variables with zero mean and variance & 2 ½k high;m =p. Furthermore, the real and imaginary parts of ';m;n ½j and ' 0 ;m;n ½j with ' Þ ' 0 will be independent since ';m;n ½j and ' 0 ;m;n ½j are constructed from disjoint bands. Also note that ';m;n ½j: The chi-squared statistic is now constructed from ';m;n ½j as follows: For pure Gaussian noise, 2 is chi-squared distributed with ¼ 2p À 2 degrees of freedom. That ¼ 2p À 2 rather than ¼ 2p results from the fact that the sample mean m;n ½j=p is subtracted from each of values of ';m;n ½j in the sum. However, this subtraction guarantees that, in the presence of a signal that exactly matches the template h 1 Mpc;m (up to an arbitrary amplitude factor and a coalescence phase), the value of 2 is unchanged. Thus, 2 is chisquared distributed with ¼ 2p À 2 degrees of freedom in Gaussian noise with or without the presence of an exactly matched signal.
If there is a small mismatch between a signal present in the data and the template, which would be expected since the templates are spaced on a grid and are expected to provide a close match but not a perfect match to a true signal, then 2 will acquire a small noncentral parameter. This is because the mismatched signal may not shift the mean value of the real parts of f ';m;n g by the same amount (for each '), and similarly the mean values of the imaginary parts of f ';m;n g may not be shifted by the same amounts. The effect on the chi-squared distribution is to introduce a noncentral parameter that is no larger than max ¼ 2 2 m =D 2 eff where D eff is the effective distance of the true signal and is the mismatch between the true signal and the template h 1 Mpc;m , which could be as large as the maximum mismatch of the template bank that is used [8].
Even for small values of (3% is a canonical value), a large value of 2 can be obtained for gravitational waves from nearby binaries. Therefore, one should not adopt a fixed threshold on 2 lest very loud binary inspirals be rejected by the veto. For a noncentral chi-squared distribution with degrees of freedom and a noncentral parameter of , the mean of the distribution is þ while the variance is between 1 and 2 times the mean (the variance equals twice the mean when ¼ 0 and the variance equals the mean for ) ). Thus it is useful to adopt a threshold on the quantity 2 =ð þ Þ, which would be expected to be on the order of unity even for very large signals. The FINDCHIRP algorithm adopts a threshold on the related quantity Ä m;n ½j ¼ 2 m;n ½j p þ 2 m;n ½j : (9.5) Sometimes the quantity r 2 ¼ 2 =p is referred to, rather than Ä, but this quantity does not include the effect of the noncentral parameter.

X. TRIGGER SELECTION
The signal-to-noise ratio threshold is the primary parameter in identifying candidate events or triggers. As we have said, the FINDCHIRP algorithm does not directly compute the signal-to-noise ratio, but rather the quantity m;n ½j given in Eq. (8.10), whose square modulus is then compared to a normalized threshold given by Eq. (8.13). The computational cost of the search is essentially the cost of OðNÞ complex multiplications plus OðN logNÞ operations to perform the reverse FFT of Eq. (8.10), and an additional OðNÞ operations to form the square modulus of m;n ½j for all j. In practice, the computational cost is dominated by the reverse complex FFT.
Triggers that exceed the signal-to-noise ratio threshold are then subjected to a chi-squared test. However, the construction of 2 m;n ½j is much more costly than the construction of m;n ½j simply because p reverse complex FFTs of the form given by Eq. (9.1) must be performed. 2 The cost of performing a chi-squared test is essentially p times as great as the cost of performing the matched filter. FINDCHIRP will only perform the chi-squared test if a threshold-crossing trigger is found. Therefore, if threshold-crossing triggers are rare then the cost of the chi-squared test is small compared to the cost of the filtering. Early LIGO science runs (S1 and S2) used a ''singlestage'' pipeline, in which the chi-squared veto is computed as described here. For reasons of computational cost, later LIGO and/or Virgo analyses used a ''two-stage'' pipeline in which triggers are required to be coincident between at least two different detectors before the chi-squared is computed. This is achieved by disabling the chi-squared test on the first pass of trigger generation on individual detectors and then only applying it on the triggers that survive the coincidence criteria. The chi-squared computation is easily parallelizable and a GPU can be used to perform the calculation. Simple GPU implementations of the chi-squared veto have been shown to reduce the cost of performing the chi-squared by a factor of $20 [54,55]. Given this reduced cost, future LSC analyses are likely to again use the single-stage pipeline for simplicity. Additional speed gains can be realized by moving the entire FINDCHIRP computation onto the GPU.
A true signal in the data is expected to produce a sharp peak in the matched-filter output at almost exactly the correct termination time t 0 (usually within one sample point of the correct time in simulations). For sufficiently loud signals, however, a signal-to-noise threshold may be crossed for several samples even though the correct termination time will have a much greater signal-to-noise ratio than nearby times. Non-Gaussian noise artifacts may produce many threshold-crossing triggers, often for a duration similar to the duration of the inspiral template that is used. In principle, a large impulse in the detector output at sample j 0 can cause triggers for samples j 0 À T spec =Át j j 0 þ ðT spec þ T chirp;m Þ=Át. Rather than record triggers for all samples in which the signal-to-noise threshold is exceeded while the chi-squared test is satisfied, FINDCHIRP has the option of maximizing over a chirp: essentially clustering together triggers that lie within a time T chirp;m . Algorithmically, whenever j m;n ½jj 2 > 2 ? and Ä m;n ½j < Ä ? , a trigger is created with a value of and 2 . If this trigger is within a time T chirp;m after an earlier trigger with a larger value of the signal-to-noise ratio , discard the current trigger (it is clustered with the previous trigger). If this trigger is within a time T chirp;m after an earlier trigger with a smaller signal-to-noise ratio , discard the earlier trigger (the previous trigger is clustered with the current trigger). The result is a set of remaining triggers that are separated by a time of at least T chirp;m . Note that this algorithm depends on the order in which the triggers are selected, i.e., a different set of triggers may arise if the triggers are examined in inverse order of j rather than in order of j. FINDCHIRP applies the conditions as j is advanced from j ¼ N=4 to j ¼ 3N=4 À 1 (i.e., forward 2 If 2 m;n ½j is only required for one particular j then there is a more efficient way to compute it. However, the FINDCHIRP algorithm does not employ the chi-squared test so much as a veto as a part of a constrained maximization of signal-to-noise for times in which the chi-squared condition is satisfied. Thus, 2 m;n ½j needs to be computed for all j if it is computed at all. in time). 3 The effect of the maximizing over a chirp is to retain any true signal without introducing any significant bias in parameters, e.g., time of arrival (which can be demonstrated by simulations), while reducing the number of triggers that are produced by noise artifacts.

XI. EXECUTION OF THE FINDCHIRP ALGORITHM
In this section, we describe the sequence of operations that comprise the FINDCHIRP algorithm and highlight the tunable parameters of each operation. The algorithm is also illustrated in Fig. 1.
Since we are only describing the FINDCHIRP algorithm itself, we assume that a bank of templates ðM; Þ m has already been constructed for a given minimal match , The main computation takes place in two loops: the outer loop (in light grey) is over templates while the inner loop (in dark grey) is over data segments. The chi-squared branch, shown with a checkerboard pattern, is within the inner loop and is computationally costly. However, this branch is only executed if a candidate event is identified. Therefore, the computational cost is dominated by the procedure shown with a thick border: the calculation of the matched filter given by Eq. (8.10). according to the methods described in [8][9][10] and we are provided with either calibrated strain data s½j or the output error signal of the interferometer e½j and instrument calibration R½k as described in [48].
The initial operation is to scale the strain data by the dynamic-range factor , and divide it into data segments s n ½j suitable for analysis. The data segment duration T, stride length Á, and number of data segments N S in a block are selected. These quantities then define a data block length according to Eq. (5.2). A sample rate 1=Át must be chosen (which must be less than or equal to the sample rate of the detector data acquisition system); the sample rate and data segment length define the number of points in a data segment N ¼ T=Át. As mentioned previously, the lengths and sample rate are chosen so that N is an integer power of 2.
The next operation is construction of an average power spectrum 2 S½k using a specified data window w½j and power spectrum estimation method (Welch's method, the median method, or the median-mean method). The number of periodograms used in the average power spectrum estimate is typically chosen to be equal to the number of data segments, although different numbers of periodograms could be chosen. An inverse spectrum duration T spec is then given in order to construct the truncated inverse power spectrum À2 Q½k, according to Eq. (7.3).
Each input data segment is Fourier transformed to obtain s n ½k. The quantity F n ½k described in Eq. (8.6) can then be constructed. All frequency components of F n ½k below a specified low-frequency cutoff f low are set to zero, as are the DC and Nyquist components. The templateindependent normalization constants &½k high described in Eq. (8.9) are also computed at this point.
The algorithm now commences a loop over the N T templates in the bank, using the specified signal-to-noiseratio and modified chi-squared thresholds, ? and Ä ? , and the method of maximizing over triggers. For each template ðM; Þ m , the FINDCHIRP template G m ½k is computed, according to Eq. (8.7). The high-frequency cutoff k high;m for the template is obtained using Eq. (3.6) and used to select the correct value of & 2 ½k high;m for the template. The normalized signal-to-noise threshold ? is then computed for this template according to Eq. (8.13).
An inner loop over the FINDCHIRP data segments is then entered. For each FINDCHIRP segment F n ½k and FINDCHIRP template G m ½k the filter output m;n is computed according to Eq. (8.10). The trigger selection algorithm described in Sec. X is now used to determine if any triggers should be generated for this data segment and template, given the supplied thresholds and trigger maximization method. If necessary, the chisquared veto is computed at this stage, according to Eq. (9.4) and the threshold given in Eq. (9.5). If any triggers are generated, the template parameters ðM; Þ m are stored, along with the termination timet 0 , signal-to-noise ratio, effective distanceD eff , termination phase 0 , chi-squared veto parameters, and the normalization constant 2 m of the trigger. The triggers are generated and stored to disk for later stages of the analysis pipeline.
The segment index n is then incremented and the loop over the data segments continues. Once all N S data segments have been filtered against the template, the template index m is incremented and the loop over templates continues until all N T templates have been filtered against all N S data segments.

XII. CONCLUSION
Profiling of the inspiral search code based on the FINDCHIRP algorithm was performed on a 3 GHz Pentium 4 CPU with seven data segments of length 256 seconds. The data was read from disk, resampled from 16384 Hz to 4096 Hz and filtered against a bank containing 474 templates using the FFTW package [50] to perform the discrete Fourier transforms; the resulting 1255 triggers were written out to disk. Of the 2909 seconds of execution time, 1088 seconds were spent performing complex FFTs required by the matched filter, and 1600 seconds performing the chi-squared veto. Of the time taken to perform the chisquared veto, 1244 seconds are spent executing inverse FFTs. In total, 2300 seconds of the 2900 are spent doing FFTs, which means that the execution of the FINDCHIRP algorithm is FFT-dominated, as desired.
In practice, the FINDCHIRP algorithm is only a part of the search for gravitational waves from binary inspiral. An inspiral analysis pipeline typically includes data selection, template bank generation, trigger generation using FINDCHIRP, trigger coincidence tests between multiple detectors, vetoes based on instrumental behavior, coherent combination of the optimal filter output from multiple detectors, and finally manual candidate followups. Pipelines vary between specific analyses and the topology of the analysis pipelines has evolved due to the specific demands of particular observing runs [14][15][16]18,19,[21][22][23][24][25].
It is simple to modify the FINDCHIRP algorithm to use restricted post-Newtonian templates higher than secondorder by adding addition terms to the construction of the FINDCHIRP template phase in Eq. (8.5b). The phasing equations used in the LSC implementation of FINDCHIRP are given in Appendix A for completeness. Appendix A described the modification of FINDCHIRP used to search for time-domain templates, e.g., those based on post-Newtonian resummation techniques, such as the effective one body formalism [37]. These modifications cannot make use of the factorization used in the stationary phase templates, but they allow efficient reuse of the search code developed and tested for the frequency-domain post-Newtonian templates. Laboratory coperative agreement PHY-0107417. Patrick Brady is also grateful to the Alfred P Sloan Foundation for support. Patrick Brady and Duncan Brown also thank and the Research Corporation Cottrell Scholars Program for support.

APPENDIX A: ALGORITHM FOR TEMPLATES GENERATED IN THE TIME-DOMAIN
The optimization of the FINDCHIRP algorithm described above is dependent on the use of frequency-domainrestricted post-Newtonian waveforms as the template. It is a simple matter, however, to modify the algorithm (and hence the code used to implement the algorithm) so that an arbitrary waveform generated in the time-domain hðtÞ may be used as the matched-filter template. This allows use of inspiral templates such as the effective one body (EOB) [37] waveforms. These waveforms have been tuned to numerical relativity simulations of coalescing black holes and have a higher overlap with higher-mass signals in the sensitive band of the LIGO detectors [38]. Other waveforms of interest are described in [56][57][58][59][60]. In this Appendix, we describe the modifications necessary to use time-domain templates in FINDCHIRP. We assume that the desired template waveform is generated in the form where t 0 and 0 are the termination time and phase, as described in Sec. III, and A m ðtÞ and m ðtÞ are the particular amplitude and phase evolution for the mth template in the bank. The bank may include parameterization over binary component spins as well as masses. The template waveform is generated from the low-frequency cutoff f low and is normalized to the canonical distance of 1 Mpc. Recall the factorization of the matched-filter output, given by Eq. (8.3):z n;m ½k ¼ ðÁfÞ À1 A 1 Mpc;m F n ½kG m ½k: (A2) Since we are now only provided with the numerical value of the waveform as a function of time, we cannot perform the same factorization of the waveform as for stationary phase templates. Instead, to compute the FINDCHIRP data segment F n ½k, we remove the template-dependent amplitude by making the replacement k À7=6 ! 1 in Eq. (8.6). Similarly, the form of A 1 Mpc;m is now much simpler, with A 1 Mpc;m ¼ 1.
To construct the FINDCHIRP template G m ½k, we construct a segment of length N sample points and populate it with the discrete samples of the template waveform h 1 Mpc;m ½j. The waveform is sampled at the sampling interval of the matched filter Át. When we place the waveform in this segment, we must ensure that the termination of the waveform is placed at the sample point j ¼ 0, i.e., t 0 ¼ 0. For example, if the template is a post-Newtonian waveform generated in the time-domain [40], then it is typical to end the waveform generation at the time which the frequency evolution of the waveform ceases to be monotonically-increasing and not the time at which the gravitational-wave frequency goes to infinity. Thus the last nonzero sample point of the generated template may not correspond to the formal time of coalescence. If the algorithm generating the waveform provides the FINDCHIRP algorithm a formal coalescence time, the waveform is placed so that the coalescence time is at the last sample point in the segment. Any subsequent signal (such as merger and ringdown, in the case of EOB) is wrapped around to the start of the segment. If a formal coalescence time is not provided, the waveform is placed so that the last nonzero sample of the template is at the end of the segment. As discussed in Sec. III, this causes some ambiguity in the meaning of the recorded arrival time, but as long as this is consistent among all detectors in the network, this is not a significant effect. After placing the waveform in the segment, we construct the discrete forward Fourier transform of the waveform, as described by Eq. (2.3)

APPENDIX B: BIAS IN MEDIAN POWER SPECTRUM ESTIMATION
Here we compute the bias of the median of a set of periodograms relative to the mean of a set of periodograms. We assume that the periodograms are obtained from Gaussian noise. In this Appendix, let us focus on one frequency bin of the periodogram, and for brevity we adopt the symbol x for the power in the frequency bin, that is, we define x ' ¼ P ' ½k for ' ¼ 1; . . . ; n. (Here n is the number of periodograms being averaged. It is either N S or N S =2 depending on the choice of method.) Let fðxÞ be the probability distribution function for x. For Gaussian noise, fðxÞ ¼ À1 e Àx= where ¼ hxi ¼ is the population mean of x. So ¼ hP½ki. The population median x 1=2 is defined by Thus the bias of the population median is ¼ ln2.
The sample mean is unbiased compared to the population mean. The sample mean is x ' : The expected value of " x is so " x is an unbiased estimator of . The sample median, however, does have a bias. The sample median is: To compute the bias, we first need to obtain the probability distribution for the sample median. For simplicity, assume now that n is odd. The probability of the sample median having a value between x med and x med þ dx med is proportional to the probability of one of the samples having a value between x med and x med þ dx med times the probability that half of the remaining samples are larger than x med and the other half are smaller than x med . Thus, the probability distribution for x med is given by where m ¼ ðn À 1Þ=2 is half the number of remaining samples after one has been selected as the median. Here, QðxÞ is the upper-tail probability of x, i.e., the probability that a sample exceeds the value x where the second equality holds for the exponential distribution function corresponding to the power in Gaussian noise. The normalization factor C is a combinatoric factor which arises from the number of ways of selecting a particular sample as the median and then choosing half of the remaining points to be greater than this value. Thus it has the value of n (number of ways to select the median sample) times n À 1 ¼ 2m choose ðn À 1Þ=2 ¼ m (number of ways of choosing half the points to be larger) where Bðx; yÞ ¼ R 1 0 t xÀ1 ð1 À tÞ yÀ1 dt is the Euler beta function. This factor can also be obtained simply by normalizing the probability distribution gðx med Þ. To do so it is useful to make the substitution t ¼ Qðx med Þ so that dt ¼ Àfðx med Þdx med and Now we can compute the expected value for the median. Note that for the exponential probability distribution, hx med i ¼ À lnt. Thus x med ¼ ÀC R 1 0 t m ð1 À tÞ m lntdt. To evaluate the integral, start with the definition (above) of the Euler beta function. Differentiate it with regard to the first argument, and use the definition of the digamma function c ðzÞ ¼ À 0 ðzÞ=ÀðzÞ to obtain  [66] to reduce c . The bias factor is therefore for odd n. This result makes sense: As n ! 1 the series approaches ln2 which is the bias for the population median. However, for n ¼ 1, ¼ 1, so there is no bias (the median is equal to the mean for one sample!).

APPENDIX C: CHI-SQUARED STATISTIC FOR A MISMATCHED SIGNAL
For simplicity we write the chi-squared statistic [13] in the equivalent form [cf. Eq. (7.1)] For brevity we have dropped the indices n and m; the explicit dependence on j will also be dropped hereafter.
In this Appendix we further simplify the notation by adopting normalized templatesũ½k ¼h 1 Mpc ½k=. In terms of these templates we define the inner products as in [13] ðs; for the p different bands, which are chosen so that ðu; uÞ ' ¼ 1=p, and the inner product With this notation, the signal-to-noise ratio is given by 2 ¼ jðs; uÞj 2 and chi-squared statistic is jðs; uÞ ' À ðs; uÞ=pj 2 1=p jðs; uÞ ' j 2 : To see how the chi-squared statistic is affected by a strong signal (considerably larger than the noise) 4 , suppose that the detector output s½j consists of the gravitational waveform av½j where a specifies the amplitude of the gravitational-wave. Here v½j is also a normalized [in terms of the inner product of Eq. (C6)] gravitational waveform that is not exactly the same as u½j. The discrepancy between the two waveforms is given by the mismatch The mismatch is the fraction of the signal-to-noise ratio that is lost by filtering the true signal av½j with the template u½j compared to if the template v½j were used. The chi-squared statistic is jðv; uÞ ' j 2 ðv; vÞ ' ðu; uÞ ' where we have used the Schwarz inequality to obtain the second line and the normalization condition ðu; uÞ ' ¼ 1=p to obtain the third. The final approximation assumes that ( 1. Thus, the chi-squared statistic is offset by an amount that is bounded by twice the squared signal-tonoise ratio observed times the mismatch factor. There is no offset for a template that perfectly matched the signal waveform. It can be shown [13] that in the presence of a signal and Gaussian noise that 2 has a noncentral chi-squared distribution [62] with ¼ 2p À 2 degrees of freedom and a noncentral parameter & 2 2 (where now might possibly be slightly greater than 2 times the measured signal-to-noise ratio squared owing to the presence of the noise). This distribution has a mean value of þ and a variance of 2 þ 4. We see then that the modified chisquared statistic Ä of Eq. (9.5) has a mean of & 2 and a variance of & ð4 or 8Þ=ðp þ 2 Þ (4 when p ) 2 and 8 when p ( 2 ) for Gaussian noise. Thus we would expect to set a threshold on Ä of Ä ? $ a few.

APPENDIX D: DISTANCE CONVENTIONS
The effective distance, D eff , is a measure of the inverse amplitude of a gravitational-wave signal produced in a detector. For a binary that is optimally located (directly above or below a two-arm interferometric detector) and oriented (where the binary's orbit is viewed face-on), the effective distance is equal to the actual distance, D. For a source that is not optimally located or oriented, the signal is weaker than one from an optimally located and oriented source at the same actual distance, and therefore the effective distance is greater than the actual distance.
If a signal from a system having a total mass M, a reduced mass , and an effective distance D eff is present in a detector's data stream, so that it induces a strain hðtÞ on the detector, then the expected signal-to-noise ratio found for a matched-filter search is if the data is filtered with a template that has exactly the same waveform parameters as the signal. A common source of confusion is the use of the term ''signal-to-noise ratio'' in reference to the matched-filter output for a particular template, m ðtÞ, as well as in reference to the measure of the strength of a signal relative to a detector's noise, %. The latter should be called the expected signal-tonoise ratio, while the actual observed signal-to-noise ratio will differ owing to the presence of random noise as well as the possible mismatch between the true signal and the template used in filtering.
The effective distance is related to other distances that have become widely used in gravitational-wave physics. The horizon distance, D hor , is the effective distance of a particular source, a 1:4M -1:4M binary neutron star system, that would produce an expected signal-to-noise ratio % ¼ 8. That is, the horizon distance is defined by where % ¼ 8, M ¼ 2:8M , and ¼ 0:7M . Horizon distance is therefore a measure of the sensitivity of a detector: it depends on the integral of the inverse power spectral density with a frequency weighting of f À7=3 .
The sense-monitor range [63], R, is related to the horizon distance by R ¼ F D hor where F À1 ' 2:2648. 5 It is called sense-monitor range after the data monitor SENSEMONITOR [64] that is run in the LIGO control rooms. The sense-monitor range represents not the farthest distance to which a detector is sensitive, but, rather, the radius of a sphere that would contain as many sources as the number that would produce an expected signal-to-noise ratio % ¼ 8 in a detector under the assumption of a homogeneous distribution of sources having random orientations. Specifically, suppose that there are N hor sources homogeneously distributed in a sphere of radius D hor . Each of these sources has an effective distance D eff (which is always greater than or equal to its actual distance D). The number N det of these sources that are detectable (i.e., which would produce an expected signal-to-noise ratio % ¼ 8) is simply the number of the sources that have D eff D hor ; the quantity F 3 ¼ N det =N hor is the fraction of the universe within the horizon distance that is actually probed by the detector. (Note that F is not equal to 5 À1=2 , as it is often mistaken to be.) While the actual distance, D, is a physical parameter of a source, the effective distance D eff is related to the inverse of the strain amplitude seen in a detector. The horizon distance D hor and sense-monitor range R, on the other hand, are not related to an actual source but rather are measures of the sensitivity of the detector. The horizon distance is a representation of the maximum distance to which a detector can see, while the sense-monitor range describes the volume that may be surveyed by a detector.
Finally, Ref. [65] defines a distance quantity which is useful in the computation of rate limits on coalescing binaries [43]. If an inspiral signal observed in the detector does not terminate in the sensitive band of the detector (as is the case for M & 10M binaries in Initial LIGO), then it follows from Eq. (3.4b) that the amplitude of the waveform is proportional to M À5=6 D eff . For the purpose of computing the efficiency of a search, it is convenient to define the chirp distance where M BNS is the chirp mass of a binary with m 1 ¼ m 2 ¼ 1:4M .

APPENDIX E: EXTENSION TO HIGHER POST-NEWTONIAN ORDERS
Extending to the FINDCHIRP algorithm to use higher post-Newtonian orders is straighforward. At 3.5 post-Newtonian order, the expression for the chirp time given by Eq. (3.5a) becomes (see Eq. (3.8b) of Ref. [