Find an Inspiral: Optimal matched filter

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.

More to read