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:
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.
zero_pad = np.zeros(data.size - template.size) template_padded = np.append(template, zero_pad) template_fft = np.fft.fft(template_padded)
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 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)
# -- 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)
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.