Find an Inspiral: Optimal matched filter
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.
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 - datafreq) 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.
More to read
- 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.