We describe an efficient method of matched filtering over long (greater than
1 day) time baselines starting from Fourier transforms of short durations
(roughly 30 minutes) of the data stream. This method plays a crucial role in
the search algorithm developed by Schutz and Papa for the detection of
continuous gravitational waves from pulsars. Also, we discuss the computational
cost--saving approximations used in this method, and the resultant performance
of the search algorithm.Comment: 4 pages, text only, accepted for publication in the proceedings of
  the 3rd Amaldi conference on gravitational wave

Schutz, B. F.

Williams, P. R.

English

arXiv

Peter R. Williams

Crossref

An Efficient Matched Filtering
Algorithm for the Detection of
Continuous Gravitational Wave Signals
Peter R. Williams and Bernard F. Schutz
MPI for Gravitational Physics, Albert Einstein Institute, Am Miihlenberg 1,
D-14476 Golm, Germany
INTRODUCTION
Neutron stars are perhaps the most promising class of gravitational wave (GW)
sources, and searches for such GW signals is particularly suited to the charac-
teristics of the GEO600 detector (see Schutz "Getting Ready for GEO600 Data"
gr-qc9910033). However, the instantaneous GW frequency of such a source will
evolve due to both intrinsic spindown effects and Doppler modulations induced by
the motion of the Earth. Thus because of the large parameter space of likely signals,
directly implemented optimal matched filtering is not computationally feasible.
In response to this problem, Schutz and Papa have developed an alternative
strategy: the Hough-Hierarchical search algorithm (see Schutz and Papa "End-
to-End Algorithm for Hierarchical Area Searches for Long-Duration GW Sources
for GEO600" gr-qc9905018). In order to carry out a blind search over a range of
intrinsic GW frequencies, the following three stages must be calculated for each
point in the parameter space of sky positions and intrinsic spindown parameters:
Stage I: Calculate demodulated Fourier transforms (DeFTs) on an intermediate
time baseline (of order 1 day) by combining FFTs of short durations (ap-
proximately 30 minutes) of the time series data. In this context demodulated
means that if there is a source at the sky position in question, and with the
intrinsic spindown parameters in question, then all spindown and modulatory
effects will have been correctly removed from the DeFTs: all signal power
will be confined to one and the same frequency bin in each DeFT. This fre-
quency is the intrinsic frequency of the source measured at the start of the
observing time. It is expected that the total observing time will be of order 4
months, and thus roughly 120 of these DeFTs will be calculated for each point
in parameter space.
Stage II: In general source parameters will not coincide exactly with those
searched for, and residual frequency evolution and modulation will remain
in the DeFTs. Thus, the peak in power associated with a given source may
CP523, Gravitational Waves: Third Edoardo Amaldi Conference, edited by S. Meshkov
© 2000 American Institute of Physics 1-56396-944-0/007$ 17.00
473
Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jsp
change frequency bins from DeFT to DeFT. Because of the relatively small
time baseline of these DeFTs and the resultant poor signal-to-noise of any
expected continuous GW signal, this evolution will be not directly apparent
in the DeFTs, but can be recovered statistically using the Hough Transform
algorithm.
Stage III: Calculate DeFTs for candidate sources with the full frequency resolution
of the total observation time, by combining the intermediate baseline DeFTs
produced in stage I.
Thus, during stage II, regions of the parameter space in which it is statistically
unlikely that there are GW sources are eliminated from the search. Thereby, in
stage III, the most computationally expensive part of the algorithm, the long time
baseline DeFTs are calculated over only a very small fraction of parameter space
and over a very small range of frequencies.
In this paper we outline the methods used in the first and third stages of this
algorithm in constructing a longer time baseline DeFT from a number of shorter
time baseline FFTs or DeFTs.
THE METHOD
Consider a time series xa of total duration T, which has been divided into M
short time series, each having N data points. Then the DeFT for a signal with a
time independent amplitude and phase 2?r<I>a&(A) is
ATM-1 _ M-lAT-l
i,(\\ — \^ ™ ,o-27"$a&(A) _ V^ V^ -2m$ajb(\) /i\
"'ox y — / CL^ — / / U/^jC , \ )
a=0 a=0 j=0
where the time indices are related by Na + j = a, and 6 is a long time baseline
frequency index. In the following discussion Latin indices j, fc, / always sum over
TV, while Greek indices sum over M. Note that $a&(A) is dependent on a vector A
of parameters which characterize the signal one is searching for. In searching for
GW signals from neutron stars these will include intrinsic spindown parameters,
and the position of the source in the sky. If xak is the matrix formed by carrying
out Fourier transforms along the short time index j in zaj, then equation 1 can be
written as
M-lN-l [ t N-l
**(*)=££*.*££
M-l AT-l
y^ Q (b A) y^ x kP k(b A)
a=0 k=0
(2)
where the product (Ja(&, X)Pak(b, A) is defined by the terms in square brackets, and
Qa(b, A) contains all parts of the square brackets independent of the short time
index j and short frequency index k.
474
Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jsp
In equation 2 we have effectively re-written equation 1, a long time baseline
DeFT in the time domain, as a sum (a index) of short time baseline DeFTs in the
frequency domain (k index), where Qa(b^ X)Pak(b, A) are these frequency domain
filters. In the presence of stationary noise with a flat spectrum, equation 2 is the
optimal detector. However, through applying various approximations, the detector
can be made "acceptably sub-optimal", in the sense that only a small fraction of
power from a signal is lost in comparison to the optimal case, while achieving vast
savings in computational cost.
To illustrate these mathematical approximations it is instructive to discuss equa-
tion 2 for a specific case of $aj&(A): a linearly varying frequency model, i.e. in the
continuum limit $(t, /0, /o) = /o^ + /o^2, where /0 and /0 are the intrinsic frequency
and spindown of the source respectively, and t is time. In the case of an actual
search for GW signals from pulsars, 3>aj&(A) will not be so simple. However, this
model is sufficiently complex to effectively demonstrate all of the approximations
to be discussed here.
In discrete form, $(*, /0, /0) can be written as $aj/3/(7) = (/3+Ml)(Na+j}/NM+
^(Na + j)2/./V2M2, where the long time baseline frequency index b = J3 + Ml.
The chosen discretization of the spindown parameter /o = 7/T2 is not practically
appropriate. However, in an actual search, a grid of points in spindown parameter
space will be chosen to ensure an acceptable loss of power from unresolved signals.
Thus, in the following discussion, only searches for resolved /o parameters will be
considered.
Approximation 1: By Taylor expanding the model phase function $(£) about the
middle of each short duration time series (i.e. about j — N/2) and discarding
terms of order (j/N)2 = t2 and higher, in the limit TV -> oo the function
Pak(P-)1">l) is R-e Pak(/3,l,~f) = sine x and Im Pa^(/3,/,7) = (1 — cosx}/x. In
the phase model considered here x = — 2?r (/3/M + / + (2a + 1)7/M2 — k) and
Qa(0, /, 7) = exp {-2m(ap/M + a^/M2)}.
Approximation 2: Consider the case where the short time baseline is chosen such
that the instantaneous model frequency f ( t ) = <J>(£,/0,/o,/o, • • •) does not
move by more than one short time baseline frequency bin over the duration of a
short time baseline data set, i.e. in the model discussed here |/o|T/M < M/T.
Then for a given a, the function Pak(b, A) will be peaked in power about the
model frequency averaged over the duration of time associated with the ath
short data set, i.e. about x = 0 (the first three terms in the above definition
of x are the index of this average model frequency). Thus only a few terms
around this model frequency will contribute significantly to the summation
over k in equation 2.
Approximation 3: The semi-periodic nature of Pak(b, A) means that this function
can be efficiently evaluated from a look-up table of values containing the peri-
odic parts, and three further operations: to calculate one instance of Pak(^i A)
will require only 8 floating point operations.
475
Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jsp
Approximation 4: If one approximates the model frequency parameter (3 in the
calculation of Pak(/3J^} as a fixed value, for example with /3 = /?0, equation
2 can be calculated as an FFT, i.e.
M-l T nterm 1
*/3l(7) - E \Q'a(t) E *«*£.*(&,', 7) e - 2 ™M, (3)
«=o L fc J
where nterm relates to approximation 2, Pafc(/3,/,7) is defined above, and
for the phase model discussed here Q'a(y} = exp{—2ni (a27/M2 — 7/4M2)}.
Thus for values of (3 sufficiently near to /?o, the loss in power due to this approx-
imation will be small. To obtain ^(7) for other values of /?, the calculation
must be repeated using another /30.
RESULTS AND DISCUSSION
Numerical tests have shown that if one chooses 10% as an acceptable loss in
power in comparison to the optimal case, then NFFT — 8 a^d nterm = 16 are the
preferred parameter combination, if the short time baseline T/M is chosen such
that in the phase model discussed here |/o|T/M < M/T. If one decides that only
a 5% loss in optimal power is acceptable, then this can be achieved with the same
parameters, but choosing T/M such that \fQ\T/M < M/2T.
The computational cost of calculating one DeFT in stage I in floating point
operations is
«>
where B is the bandwidth of the search. This is comparable to the computational
cost of the corresponding steps in the Hierarchical Stack / Slide algorithm of Brady
and Creighton ("Searching for Periodic Sources with LIGO: Hierarchical Searches"
gr-qc9812014). The Hough-Hierarchical search algorithm also has a number of
computational advantages. To calculate a given bandwidth of a DeFT requires
only the FFT data from this bandwidth and an additional small overlap. Thus the
algorithm can be easily parallelized by distributing data and work by bandwidth;
and no communication between processors is required. Also, the complete three
stage algorithm can be arranged in such a way that once a bandwidth of FFT
data is read from disk by a processor, all computation required on this data can be
carried out while this data is held in memory, thus time spent reading data from
disk is a negligible fraction of the total computational time: each processor will
need to read roughly 40 Mb from disk once every two weeks. Furthermore, little
additional memory is required as workspace for stages I and III: less than 100 kb.
The GE0600 data analysis team are currently working on coding this algorithm
in a computationally optimal manner, as well as integrating this with the Hough
Transform part of the procedure.
476
Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jsp


An efficient matched filtering algorithm for the detection of continuous gravitational wave signals

We describe an efficient method of matched filtering over long (greater than 1 day) time baselines starting from Fourier transforms of short durations (roughly 30 minutes) of the data stream. This method plays a crucial role in the search algorithm developed by Schutz and Papa for the detection of continuous gravitational waves from pulsars. Also, we discuss the computational cost--saving approximations used in this method, and the resultant performance of the search algorithm

Williams, P.

Schutz, B.

MPG.PuRe

An efficient Matched Filtering Algorithm for the Detection of Continuous Gravitational Wave Signals

Schutz, Bernard F.

Williams, Peter R.

Max Planck Society eDoc Server

An Efficient Matched FilteringAlgorithm for the Detection ofContinuous Gravitational Wave SignalsPeter R. Williams and Bernard F. SchutzMPI for Gravitational Physics, Albert Einstein Institute, Am Miihlenberg 1,D-14476 Golm, GermanyINTRODUCTIONNeutron stars are perhaps the most promising class of gravitational wave (GW)sources, and searches for such GW signals is particularly suited to the charac-teristics of the GEO600 detector (see Schutz "Getting Ready for GEO600 Data"gr-qc9910033). However, the instantaneous GW frequency of such a source willevolve due to both intrinsic spindown effects and Doppler modulations induced bythe motion of the Earth. Thus because of the large parameter space of likely signals,directly implemented optimal matched filtering is not computationally feasible.In response to this problem, Schutz and Papa have developed an alternativestrategy: the Hough-Hierarchical search algorithm (see Schutz and Papa "End-to-End Algorithm for Hierarchical Area Searches for Long-Duration GW Sourcesfor GEO600" gr-qc9905018). In order to carry out a blind search over a range ofintrinsic GW frequencies, the following three stages must be calculated for eachpoint in the parameter space of sky positions and intrinsic spindown parameters:Stage I: Calculate demodulated Fourier transforms (DeFTs) on an intermediatetime baseline (of order 1 day) by combining FFTs of short durations (ap-proximately 30 minutes) of the time series data. In this context demodulatedmeans that if there is a source at the sky position in question, and with theintrinsic spindown parameters in question, then all spindown and modulatoryeffects will have been correctly removed from the DeFTs: all signal powerwill be confined to one and the same frequency bin in each DeFT. This fre-quency is the intrinsic frequency of the source measured at the start of theobserving time. It is expected that the total observing time will be of order 4months, and thus roughly 120 of these DeFTs will be calculated for each pointin parameter space.Stage II: In general source parameters will not coincide exactly with thosesearched for, and residual frequency evolution and modulation will remainin the DeFTs. Thus, the peak in power associated with a given source mayCP523, Gravitational Waves: Third Edoardo Amaldi Conference, edited by S. Meshkov© 2000 American Institute of Physics 1-56396-944-0/007$ 17.00473Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jspchange frequency bins from DeFT to DeFT. Because of the relatively smalltime baseline of these DeFTs and the resultant poor signal-to-noise of anyexpected continuous GW signal, this evolution will be not directly apparentin the DeFTs, but can be recovered statistically using the Hough Transformalgorithm.Stage III: Calculate DeFTs for candidate sources with the full frequency resolutionof the total observation time, by combining the intermediate baseline DeFTsproduced in stage I.Thus, during stage II, regions of the parameter space in which it is statisticallyunlikely that there are GW sources are eliminated from the search. Thereby, instage III, the most computationally expensive part of the algorithm, the long timebaseline DeFTs are calculated over only a very small fraction of parameter spaceand over a very small range of frequencies.In this paper we outline the methods used in the first and third stages of thisalgorithm in constructing a longer time baseline DeFT from a number of shortertime baseline FFTs or DeFTs.THE METHODConsider a time series xa of total duration T, which has been divided into Mshort time series, each having N data points. Then the DeFT for a signal with atime independent amplitude and phase 2?r<I>a&(A) isATM-1 _ M-lAT-li,(\\ — \^ ™ ,o-27"$a&(A) _ V^ V^ -2m$ajb(\) /i\"'ox y — / CL^ — / / U/^jC , \ )a=0 a=0 j=0where the time indices are related by Na + j = a, and 6 is a long time baselinefrequency index. In the following discussion Latin indices j, fc, / always sum overTV, while Greek indices sum over M. Note that $a&(A) is dependent on a vector Aof parameters which characterize the signal one is searching for. In searching forGW signals from neutron stars these will include intrinsic spindown parameters,and the position of the source in the sky. If xak is the matrix formed by carryingout Fourier transforms along the short time index j in zaj, then equation 1 can bewritten asM-lN-l [ t N-l**(*)=££*.*££M-l AT-ly^ Q (b A) y^ x kP k(b A)a=0 k=0(2)where the product (Ja(&, X)Pak(b, A) is defined by the terms in square brackets, andQa(b, A) contains all parts of the square brackets independent of the short timeindex j and short frequency index k.474Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jspIn equation 2 we have effectively re-written equation 1, a long time baselineDeFT in the time domain, as a sum (a index) of short time baseline DeFTs in thefrequency domain (k index), where Qa(b^ X)Pak(b, A) are these frequency domainfilters. In the presence of stationary noise with a flat spectrum, equation 2 is theoptimal detector. However, through applying various approximations, the detectorcan be made "acceptably sub-optimal", in the sense that only a small fraction ofpower from a signal is lost in comparison to the optimal case, while achieving vastsavings in computational cost.To illustrate these mathematical approximations it is instructive to discuss equa-tion 2 for a specific case of $aj&(A): a linearly varying frequency model, i.e. in thecontinuum limit $(t, /0, /o) = /o^ + /o^2, where /0 and /0 are the intrinsic frequencyand spindown of the source respectively, and t is time. In the case of an actualsearch for GW signals from pulsars, 3>aj&(A) will not be so simple. However, thismodel is sufficiently complex to effectively demonstrate all of the approximationsto be discussed here.In discrete form, $(*, /0, /0) can be written as $aj/3/(7) = (/3+Ml)(Na+j}/NM+^(Na + j)2/./V2M2, where the long time baseline frequency index b = J3 + Ml.The chosen discretization of the spindown parameter /o = 7/T2 is not practicallyappropriate. However, in an actual search, a grid of points in spindown parameterspace will be chosen to ensure an acceptable loss of power from unresolved signals.Thus, in the following discussion, only searches for resolved /o parameters will beconsidered.Approximation 1: By Taylor expanding the model phase function $(£) about themiddle of each short duration time series (i.e. about j — N/2) and discardingterms of order (j/N)2 = t2 and higher, in the limit TV -> oo the functionPak(P-)1">l) is R-e Pak(/3,l,~f) = sine x and Im Pa^(/3,/,7) = (1 — cosx}/x. Inthe phase model considered here x = — 2?r (/3/M + / + (2a + 1)7/M2 — k) andQa(0, /, 7) = exp {-2m(ap/M + a^/M2)}.Approximation 2: Consider the case where the short time baseline is chosen suchthat the instantaneous model frequency f ( t ) = <J>(£,/0,/o,/o, • • •) does notmove by more than one short time baseline frequency bin over the duration of ashort time baseline data set, i.e. in the model discussed here |/o|T/M < M/T.Then for a given a, the function Pak(b, A) will be peaked in power about themodel frequency averaged over the duration of time associated with the athshort data set, i.e. about x = 0 (the first three terms in the above definitionof x are the index of this average model frequency). Thus only a few termsaround this model frequency will contribute significantly to the summationover k in equation 2.Approximation 3: The semi-periodic nature of Pak(b, A) means that this functioncan be efficiently evaluated from a look-up table of values containing the peri-odic parts, and three further operations: to calculate one instance of Pak(^i A)will require only 8 floating point operations.475Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jspApproximation 4: If one approximates the model frequency parameter (3 in thecalculation of Pak(/3J^} as a fixed value, for example with /3 = /?0, equation2 can be calculated as an FFT, i.e.M-l T nterm 1*/3l(7) - E \Q'a(t) E *«*£.*(&,', 7) e - 2 ™M, (3)«=o L fc Jwhere nterm relates to approximation 2, Pafc(/3,/,7) is defined above, andfor the phase model discussed here Q'a(y} = exp{—2ni (a27/M2 — 7/4M2)}.Thus for values of (3 sufficiently near to /?o, the loss in power due to this approx-imation will be small. To obtain ^(7) for other values of /?, the calculationmust be repeated using another /30.RESULTS AND DISCUSSIONNumerical tests have shown that if one chooses 10% as an acceptable loss inpower in comparison to the optimal case, then NFFT — 8 a^d nterm = 16 are thepreferred parameter combination, if the short time baseline T/M is chosen suchthat in the phase model discussed here |/o|T/M < M/T. If one decides that onlya 5% loss in optimal power is acceptable, then this can be achieved with the sameparameters, but choosing T/M such that \fQ\T/M < M/2T.The computational cost of calculating one DeFT in stage I in floating pointoperations is«>where B is the bandwidth of the search. This is comparable to the computationalcost of the corresponding steps in the Hierarchical Stack / Slide algorithm of Bradyand Creighton ("Searching for Periodic Sources with LIGO: Hierarchical Searches"gr-qc9812014). The Hough-Hierarchical search algorithm also has a number ofcomputational advantages. To calculate a given bandwidth of a DeFT requiresonly the FFT data from this bandwidth and an additional small overlap. Thus thealgorithm can be easily parallelized by distributing data and work by bandwidth;and no communication between processors is required. Also, the complete threestage algorithm can be arranged in such a way that once a bandwidth of FFTdata is read from disk by a processor, all computation required on this data can becarried out while this data is held in memory, thus time spent reading data fromdisk is a negligible fraction of the total computational time: each processor willneed to read roughly 40 Mb from disk once every two weeks. Furthermore, littleadditional memory is required as workspace for stages I and III: less than 100 kb.The GE0600 data analysis team are currently working on coding this algorithmin a computationally optimal manner, as well as integrating this with the HoughTransform part of the procedure.476Downloaded 26 Feb 2008 to 194.94.224.254. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jsp

http://pubman.mpdl.mpg.de/pubman/item/escidoc:39965/component/escidoc:39966/AIP523_473.pdf

An efficient Matched Filtering Algorithm for the Detection of Continuous Gravitational Wave Signals

Abstract

Similar works

Full text

Available Versions

Crossref

MPG.PuRe

Crossref

Max Planck Society eDoc Server

Crossref

Crossref

Crossref

Crossref

Crossref

Crossref