Dirty details of optimal filtering: wraparound and windowing

To carry out optimal filtering, we need to break the data set (which might be hour, days, or weeks in length) into shorter stretches of $N$ points (which might be seconds or minutes in length). We can understand the effects of ``chopping up" the data most easily in the case for which (1) the instrument noise is white, so that $S_h(f)=1$; (2) the source is so close that its signal overwhelms the noise in the IFO, and (3) we are looking for a signal with a given phase (not a linear combination of the two orbital phases).

We want to calculate a signal $S$ as a function of lag $t_0$ using an FFT.

$$S(t_0) = \int h(t) T(t-t_0) dt \approx S(i_0) = \sum_j h_j T_{j-i_0},$$ (6.19.100)

where we have written both the continuous-time and discrete-time version of the same equation. Using the definition of the discrete Fourier transform, and writing

$$h_j = \sum_{k=0}^{N-1} {\rm e}^{-2 \pi i j k/N} \tilde h_k \... ...'=0}^{N-1} {\rm e}^{-2 \pi i (j-i_0) k'/N} \tilde T_{k'} \quad$$ (6.19.101)

one can easily compute that the signal as a function of lag $i_0$ is

$$S(i_0) = \sum_{j=0}^{N-1} \sum_{k=0}^{N-1} \sum_{k'=0}^{N-1} {\rm e}^{-2 \pi i j k/N} \tilde h_k {\rm e}^{-2 \pi i (j-i_0) k'/N} \tilde T_{k'}$$ (6.19.102)

$$= \sum_{k=0}^{N-1} \sum_{k'=0}^{N-1} N \delta_{k,-k'} {\rm e}^{2 \pi i i_0 k'/N} \tilde h_k \tilde T_{k'}$$ (6.19.103)

$$= \sum_{k=0}^{N-1} N {\rm e}^{-2 \pi i i_0 k/N} \tilde h_k \tilde T^*_{k}.$$ (6.19.104)

Thus, if the data is treated as periodic, and the template is treated as periodic, one can compute the correlation as a function of time using only an FFT. In particular, the use of rectangular windowing does create sidelobes of the template's frequency components. However it also creates identical sidelobes of the signal's frequency components - so in effect the correlation in the time domain can be calculated exactly, without any windowing of the signal being necessary.

The only complication arises from the fact that the FFT treats the data as being periodic. Let's consider some simple examples to illustrate the effects of this. In all of our examples, the number of data points is $N=65,536=2^{16}$ and the (schematic) chirp filter has length $m=13,500$ and is zero-padded after that time. Please remember, in all the figures that follow, to identify the far right hand side of the graph ($i=65535$) with the far left hand side ($i=0$). Figure  shows $S(i_0)$ for a schematic chirp which begins at the first data point in the rectangular window. You will notice that the filter output peaks at $i=0$.

Figure: A chirp starting at initial time $i=0$ and ending at time $i=13500$ is processed through a chirp filter, whose output peaks at time $i=0$. Notice that because of wraparound, the (non-causal) filter output begins ``earlier" than $i=0$.

If the incoming chirp arrives somewhat later (it starts at $i=15,000$) as shown in Figure  then the filter output peaks at the start time, as shown.

Figure: A chirp starting at initial time $i=15,000$ and ending at time $i=28,500$ is processed through a chirp filter, whose output peaks at time $i=15,000$.

A chirp in the signal which starts at the $i=65,535-13,500$ as shown in Figure  causes the filter output to peak at $i=52,035$.

Figure: A chirp starting at initial time $i=52,035$.

Thus, in order to find chirps, we need to find the maxima of the filter output over the interval $i=0\cdots,N-m$.

Chirp filters can be ``stimulated" or ``triggered" by events that are not chirps. We will shortly discuss some techniques that can be used to distinguish triggering events that are chirps from those that are simply noise spikes or other transient (but non-chirp) varieties of non-stationary interferometer noise. Suppose that a chirp filter is triggered by some kind of transient event in the IFO output. At what time did this transient event ocurr? The answer to this question can be seen by examining the impulse response of the ``periodic filter" scheme, as shown in the following figures.

Figure: An impulse shortly after $i=0$.

Figure: An impulse at $i=15,000$.

Figure: An impulse at $i=28,500$.

Figure: An impulse shortly before $i=65535$

Thus, by searching for maxima in the filter output over the range $i=0,\cdots,N-m-1$ we can detect either true chirps in the data stream, starting in the time interval $i=0,\cdots,N-m-1$ and coalescing (roughly speaking) in the time interval $i=m,\cdots, N-1$, or we can detect transient impulse-like events in the data stream, which take place in the time interval $i=m,\cdots, N-1$. In the GRASP optimal filtering code, after examining the stretch of $N$ data points, we then shift the data points $i=N-m,\cdots,N-1$ into the range $i=0,\cdots,m-1$ and acquire a new additional set of $N-m$ data points covering remaining (new) time interval.

To indicate the time at which the filter output reached its maximum, several different conventions can be used. First, we can indicate the peak offset. This is the offset from the start of the filter output at which the filter output reaches a maximum value. Alternatively, we can use the impulse offset. This is the offset at which the filter would have peaked if the maximum were due to a delta-function like impulse at the input. These quantities are defined in equation ().

Note that in practice, because the chirp signal has to be convolved with the response function $R(f)$ of the detector, the impulse response of the filter is typically a few points longer than the actual chirp signal. For this reason it is smart to assume that the impulse response of your optimal filter is slightly longer (say a hundred points longer) than the actual time-domain length of the corresponding chirp. This safety margin is set with the #define SAFETY statement in the optimal filtering example. You lose a tiny bit of efficiency but reduce the likelihood that boundary effects from the data discontinuity at the start/end of the rectangular window will significantly stimulate the optimal filter output for $i=0,\cdots,N-m-1$. (See Figs.  and to see an illustration of how this windowing discontinuity will corrupt the filter's output.)

We have demonstrated explicitly that with no windowing (or rather, rectangular windowing) of the data, one can find the appropriate correlation between the signal and a filter exactly: the rectangular window has the same effect on the signal as it does on the template (shifting energy into sidelobes in identical fashion). The only complication was that because of the periodic nature of the FFT one has to be caseful about wrap-around errors in relating the output of a filter to the time of occurrence of a signal or impulse.

There is one remaining ugly question. The optimal filter $\tilde Q$ depends upon the noise power spectrum of the detector. In real-world filtering, should this noise power spectrum be calculated with windowed, or non-windowed data? We can determine the correlation between signal and template exactly, with only rectangular windowing, because energy in either of these functions is shifted into sidelobes in identical fashion. However a ``quiet" part of the IFO spectrum can be corrupted by sidelobes of a nearby noisy region. The effect of this is that the signal get rather less weight from this region of frequency space than it ought, in theory, to receive. This would argue for using only properly-windowed data to find the noise power spectrum to use in determining an optimal filter.

In fact, in our experience, it does not make any difference, at least not when you are searching for binary inspiral chirps. The reason is that the SNR obtained in an optimal filter is only sensitive at second order to errors in the optimal filter function. Thus, the errors due to noise sidelobes which appear if you fail to window the data to calculate an optimal filter are typically not large.