This tutorial assumes that you have already
read in the data and template from the
**Find an Inspiral** tutorial.

This tutorial implements a matched filter, as described in
Allen et. al 2011, arXiv:gr-qc/0509116.

Taking a cross correlation in the time domain is a good technique,
but we know the gravitational wave data has more noise in some
frequency bins than others. To take advantage of this, we can
perform the cross-correlation in the frequency domain, and weight
each frequency bin by the inverse of the noise power. This should
give more weight to frequencies where the noise is low, and so improve
the signal-to-noise ratio of the output. This approach is known
as a "matched filter", and is a well-known technique in signal
processing.
First, take the Fourier Transform (FFT) of the data and the template.

`data_fft=np.fft.fft(data)`

For this to work, we need the template and the data to be the
same length. So, we'll zero-pad the template before we take the FFT:
```
zero_pad = np.zeros(data.size - template.size)
template_padded = np.append(template, zero_pad)
template_fft = np.fft.fft(template_padded)
```

Next, we need an estimate of the noise power in each FFT bin. We'll
take a PSD of the last 4 seconds of data, where we'll assume there's
no signal. The numpy FFT function returns an array with a particular
convention for the order of the frequency bins, so we'll need to
interpolate the PSD to estimate values at each FFT frequency.
```
# -- Calculate the PSD of the data
power_data, freq_psd = plt.psd(data[12*fs:], Fs=fs, NFFT=fs, visible=False)
# -- Interpolate to get the PSD values at the needed frequencies
datafreq = np.fft.fftfreq(data.size)*fs
power_vec = np.interp(datafreq, freq_psd, power_data)
```

Finally, we can multiply the Fourier Space template and data, and divide
by the noise power in each frequency bin. Taking the Inverse
Fourier Transform (IFFT) of the filter output puts it back in the time
domain, so the result will be plotted as a function of time off-set between
the template and the data:
```
# -- Calculate the matched filter output
optimal = data_fft * template_fft.conjugate() / power_vec
optimal_time = 2*np.fft.ifft(optimal)
```

Finally, we can normalize the matched filter output so that we expect a value of 1 at times of just noise. Then, the peak of the matched filter output will tell us the signal-to-noise ratio (SNR) of the signal.

```
# -- Normalize the matched filter output
df = np.abs(datafreq[1] - datafreq[0])
sigmasq = 2*(template_fft * template_fft.conjugate() / power_vec).sum() * df
sigma = np.sqrt(np.abs(sigmasq))
SNR = abs(optimal_time) / (sigma)
# -- Plot the result
plt.figure()
plt.plot(time, SNR)
```

See the normalized matched filter output

Do you see the signal now? When does it occur? Note the matched filter technique expresses time in terms of the off-set between the filter and the data, and so shows the start time of the signal, much like the cross-correlation. This signal was really loud - it has an SNR of around 77. With an optimal matched filter, we can typically detect inspiral signals with an SNR as low as 6.

You can
also download the entire script
described in this tutorial.

- This tutorial is an introduction to digital signal processing. There are many books on digital signal processing available.
- Peter Shawhan's summer school course on gravitational wave data analysis.
- This prescription for a matched filter is described in detail in Allen et. al 2011, arXiv:gr-qc/0509116
- Go back to the Find an Inspiral tutorial.
- Go back to the tutorials index page.