.. _gwpy-example-timeseries-inject:

.. sectionauthor:: Alex Urban <alexander.urban@ligo.org>

.. currentmodule:: gwpy.timeseries

Inject a known signal into a `TimeSeries`
#########################################

It can often be useful to add some known signal to an inherently random
or noisy timeseries. For example, one might want to examine what
would happen if a binary black hole merger signal occured at or near
the time of a glitch. In LIGO data analysis, this procedure is referred
to as an _injection_.

In the example below, we will create a stream of random, white Gaussian
noise, then inject a simulation of GW150914 into it at a known time.

First, we prepare one second of Gaussian noise:


.. plot::
   :context: reset
   :nofigs:
   :include-source:

   from numpy import random
   from gwpy.timeseries import TimeSeries
   noise = TimeSeries(random.normal(scale=.1, size=16384), sample_rate=16384)

Then we can download a simulation of the GW150914 signal from GWOSC:


.. plot::
   :context:
   :nofigs:
   :include-source:

   from astropy.utils.data import get_readable_fileobj
   url = ("https://www.gw-openscience.org/s/events/GW150914/P150914/"
          "fig2-unfiltered-waveform-H.txt")
   with get_readable_fileobj(url) as f:
       signal = TimeSeries.read(f, format='txt')
   signal.t0 = .5  # make sure this intersects with noise time samples

Note, since this simulation cuts off before a certain time, it is
important to taper its ends to zero to avoid ringing artifacts.
We can accomplish this using the
:meth:`~gwpy.timeseries.TimeSeries.taper` method.


.. plot::
   :context:
   :nofigs:
   :include-source:

   signal = signal.taper()

Since the time samples overlap, we can inject this into our noise data
using :meth:`~gwpy.types.series.Series.inject`:


.. plot::
   :context:
   :nofigs:
   :include-source:

   data = noise.inject(signal)

Finally, we can visualize the full process in the time domain:


.. plot::
   :context:
   :include-source:

   from gwpy.plot import Plot
   plot = Plot(noise, signal, data, separate=True, sharex=True, sharey=True)
   plot.gca().set_epoch(0)
   plot.show()

We can clearly see that the loud GW150914-like signal has been layered
on top of Gaussian noise with the correct amplitude and phase evolution.