BUG: Fix matched filter SNR calculation in time domain injection (#957)

* TST: add test comparing time domain and frequency domain injections

* BUG: fix calculation of SNR in time domain injection

* FMT: flake fixes for injection test

* BUG: make sure time domain injection metadata is correct

* TST: loosen overlap test

* TST: mark test as flaky

Co-authored-by: Michael J. Williams <michaeljw1@googlemail.com>

* BUG: add missing import

---------

Co-authored-by: Michael J. Williams <michaeljw1@googlemail.com>

Author: ColmTalbot

bilby/gw/detector/__init__.py

@@ -106,13 +106,23 @@ def inject_signal_into_gwpy_timeseries(
         fig.savefig('{}/{}_{}_time_domain_injected_signal'.format(
             outdir, ifo.name, label))
 
-    meta_data = dict()
-    frequency_domain_signal, _ = utils.nfft(signal, waveform_generator.sampling_frequency)
+    meta_data = dict(name=det)
+    frequency_domain_signal, _ = utils.nfft(
+        signal_shifted.value, waveform_generator.sampling_frequency
+    )
+    frequency_domain_data, _ = utils.nfft(
+        signal_and_data.value, waveform_generator.sampling_frequency
+    )
     meta_data['optimal_SNR'] = (
         np.sqrt(ifo.optimal_snr_squared(signal=frequency_domain_signal)).real)
-    meta_data['matched_filter_SNR'] = (
-        ifo.matched_filter_snr(signal=frequency_domain_signal))
-    meta_data['parameters'] = parameters
+    from ..utils import matched_filter_snr
+    meta_data['matched_filter_SNR'] = matched_filter_snr(
+        frequency_domain_signal,
+        frequency_domain_data,
+        ifo.power_spectral_density_array,
+        ifo.duration,
+    )
+    meta_data["parameters"] = parameters
 
     logger.info("Injected signal in {}:".format(ifo.name))
     logger.info("  optimal SNR = {:.2f}".format(meta_data['optimal_SNR']))

test/gw/detector/injection_test.py

@@ -0,0 +1,58 @@
+import bilby
+import pytest
+from gwpy.frequencyseries import FrequencySeries
+
+
+@pytest.mark.flaky(reruns=3)
+def test_injection_into_timeseries_matches_ifo_injections():
+    """
+    Test that injecting into a gwpy timeseries agree with the frequency-domain
+    injection.
+
+    Caveats:
+    - there is a discrete alignment of the maximum strain in the time domain
+      injection that means the overlap and SNRs don't exactly match. This is
+      especially pronounced when the sampling frequency is low.
+    """
+    duration = 8
+    sampling_frequency = 16384
+    ifo = bilby.gw.detector.get_empty_interferometer("H1")
+    ifo.set_strain_data_from_zero_noise(
+        sampling_frequency=sampling_frequency, duration=duration
+    )
+    wfg = bilby.gw.waveform_generator.WaveformGenerator(
+        frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+        duration=duration,
+        sampling_frequency=sampling_frequency,
+    )
+    priors = bilby.gw.prior.BBHPriorDict()
+    priors["geocent_time"] = 6.0
+    tseries = ifo.strain_data.to_gwpy_timeseries()
+    parameters = priors.sample()
+    ifo.inject_signal(parameters, waveform_generator=wfg, raise_error=False)
+    signal_and_data, meta_data = bilby.gw.detector.inject_signal_into_gwpy_timeseries(
+        tseries,
+        waveform_generator=wfg,
+        parameters=parameters,
+        det="H1",
+    )
+
+    assert meta_data.keys() == ifo.meta_data.keys()
+    assert meta_data["parameters"] == ifo.meta_data["parameters"]
+    for key in ["optimal_SNR", "matched_filter_SNR"]:
+        assert abs(meta_data[key] - ifo.meta_data[key]) < 1e-5
+
+    asd = FrequencySeries(
+        ifo.amplitude_spectral_density_array,
+        frequencies=ifo.frequency_array,
+        unit="1/Hz",
+    )
+    whitened_data_1 = signal_and_data.whiten(asd=asd).value
+    whitened_data_2 = ifo.whitened_time_domain_strain
+
+    mismatch = 1 - (
+        sum(whitened_data_1 * whitened_data_2)
+        / sum(whitened_data_1**2)**0.5
+        / sum(whitened_data_2**2)**0.5
+    )
+    assert mismatch < 3e-3