Implementation of barycentric resampling for continuous wave searches in gravitational wave data

0.25 MB·English

by  P. Patel

Implementation of barycentric resampling for continuous wave searches in
gravitational wave data
Pinkesh Patel,1 Xavier Siemens,2 Rejean Dupuis, and Joseph Betzwieser1
1California Institute of Technology∗
2University of Wisconsin Milwaukee†
(Dated: January 26, 2010)
We describe an efficient implementation of a coherent statistic for continuous gravitational wave
searchesfromneutronstars. Thealgorithmworksbytransformingthedatatakenbyagravitational
wave detector from a moving Earth bound frame to one that sits at the Solar System barycenter.
Manypracticaldifficultiesariseintheimplementationofthisalgorithm,someofwhichhavenotbeen
discussed previously. These difficulties include constraints of small computer memory, discreteness
of the data, losses due to interpolation and gaps in real data. This implementation is considerably
0 more efficient than previous implementations of these kinds of searches on Laser Interferometer
1 Gravitational Wave(LIGO) detector data.
0
2 PACSnumbers:

n
a
I. INTRODUCTION which keep track of the phase of the gravitational wave
J
signals over long periods of time, and (2) semi-coherent
6
methods [18], which combine shorter periods of data
2 Rapidly rotating neutron stars are among the most
promising sources of continuous gravitational waves. without tracking the phase (for example, taking Fourier

] They can emit gravitational waves through a variety of transforms of short segments of data and then summing
c
the power).
q mechanisms, including unstable oscillation modes [1, 2]
- anddeformationsofthecrust[1,3–6]. Neutronstarscan When the sky location and phase evolution of a neu-
gr radiate powerful beams of radio waves from their mag- tron star are known, a coherent search for continuous
[ netic poles. If a neutron star’s magnetic poles are not gravitationalwavesisrelativelystraightforward[17]. As-
alignedwithitsrotationalaxis,thebeamssweepthrough suming that the noise in a gravitational wave detector
2 space,andiftheEarthlieswithinthesweepofthebeams, follows Gaussian statistics, in the presence of a signal,
v
the star is observed as a point source in space emitting the signal to noise recovered in a search increases with
5
5 bursts of radio waves. Such a neutron star is called a thesquarerootoftheamountofdatausedinthe search.
2 pulsar [7, 8]. Since the first discovery [9], around 2000 Thisisbecausethesignalamplitudegrowslinearlywhile
4 pulsars have been detected [10–12]. the noise follows a random walk. Thus, with enough
. data, it is possible to recover any continuous signal out
2 Duetomagneticdipoleradiationandgravitationalra-
of noisy data.
1 diation,therotationalfrequenciesofneutronstarsslowly
9 decrease in time. Other than this effect, gravitational If certain parameters of the signal (sky location, fre-
0 wavesfromisolatedrotatingneutronstarsareessentially quency,spindownsandbinaryparameters)arenotknown
: the search becomes much more involved. The reason is
v monochromatic in the rest frame of the star. The waves
that the number of points needed to cover the search
i are continuous and their frequency is determined by the
X parameter space (and ensure no signals are lost) grows
rotational frequency of the star. The motion of the de-
r tector as the earthrotates aboutits axis andaroundthe like a large power of the amount of data used [19]. This
a
makesthesensitivityofgravitationalwavesearchescom-
sun, however, modulates the phase as well as the ampli-
putationally bound: One cannot simply integrate arbi-
tude of the received signal. In order to recover the sig-
trary amounts of data to gain sensitivity because there
nalfrominterferometricdataoptimally,bothoftheseef-
is not enough computational power available to perform
fectsmustbetakenintoaccount. Detectinggravitational
the search. Thus, more efficient code and greater com-
wavesfromneutronstarscouldrevealinformationabout
puting power are highly desirable, since they translate
the strength of neutron star crusts and the equation of
into more data being analyzed and therefore an increase
state of the nuclear matter that makes up the star [6].
in the sensitivity of gravitational wave searches.
Continuousgravitationalwavesmayalsobeproducedby
other sources, such as cosmic strings [13, 14]. A promising method for blind searches involves ex-
ploitinglarge-scalecorrelationsinthe coherentdetection
There are a number of techniques available for contin-
statistic [23]. Another method that has been success-
uous wave searches. These techniques can be loosely di-
ful in these kinds of searches is the hierarchical scheme
vided into two categories: (1) coherentmethods [15, 17],
of incoherently combining coherent sets of data. Some
of the methods currently under use include the Hough
transform and stack-slide [18].
∗[email protected] In this paper we focus on an efficient implementation
†[email protected] of coherent techniques. The method we present here is
2
similar to several previous implementations [25–27], in is90◦). The functions a(t)andb(t)bothdepend ontime
that it uses fast Fourier transforms (FFTs) to calculate andlocationofsourceanddetector,butareindependent
the so-called F-statistic [15], the logarithm of the like- of the polarization angle ψ.
lihood function maximized over the intrinsic (and un- Inthedetectorframethephaseofagravitationalwave
known) parameters of the gravitational wave produced produced by an isolated neutron star can be written
by a neutron star, but with one very important differ- as [15]
ence. We resample the time domain data to the Solar
System barycenter before taking a FFT. This allows us s tk+1 2π s tk
Ψ(t)=Φ +2π f(k) + n ·r (t) f(k) ,
touseasingleFFTtocalculatethedetectionstatisticfor 0 0 (k+1)! c 0 d 0 k!
arbitrarily many frequencies and an arbitraryamount of Xk=0 kX=0
(4)
observationtime, while previousimplementationshavea
whereΦ isthephaseatthestarttimeoftheobservation,
maximumfrequencybandandobservationtimethatcan 0
f(k) is the kth derivative of the frequency, c is the speed
be calculated with a single FFT, which are determined 0
oflight,αandδaretherightascensionanddeclinationof
by losses due to phase mismatch. Another set of tech-
the source, n =n (α,δ) is the unit vector of the source
niques, described in [28, 33, 34], implement stroboscopic 0 0
intheSolarSystembarycenter(SSB)referenceframe,r
resampling methods described in [32]. The stroboscopic d
is the position vector of the detector in the same frame,
method requiresdata at full bandwidth, and is therefore
and s is the order of the expansion. Neglecting changes
not suitable for distributed computing applications such
inthepropermotionofthestar,thethirdterminEq.(4)
as Einstein@Home.
is a correction to the phase due to the detector motion
In Section II we review the signal properties and the
relative to the neutron star.
nearly-optimal coherent statistic that can be used to
We can define Φ(t)=Ψ(t)−Φ (t), as well as defining
extract continuous signals from interferometric gravita- 0
tional wave data. In Section III we discuss how to
implement the calculation of the nearly-optimal coher- Φ (t)=2π s f(k) tk+1 +2πn ·r (t) s f(k)tk (5)
ent statistic in a computationally efficient way in both s 0 (k+1)! c 0 d 0 k!
k=1 k=1
the time and frequency domains. In Section IV we de- X X
scribetheresultsofacomputercodeusingthisalgorithm and
on software injections of gravitational waves into gaus-
sian noise. Lastly, in Section V we address important n0·rd(t)
t = . (6)
m
technical issues to do with practical implementations of c
barycentric resampling.
Equations (5) and (6) let us write
II. PRELIMINARIES
Φ(t)=2πf[t+t (t;α,δ)]+Φ (t;f(k),α,δ), (7)
m s
In this section we closely follow the method of Jara-
which has the modulation due to the detector’s motion
nowski, Krolak, and Schutz [15] to provide the back-
around the SSB clearly separated from the modulation
ground on the signal and the detection statistic. Power-
due to the gravitational wave’s intrinsic frequency, al-
recycled Fabri-Perot Michelson interferometers such as
though not the derivatives of the frequency.
those used by the Laser Interferometer Gravitational
An almost optimal statistic for the detection of con-
Wave Observatory (LIGO) are sensitive to the strain
tinuous gravitational wave signals is called the F-
caused by gravitational waves passing through it. The
statistic [15, 16]. It is the logarithm of the likelihood
strain measured at a detector can be written as [15]
function maximized over the intrinsic and unknown sig-
nal parameters. The F-statistic is given by
h(t)=F+(t)h+(t)+F×(t)h×(t), (1)
4 B|F |2+A|F |2−2CR(F F∗)
wheretisthetimeinthedetectorframe,andh+ andh× F = a b a b . (8)
S (f)T D
are the “plus” and “cross” polarizations of gravitational h 0
wave. F+(t) and F×(t) are the beam-pattern functions where S (f) is the one-sided spectral density of the de-
of the interferometer and are given by h
tector’s noise at frequency f and T is the observation
0
time. A, B, C, and D are given by
F (t)=sinζ[a(t)cos2ψ+b(t)sin2ψ], (2)
+
A=(aka);B =(bkb);C =(akb);D =A·B−C2 (9)
and
with
F×(t)=sinζ[b(t)cos2ψ−a(t)sin2ψ], (3)
T0
whereψ isthepolarizationangleofthewaveandζ isthe 2 2
(xky)= x(t)y(t)dt. (10)
angle between detector arms (which in the case of LIGO T0 Z−2T0
3
F and F are integrals defined as A. Time Domain Analysis
a b
T0
F (f)= 2 a(t)x(t)e−iΦ(t)dt (11) TheF-statisticcanbecalculatedfromatimeseriesdi-
a
−T0 rectlybyfollowingthestepsoutlinedinSectionII. How-
Z 2
ever,due to the large amounts of data involved,it is im-
and
practicaltodothisfortheentiredataset. Onewaytoad-
T0
F (f)= 2 b(t)x(t)e−iΦ(t)dt. (12) dressthisproblemistodividethedataintoband-limited
b −T0 time series, making it possible to analyze one small sub-
Z 2 bandatatime. Timeseriesspanningdifferentfrequency
We define a new time variable called tb as follows: bands are then analyzed in parallel on a Beowulf cluster
or a distributed computing system. In this section we
t =t+t . (13)
b m
provide details on how this is accomplished in the time
Taking a derivative with respect to t on both sides of domain, and address some of the difficulties that arise.
Eq. (13), we get
dt dt 1. Heterodyning, low-pass filtering, and downsampling
b m
=1+ (14)
dt dt
From Eqs. (6) and (14), we get Let the output of the instrument be the time series
x(t), and its Fourier transform be
dt n ·v (t) ∞
m 0 d
dt = c (15) x˜(f)= x(t)e−2πiftdt. (18)
−∞
where v (t) is the velocity of the detector in the SSB Z
d
frame andthus n0·vd(t) is the Doppler shift ofthe source Ifwe considerthe Fouriertransformofthe complex time
with respect to thecdetector. For a detector located on series xh(t)=x(t)e2πifht,
Earth, the maximum Doppler shift experienced is of the
∞
order of 10−4. Using this fact and equation 14 we get x˜ (f) = x(t)e2πifhte−2πiftdt
h
δtb ≈δt. Z−∞
We can thus rewrite the Eqs. for F and F as ∞
a b = x(t)·e−2πi(f−fh)tdt
T0 −∞
Fa(f)= −T20 a(tb)x(tb)e−2πiftbeiΦs(tb)dtb, (16) = xZ˜(f −fh), (19)
Z 2
and it is obvious that multiplying the time series x(t) by
e2πifht has shifted all the frequencies in the time series
T0
F (f)= 2 b(t )x(t )e−2πiftbeiΦs(tb)dt (17) x(t) by fh. This procedureis referredto as complexhet-
b b b b
−T0 erodyning.
Z 2
IfjustasmallfrequencybandB ofdataaroundf isof
which are just the Fourier transforms of the resampled h
interest,low-passfilteringfollowedbydownsamplingcan
data and the detector response, multiplied by a phase
be used to reduce the bandwidth of the data appropri-
eiΦs(tb) [15]. Eqs. (16) and (17) can be efficiently eval-
ately. Specifically, if we wish to downsample by a factor
uated using FFTs. Details of the resampling procedure
D, the new Nyquist frequency of our time series will be
can be found in Sec. IIIA2.
given by
f B
Nyq,old
III. IMPLEMENTATION OF BARYCENTRIC fNyq,new = = . (20)
D 2
RESAMPLING
A simple but effective downsampling technique involves
Gravitationalwavedetectorscollectdataattherateof pickingeveryDth pointinthetimeseries. Toavoidalias-
about16-20kHz forspansoftime onthe orderofayear. ing effects however, prior to downsampling a low pass
This means that typical searches for gravitationalwaves filter must be applied to the data with a sharp fall-off
willinvolveontheorderofaterabyte(TB)ofdata. Com- around the new Nyquist frequency. The heterodyned,
puters currentlyhavememories ofa few gigabytes(GB), band-limited,downsampledcomplextimeserieswillhave
makingitnecessarytobreakupthedataintopiecesthat asamplingtime∆t=1/B. Forexample,supposeweare
can fit in the memory of a single computer. To analyze only interested in analyzing data between 990 Hz and
the full datasethundreds tothousands ofthese comput- 1 kHz. By multiplying the data with the phase factor
ers can then be used together in the form of a Beowulf e2π(995)it, data at 995 Hz moves to 0 Hz (DC), 990 Hz
cluster,ortenstohundredsofthousandswithdistributed moves to -5 Hz, and 1 kHz to +5 Hz (we have taken t
computing systems such as Einstein@Home [22]. to be measured in seconds). To avoid aliasing problems
4
whenwedownsample,welow-passfilterthedataat5Hz,
the newNyquistfrequency. Wecanthendownsampleby
picking one point out of every 100. The resulting com-
plex time series willbe sampledat10 Hz and containall
the information in the original time series between 990
Hz and 1 kHz.
2. Barycentric resampling and heterodyne correction
In this section we explain how to use the low band-
width heterodyned complex time series to compute the
F-statistic given by Eq. (8).
InthefollowingwewillworkonlywithF . Theproce-
a
dureforF iscompletelyanalogous. Itiseasiesttobegin
b
with the integral definition for F in Eq. (11) with the
a
phase explicitly written out, namely,
FIG. 1: Graphical description of theresampling procedure
T0
F (f)= 2 a(t)x(t)e−2πif(t+tm)eiΦs(t)dt, (21)
a
−T0
Z 2
and a similar expression holds for F :
b
and a similar expression holds for F . The heterodyned
b
version of F is N
a F (f −f )= b(tk)z(tk)e−2πiftkbeiΦs(tkb)dt , (28)
b h b b b
F (f −f )= T20 a(t)x(t)e−2πi(f−fh)(t+tm)eiΦs(t)dt. kX=1
a h
−T0 where tk is the kth datum in the time seriesas measured
Z 2 b
(22) in the barycentric frame and dt = tk+1−tk. The rela-
If we already have a complex heterodyned time series b b b
tionship between t and t can be written as
b
x (t) (heterodyned in the detector frame), we can use it
h
to absorbsome (but not all) of the heterodyne exponent tk =tk+t (tk;α,δ). (29)
b m
in Eq. (22) as follows:
This relationship between tk and tk can be used to cal-
x(t)e−2πi(f−fh)(t+tm) =x (t)e2πifhtme−2πif(t+tm). culate z(tk) from the time series z(btk). In practice, one
h b
(23) starts out with z(tk), i.e. data sampled regularly in the
detector frame. Then we calculate Tk(tk), which are de-
b
ThismeansthatratherthanEq.(22),weshouldevaluate tector times corresponding to regularly spaced samples
in the barycentric frame. These Tk(tk) are irregularly
b
T0 sampled in the detector frame, but since we have z(tk),
Fa(f −fh)= 2 a(t)z(t)e−2πif(t+tm)eiΦs(t)dt, (24) we can calculate z(Tk(tk)) by using interpolation. The
Z−2T0 interpolatedtimeserieszb(Tk(tk))isthez(tk)ofEqs.(27)
b b
and (28). A similar procedure may be used to calculate
where
thea(tk)froma(tk),andtheb(tk)fromb(tk). Thefactor
b b
z(t)=xh(t)e2πifhtm. (25) of eiΦs(tkb) in Eqs.(27) and(28) is calculated using equa-
tion (5). In this case, instead of calculating Φ (tk), we
s
At this point we have an expression which looks like calculate Φ (Tk(tk)), which is equivalent to calculating
s b
Eqs. (11) and (12), and we can write the integral over Φ (tk). Whileintheoryonehastocalculatethequantity
s b
t instead as an integral over tb: n0·rd(t) in equation (5), in practice this information is
already encoded in Tk(tk) as
T0 b
Fa(f −fh)= −T20 a(tb)z(tb)e−2πiftbeiΦs(tb)dtb, (26) n0·rd(t)=tm·c=(tkb −Tk(tkb))·c. (30)
Z 2
with a similar expression for F . With all the parts of Eqs. (27) and (28) in hand, we can
b
The discrete version of Eq. (26) for a time series with compute Fa(f −fh) and Fb(f −fh).
N points reads In summary, the procedure is the following:
1. Start with a heterodyned, band-limited, downsam-
N
Fa(f −fh)= a(tkb)z(tkb)e−2πiftkbeiΦs(tkb)dtb, (27) pledxh(tk)with tk regularlyspacedintime, inthe
frame of reference of the detector.
k=1
X
5
2. Correctthex (tk)fortheheterodyningdoneinthe 1. Dealing with non-stationary and colored data
h
detectorframebymultiplyingwithe2πifhtm topro-
duce the z(tk). To deal with non-stationarities,variations in the noise
floor from SFT to SFT, and colored data, we can nor-
3. The z(tk) correspond to data irregularly spaced in malize our SFT data to absorb the 1/S (f) term in the
h
the barycentricframe. Calculate Tk(tk), whichare definition of the F-statistic in Eq. (8). If X is the
b α,k
times in the detector frame corresponding to regu- kth frequencybin ofthe αth SFT, thenwe canredefine a
larly sampled solar system barycenter times. normalized data point Xˆ as
α,k
X
4. Usinginterpolation,calculatez(Tk(tk))fromz(tk), Xα,k −→Xˆα,k = α,k , (31)
b S
which is the z(tk) used in Eqs. (27) and (28). α,k
b
where S is an estimate of thepone-sided power spec-
α,k
5. Similarly, from a(tk) and b(tk) calculate a(tk) and tral density for the kth frequency bin of the αth SFT.
b
b(tk) respectively. Estimatorsusedfor this purpose should be robustin the
b
presence of spectral features in the data, such as a run-
ning median.
6. Using FFTs, evaluate Eqs. (27) and (28) to calcu-
late F (f −f ) and F (f −f ).
a h b h
2. Merging SFTs into long time-baseline Fourier
7. Use Eq. (8) to calculate the F-statistic. transforms
There are many practical difficulties that arise when
dealingwithSFTs. Oftencontiguouschunksofdatahave
B. Frequency Domain Analysis
to be divided up into multiple SFTs and it is necessary
to coherently combine them into one long time-baseline
In the previous section we describe a practical way of
SFT. This is done using the Dirichlet kernel, which is
calculating the F-statistic from time series data. How-
theequivalentofasincinterpolation(idealinterpolation)
ever, in practice the calculation is done in the frequency
done in the time domain. In order to keep the compu-
domain for a couple of reasons. One is that much of
tational cost down, the Dirichlet kernel is truncated at
the code written in the LIGO Scientific Collaboration’s
a finite number of points (usually around 16). This in-
(LSC)ContinuousWavesworkinggroupistailoredtoan
troduces a slight interpolation error, which cannot be
analysis performed in the frequency domain and hence
avoided without sacrificing a large amount of computa-
there exist many data processing and validation tools to
tional power.
process the data that are useful to this code. Another
Suppose we divide the data x(t) of length T into M
0
reasonisthatgravitationalwavedetectorsaresubjectto
short chunks of length T each with N points, so that
SFT
manysourcesofnoise,someofwhichchangedailyoreven
T = MT . The discrete Fourier transform (DFT) of
0 SFT
hourly, such as wind, microseism, earthquakes, anthro-
the data is
pogenic noise, etc. These change the noise floor of any
NM−1
analysis as a function of time. Working in the frequency
X = x e−2πilb/NM, (32)
domain is a natural way to deal with this problem. b l
l=0
We begin a frequency domainanalysisby taking short X
time-baseline Fourier transforms of the time domain where x = x(l∆t), ∆t is the sampling time, and b is
l
data, called short Fourier transforms (SFTs). When we a long time-baseline frequency index. We can write the
calculate the F-statistic, we divide by the noise in the Fourier transform in terms of two sums:
instrument at that frequency, as shownin Eq. (8). How- M−1N−1
ever,Eq.(8) assumesthe noise is stationary. To account X = x e−2πib(j+Nα)/NM , (33)
b α,j
for the non-stationarity of the noise we need to weight
α=0 j=0
by the noise over time, which is done on a per SFT ba- X X
sis. This normalization process is described in the next wherexα,j =x((j+Nα)∆t). We canexpressthe xα,j in
section. terms of an inverse DFT of a short chunk of data,
The computational cost of estimating the noise per N−1
1
SFTscaleswiththenumberofSFTsandthusforafixed x = X e2πijk/N, (34)
α,j α,k
observation time scales inversely with the time-baseline. N
k=0
X
Acompromiseisneededbetweenthedemandsofcompu-
where the X are the starting SFT data,
tationaltime andrelativestationarityofthe detector for α,k
a given time-baseline. In LIGO, SFTs are usually 1800 N−1
seconds long, since the detector is reasonably stationary X = x e−2πijk/N. (35)
α,k α,j
for that time.
j=0
X
6
Replacing x with Eq. (34) in Eq. (33) gives and take a sub-band of Xˆ , inverse Fourier transform it,
α,j b
and produce the heterodyned time series, and correct it
M−1N−1 N−1
X = 1 X e2πijk/N e−2πib(j+Nα)/NM to produce zˆ(tkb). In terms of this time series, we can
b N α,k write
!
α=0 j=0 k=0
X X X
M−1 N−1 N−1 N
= 1 e−2πibα/M Xα,k e−2πij(b/M−k)/N . Fˆa(f −fh)= zˆ(tkb)a(tkb)e−2πiftkbeiΦs(tkb) (44)
N
αX=0 kX=0 Xj=0 Xk=1
(36)
and
The last sum in this expressioncan be evaluated analyt-
N
ically. In particular, Fˆ (f −f )= zˆ(tk)b(tk)e−2πiftkbeiΦs(tkb) , (45)
b h b b
N−1 1−zNc kX=1
zcj = . (37)
1−zc and thus
j=0
X
We take z = e, c = −iy/N, with y = 2π(b/M −k), so F = 4 B|Fˆa|2+A|Fˆb|2−2Cℜ(FˆaFˆb∗). (46)
that the sum is given by T D
0
N−1 1−e−iy
e−iyj/N = . (38)
1−e−iy/N 4. Heterodyning
j=0
X
InthelargeN limittheexponentofthedenominatorwill As shown before in Eqs. (18) and (19), heterodyning
be small so that is a procedure by which the frequency of interest can be
1−e−iy 1−e−iy iN shifted arbitrarily. When one applies the kind of correc-
≈ = (e−iy−1) tion in Eq. (18), we effectively move all the frequencies
1−e−iy/N 1−(1−iy/N) y
byasetamount. By doingso,weconvertthe time series
siny 1−cosy
= N( −i ). (39) fromarealtime seriesto acomplextime series,withthe
y y same amount of information content.
Heterodyning in the frequency domain can be done in
This means we can write Eq. (36) as
twoways,oneinwhichthetime seriesproducedafterin-
M−1 N−1 verseFouriertransformingisrealandanotherinwhichit
Xb = e−2πibα/M Xα,kPb,k, (40) iscomplex. Acosinetransformusedtoheterodynewould
α=0 k=0 producearealtimeseries,butthismethodisnotusedin
X X
with the Dirichlet kernel an implementation of the technique (see section IV). A
complex heterodyned time series is produced by inverse
siny 1−cosy
P = −i , (41) Fourier transforming a relabelled band of the frequen-
b,k
y y
cies. Since in Eq. (18), all frequencies are shifted by a
andy =2π(b/M−k). The functionP is verystrongly fixed amount, the equivalent procedure in the frequency
b,k
peaked around y = 0, which is near a value of the fre- domainisjustrelabellingtheheterodynefrequencyfh as
quency index k∗ = floor(b/M). This means one only DCandsubsequentlyalltheotherfrequenciesrelativeto
needs to evaluate the sum over k for a few terms ∆k this new DC.
around k∗. With this in mind we write Taking the example from Section IIIA1, we can just
internally change the labels of the 995 Hz frequency bin
M−1 k∗+∆k
X ≈ e−2πibα/M X P . (42) to DC and 1000 Hz to 5 Hz. Once this relabelling is
b α,k α,k done, the original data will have all shifted by 995 Hz,
α=0 k=k∗−∆k
X X with the 10 Hz from -5 Hz to +5 Hz containing all the
To produce a heterodyned time series a sub-band of the relevant information. If one were using the whole band
Xb may be selected and inverse Fourier transformed. without downsampling or filtering, then this relabelling
would have to wrap around the Nyquist frequency edge,
but since the whole purpose of heterodyning is to down-
3. Normalized long time-baseline Fourier transforms sample, it is never necessary to do so.
With the normalized SFT data Xˆ from Eq. (31)
α,k
we can construct a normalized version of the long time- 5. Downsampling and low-pass filtering
baseline Fourier transform
M−1 k∗+∆k Following the time domain algorithm, after hetero-
Xˆ ≈ e−2πibα/M Xˆ P , (43) dyning the data, it needs to be downsampled and low-
b α,k α,k
α=0 k=k∗−∆k passfiltered. Thedownsamplingandlow-passfilteringis
X X
7
achieved by simply throwing out the data that is not in
the band of interest. The heterodyning is done in such a
wayastokeepthecenterofthebandofinterestatDC.A
Tukey window applied to the band of interest, keeping a
littlebitofdataonbothedgestofacilitatetheriseofthe
window from 0 to 1, is a good choice of a low-passfilter.
Once an inverse Fourier transform is performed on this
smaller subset in the frequency domain, it generates the
same heterodyned, downsampled, and low-pass filtered
time series as the time domain algorithm.
FIG. 2: Pictorial description of data pre-processing
6. Gaps in the data
Datacollectedbyaninterferometerwillhavegapsdue
to periods of downtime. These gaps need to be dealt 2. Normalize these SFTs and assign them weights.
with in a manner that preserves the phase coherence of
the segments around the gaps. The gaps increase the 3. Identify contiguous sets of SFTs.
analysis time without contributing any power to the F-
4. Combine each contiguous chunk of SFTs into one
statistic, and thus act like a zero padding.
long time-baseline Fourier transform (FT).
The data is divided up into a series of contiguous
chunksandgaps. ForeachcontiguouschunktheSFTsin
5. Createadownsampled,heterodyned,andlow-pass-
that chunk are normalized, patched up and then a het-
filtered time series by inverse Fourier transforming
erodyned,downsampledandlow-pass-filteredtime series
the desired frequencies from the FT.
is calculated from it. Heterodyning done by relabelling
is equivalent to multiplying with e2πifh(t−tc), where tc 6. Stitchallthesetimedomainchunkstogether,filling
is the start time of the data chunk being heterodyned gaps with zeros.
and f is the heterodyne frequency. If we have multiple
h
chunks that are separately being heterodyned, then t is
c
differentforeachchunk. Inthetimedomainanalysis,we
assumed that the heterodyne reference time is the same IV. RESULTS
as the start time of the analysis. In order to achieve the
same kind of heterodyning, one needs to multiply each A. Speed
newly createdtime series with a correcting phase factor,
namely
The scheme previously used to compute the F-
e2πifh(tc−ts), (47) statistic, involvedthe use of the Dirichlet kernelto com-
bine a series of SFTs [31] [24], which were calcu-
where t is the start time of the overall analysis.
s lated for 30 minutes of data taken at 16 kHz. The
A Tukey window can then be applied to each of these
30 minute window was set by the maximum Doppler
time series to smoothly bring the data to zero at the
shift due to the motion of the Earth. A C code
edges, which correspond to the gaps. The gaps are then
called ComputeFStatistic v2 [29] was written in the
filled with zeros, as no data was collected during those
LIGO Analysis Library (LAL) to calculate the F-
times. This procedure is repeated for all the gaps and
statistic using this algorithm. The code which imple-
contiguouschunks. Attheend,atimeseriesisproduced,
ments our method is also written in C and is called
which is contiguous and spans the time of the analysis.
ComputeFStatistic resamp[30]. Henceforthwewillrefer
By ensuring that the timestamps of the first datum of
to the previous implementation as the LAL implementa-
each contiguous chunk correspondwith the start time of
tion and our implementation as Resampling.
that chunk, we ensure that the phase coherence is main-
The F-statistic is calculated for a series of templates
tained throughtout.
looping over various parameters such as sky location, α
andδ,spin-downsfk,andvariousfrequenciesf. We can
ignore the way the two implementations deal with loops
7. Summary
over α, δ, and fk, since they both loop over them in the
same manner. The speed of computation for a loop over
To summarize, a simple algorithm to produce a time
frequencies f is worth comparing, however.
series equivalent to the one used for the time domain
Assume that we have N data points (take for example
analysis is as follows:
106 secondsofdataat100Hz,i.e. 108datapoints). Now
1. DividethedataintotimechunksandFouriertrans- assume that the number of operations per sky location
form them to create SFTs. and per spin-down is N . If the number of Dirichlet
ops
8
kernel points used is N , then the total number of Histogram of F−Stat values at injected frequecy
DirKer 0.03
operations used by the LAL implementation is: Resampling
n LAL Implementation
o0.025
cti Theoretical Prediction
n
u
F 0.02
NTLoAtL =Nops·NDirKer·NSFTs·N , (48) ensity 0.015
D
y
Where Nops is defined as the number of operations con- abilit 0.01
ducted in the innermost loop and is approximately of Prob0.005
order 10, N is the number of times the Dirichlet
DirKer
Kernel loop is repeated, N = Tobs is the number of 00 20 40 60 80 100 120 140
SFTs TSFT 2F
SFTs, and N is the number of data points.
Compare this to the resampling method, which con- FIG. 3: Histogram of results of Monte Carlo simulation with
sists of 4 major steps: signals injected in different instances of noise
1. Calculating t (t), given a sky location and time.
b
2. Calculating the integrands of F and F .
a b
3. Interpolating and calculating the beam patterns. B. Validations
4. Taking the Fourier transform.
The probability density distribution of the F-statistic
Each of these steps involves order 10 operations, but for Gaussian noise of zero mean and unity standard de-
all of these steps are sequential, therefore they only viationis a χ2 distribution with four degreesof freedom.
add, resulting in a total number of operations per data In the presence of a signal, the distribution is a χ2 of
point,NResamp,ofapproximately30operations. Thelast four degrees of freedom with a non-centrality parameter
ops
step is the Fourier transform, which is of order NlogN, given by the F-statistic in the absence of noise for the
therefore the total number of steps is: particular signal.
NResamp =(NResamp+logN)·N . (49) Resampling uses various approximate methods in the
Tot ops
calculation of the F-statistic, and this can lead to dis-
Therefore the ratio of operations between the two meth- agreements between the theoretical F-statistic probabil-
ods is ity density function and the output of the code. These
changes are of the order of a few percent and are within
NLAL N ·N ·N acceptable limits. The validity of the code can be tested
Tot = ops DirKer SFTs (50)
NResamp NResamp+logN byusingaMonteCarlosimulationofaboutamilliondif-
Tot ops ferent injections of the same signal in different instances
of noise. The noise is generated as a Gaussian noise of
To first order, we have
zero mean and unity standard deviation, and the signal
is added into this noise. For each individual injection
NLAL N
Tot ≈ SFTs . (51) the signal is chosen with a given set of amplitude pa-
NTRoetsamp logN rameters and a fixed sky location and spindowns, and
the search is conducted over these exact chosen param-
Therefore for large observation times, this method of eters in order to avoid any mismatches. These Monte
calculating the F-Statistic is faster and, in the case of a Carlosarethenrepeatedwithanothersetofparameters,
targeted search, it allows for a large parameter space in which are themselves chosen randomly. While it is not
F(k)’s. an exhaustive test, randomly chosen parameters ensure
The speed-up in practice is reduced by a few practi- that we are not biased in the validation test. The plot
cal issues as seen in section V. However, Resampling in figure 3 is produced by performing one such Monte
is still considerably more efficient than the LAL imple- Carlo simulation. In this case, both the LAL implemen-
mentation. ForEinstein@Home,becauseoftherelatively tationandResamplingwererunonthe samesetofdata.
small coherent integration time, the speed-up is around The F-statistic was picked out at the appropriate fre-
10. But for targeted searchesthat span multiple months quency and this was repeated about a million times. A
or years, the improvement can be as high as a factor of histogram of these F-statistic values was then plotted.
2000. Thus,whilesometargetedsearcheswhichintegrate As one can see, there is very good agreementin between
over a couple of years were impossible to do previously, the expected distribution of the F-statistic and the two
they are now possible. implementations.
9
V. PRACTICAL CONSIDERATIONS to discard the higher frequencies.
VI. SUMMARY AND CONCLUSIONS
A. Discreteness
In this paper, we describe an efficient implementation
In the implementation of the algorithm explained ofthebarycentricresamplingtechnique,whichdealswith
above, one major obstacle is the fact that the data col- the non-stationarity of the detector and calculates the
lected by any physical instrument is discrete and thus F-statistic. Although the calculation of the F-statistic
must be handled appropriately. Take, for example, the has been targeted, this technique can be used for many
heterodyne frequency used in the calculation. This fre- other kinds of searches. The major contribution of this
quency cannot be chosen arbitrarily, as only certain fre- technique is to remove the Doppler shift of the Earth’s
quenciesaresampledandthusthereareonlycertainper- motion in a gravitational wave signal. Thus, once this
mitted choices. Doppler shift is removed, both frequentist and Bayesian
Most major FFT computation algorithms output the techniques can be applied to the data. In the process of
frequencyseriesinaspecificformat,whichsplitthe data implementing this algorithm, a series of practical issues
into two parts. The first bin output by these algorithms aredealtwith,includingconstraintsofmoderncomputer
is the DC followed by the first positive frequency bin memory, discreteness of the data taken, losses due to in-
up to positive Nyquist and then follows this up with terpolation, and gaps in real data.
the negativefrequenciesstartingatthe negativeNyquist The computational savings due to this technique can
frequency. This order of placing frequency bins speeds be used in various ways. One such use is to increase
up computation and is necessary for the internal work- thecoherentintegrationtime forall-skysearcheslikethe
ings of these algorithms. Thus when an inverse FFT is Einstein@Homesearches. CurrentlyEinstein@Home[22]
performed on the frequency domain data in the form of uses 40 hour long coherent integrationtime. The resam-
SFTs, a simple reshuffling needs to be done. The fre- pling code will be about 10times faster for such integra-
quency selected to be the first bin will become the new tion times, and for the same computational power and
DC and thus the data will have been heterodyned by keeping the same scaling for the search, we can coher-
that said frequency. In order to ensure that the same ently integrate 64 hours instead, which corresponds to a
frequency bin is chosen as DC, one needs an odd num- sensitivity increase of about 25%.
berofbinsperSFT.Ifthenumberofbinsareeven,then The resampling technique is mosteffective for long in-
uponincreasingtheamountofdataitcanshiftthisnum- tegration times, which are feasible for targeted searches
ber to an odd number as the increase is always done by likethesearchforgravitationalwavesfromtheCrabpul-
changingthe number ofSFTs. Butifthe number ofbins sar[21]. Thecomputationalsavingscanbeusedtosearch
perSFTisodd,thenitwillremainoddforanynumberof over wider parameter spaces like more spindown param-
SFTs. Thisensuresthatthereisnomismatchinchoosing eters or to search over binary systems.
the appropriate bin as the heterodyne frequency.
B. Interpolation Issue Acknowledgments
When the resampling algorithm is used on discrete We would like to thank the membership of the LIGO
data, one needs to interpolate between data points to VirgoCollaboration’scontinuouswavesgroup,especially
gofromthedetectorframetotheSSBframe. Thisinter- Stefano Braccini, Vladimir Dergachev, Greg Mendell,
polationactslikeanonlinearlow-passfilteranddestroys Chris Messenger, Marialessandra Papa and Reinhard
the power at higher frequencies in the band of analysis. Prix for helpful discussions. We are also grateful to
Sincethefilterisnonlinear,thefrequencyresponseisnot Patrick Brady and Alan Weinstein for useful conversa-
well-defined, and thus there is no way to compensate for tions. LIGO was constructed by the California Insti-
the power loss at high frequencies. The power loss can tute of Technology and Massachusetts Institute of Tech-
besignificant(oforder30%)andisunacceptableinmost nology with funding from the National Science Founda-
analyses. The exact nature of the filter depends on the tion and operates under cooperative agreement PHY-
type of the interpolation routine used and the sky loca- 0757058. This paper has LIGO Document Number
tion that one resamples to. The only work around is to ligo-p0900301-v1. XavierSiemensissupportedinpartby
performthecomputationoveralargerbandthantheone NSF GrantNo. PHY-0758155andthe ResearchGrowth
desired. Inpracticeitissufficienttodoublethebandand Initiative at the University of Wisconsin-Milwaukee.
[1] L. Bildsten, Astrophys. J. 501 L89, 1998. trophys. J. 516 307, 1999.
[2] N. Andersson, K. D. Kokkotas, and N. Stergioulas, As- [3] G. Ushomirsky, C. Cutler, and L. Bildsten, Mon. Not.
10
Roy. Astron. Soc. 319 902, 2000. [21] B.Abbottetal.[TheLIGOScientificCollaboration],As-
[4] C. Cutler, Phys. Rev. D66 084025, 2002. trophys. J. 683, L45 (2008)
[5] A. Melatos and D. J. B. Payne, Astrophys.J. 623 1044, [22] B. Abbott et al. [LIGO Scientific Collaboration], Phys.
2005. Rev. D79, 022001 (2009)
[6] B. J. Owen, Phys. Rev. Lett. 95 211101, 2005. [23] H. Pletsch and B. Allen Phys. Rev. Lett. 103, 181102
[7] T. Gold, Nature 218 731, 1968. (2009)
[8] F. Pacini, Nature 219 145, 1968. [24] P.R.WilliamsandB.F.SchutzarXiv:gr-qc,9912029v1
[9] A. Hewish, S. J. Bell, J. D. H. Pilkington, P. F. Scott, [25] P. Astone, K.M. Borkowski, P. Jaranowski, and
and R.A. Collins, Nature 217 709, 1968. A. Kr´olak Phys. Rev. D65, 042003 (2002)
[10] R.N.Manchester,G.B.Hobbs,A.Teoh,andM.Hobbs, [26] P. Astone et al. Class. Quant. Grav. 20, s665 (2003)
TheATNFPulsarCatalogue.http://www.atnf.csiro.au/ [27] A.Kr´olak, M. Tinto, and M. Vallisneri Phys. Rev. D70,
research/pulsar/psrcat/. 022003 (2004)
[11] R.N.Manchester,G.B.Hobbs,A.Teoh,andM.Hobbs, [28] P. Astone et al. Class. Quant. Grav. 25, 184012 (2008)
The ATNFPulsar Catalogue 2004. [29] ComputeFStatistic v2, LIGO Analysis Library(LAL),
[12] D. Lorimer, Binary and Millisecond Pulsars. http:// https://www.lsc-group.phys.uwm.edu/daswg/projects/
relativity.livingreviews.org/Articles/lrr-2005-7/. lal.htm
[13] F. Dubath and J. V. Rocha, Phys. Rev. D76, 024001 [30] ComputeFStatistic resamp, LAL, https://www.lsc-
(2007) group.phys.uwm.edu/daswg/projects/lal.htm
[14] M. R. DePies and C. J. Hogan, arXiv:0904.1052 [astro- [31] R. Prix, LIGO Document T0900149-v1,
ph.CO]. https://dcc.ligo.org/cgi-bin/DocDB/ShowDocument?
[15] P. Jaranowski, A. Krolak, and B. F. Schutz, Phys. Rev. docid=1665
D58063001, 1998. [32] B. F. Schutz, in The Detection of Gravitational Waves,
[16] R. Prix,B. Krishnan Class. Quant. Grav. 26, 204013 editedbyD.G.Blair(CambridgeUniversityPress,Cam-
(2009) bridge, 1991), Chap.16, pp. 406-451.
[17] B. Abbott et al., [The LIGO Scientific Collaboration], [33] O. Torre, Thesis, University of Pisa(2009), Published as
Phys. Rev. D69 082004, 2004. LSC-Virgo Technical Report(VIR-0013A-10)
[18] B. Abbott et al. [The LIGO Scientific Collaboration], [34] S. Braccini et al. ”Geometrical signal resampling in
Phys. Rev. D77, 022001 (2008) continuous gravitational wave search to compensate for
[19] P. R. Brady, T. Creighton, C. Cutler and B. F. Schutz, Doppler and Spin Down effect”, to be submitted to
Phys. Rev. D57, 2101 (1998) Class.Quantum.Grav.
[20] B. Abbott et al., [The LIGO Scientific Collaboration],
Phys. Rev. D76 082001, 2007.