8. Calculating the SNR associated with a given astrophysical signal model¶
The example Filtering a TimeSeries to detect gravitational waves showed us we can visually extract a signal from the noise using basic signal-processing techniques.
However, an actual astrophysical search algorithm detects signals by calculating the signal-to-noise ratio (SNR) of data for each in a large bank of signal models, known as templates.
Using PyCBC (the actual search code), we can do that.
First, as always, we fetch some of the public data from the LIGO Open Science Center:
from gwpy.timeseries import TimeSeries data = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)
and condition it by applying a highpass filter at 15 Hz
high = data.highpass(15)
This is important to remove noise at lower frequencies that isn’t accurately calibrated, and swamps smaller noises at higher frequencies.
For this example, we want to calculate the SNR over a 4 second segment, so
we calculate a Power Spectral Density with a 4 second FFT length (using all
of the data), then
crop() the data:
psd = high.psd(4, 2) zoom = high.crop(1126259460, 1126259464)
In order to calculate signal-to-noise ratio, we need a signal model
against which to compare our data.
For this we import
pycbc.waveform.get_fd_waveform and generate a template as a
from pycbc.waveform import get_fd_waveform hp, _ = get_fd_waveform(approximant="IMRPhenomD", mass1=40, mass2=32, f_lower=20, f_final=2048, delta_f=psd.df.value)
At this point we are ready to calculate the SNR, so we import the
pycbc.filter.matched_filter method, and pass it
our template, the data, and the PSD:
from pycbc.filter import matched_filter snr = matched_filter(hp, zoom.to_pycbc(), psd=psd.to_pycbc(), low_frequency_cutoff=15) snrts = TimeSeries.from_pycbc(snr).abs()
Here we have used the
to_pycbc() methods of the
objects to convert from GWpy objects to something that PyCBC functions
can understand, and then used the
to convert back to a GWpy object.
We can plot the SNR
TimeSeries around the region of interest:
plot = snrts.plot() ax = plot.gca() ax.set_xlim(1126259461, 1126259463) ax.set_epoch(1126259462.427) ax.set_ylabel('Signal-to-noise ratio (SNR)') ax.set_title('LIGO-Hanford signal-correlation for GW150914') plot.show()
We can clearly see a large spike (above 17!) at the time of the GW150914 signal! This is, in principle, how the full, blind, CBC search is performed, using all of the available data, and a bank of tens of thousand of signal models.