.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/signal/gw150914.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_signal_gw150914.py: .. sectionauthor:: Duncan Macleod .. currentmodule:: gwpy.timeseries Filtering a `TimeSeries` to detect gravitational waves ###################################################### The raw 'strain' output of the LIGO detectors is recorded as a `TimeSeries` with contributions from a large number of known and unknown noise sources, as well as possible gravitational wave signals. In order to uncover a real signal we need to filter out noises that otherwise hide the signal in the data. We can do this by using the :mod:`gwpy.signal` module to design a digital filter to cut out low and high frequency noise, as well as notch out fixed frequencies polluted by known artefacts. .. GENERATED FROM PYTHON SOURCE LINES 37-40 Data access ----------- First we download the raw strain data from the GWOSC public archive: .. GENERATED FROM PYTHON SOURCE LINES 40-44 .. code-block:: Python from gwpy.timeseries import TimeSeries hdata = TimeSeries.fetch_open_data("H1", 1126259446, 1126259478) .. GENERATED FROM PYTHON SOURCE LINES 45-52 Data conditioning ----------------- We can design a zero-pole-gain filter to remove the extranious noise. First we import the `gwpy.signal.filter_design` module and create a :meth:`~gwpy.signal.filter_design.bandpass` filter to remove both low and high frequency content .. GENERATED FROM PYTHON SOURCE LINES 52-56 .. code-block:: Python from gwpy.signal import filter_design bp = filter_design.bandpass(50, 250, hdata.sample_rate) .. GENERATED FROM PYTHON SOURCE LINES 57-60 Now we want to combine the bandpass with a series of :meth:`~gwpy.signal.filter_design.notch` filters, so we create those for the first three harmonics of the 60 Hz AC mains power: .. GENERATED FROM PYTHON SOURCE LINES 60-68 .. code-block:: Python notch_frequencies = [ 60, 120, 180, ] notches = [filter_design.notch(f, hdata.sample_rate) for f in notch_frequencies] .. GENERATED FROM PYTHON SOURCE LINES 69-70 and concatenate each of our filters together to create a single ZPK: .. GENERATED FROM PYTHON SOURCE LINES 70-73 .. code-block:: Python zpk = filter_design.concatenate_zpks(bp, *notches) .. GENERATED FROM PYTHON SOURCE LINES 74-77 Now, we can apply our combined filter to the data, using `filtfilt=True` to filter both backwards and forwards to preserve the correct phase at all frequencies .. GENERATED FROM PYTHON SOURCE LINES 77-80 .. code-block:: Python hfilt = hdata.filter(zpk, filtfilt=True) .. GENERATED FROM PYTHON SOURCE LINES 81-93 .. note:: The :mod:`~gwpy.signal.filter_design` methods return digital filters by default, so we apply them using `TimeSeries.filter`. If we had analogue filters (perhaps by passing `analog=True` to the filter design method), the easiest application would be `TimeSeries.zpk` The :mod:`~gwpy.signal.filter_design` methods return infinite impulse response filters by default, which, when applied, corrupt a small amount of data at the beginning and the end of our original `TimeSeries`. We can discard those data using the :meth:`~TimeSeries.crop` method (for consistency we apply this to both data series): .. GENERATED FROM PYTHON SOURCE LINES 93-97 .. code-block:: Python hdata = hdata.crop(*hdata.span.contract(1)) hfilt = hfilt.crop(*hfilt.span.contract(1)) .. GENERATED FROM PYTHON SOURCE LINES 98-102 Visualisation ------------- Finally, we can :meth:`~TimeSeries.plot` the original and filtered data, adding some code to prettify the figure: .. GENERATED FROM PYTHON SOURCE LINES 102-132 .. code-block:: Python from gwpy.plot import Plot plot = Plot( hdata, hfilt, figsize=[12, 6], separate=True, sharex=True, color="gwpy:ligo-hanford", ) ax1, ax2 = plot.axes ax1.set_title("LIGO-Hanford strain data around GW150914") ax1.text( 1.0, 1.01, "Unfiltered data", transform=ax1.transAxes, ha="right", ) ax1.set_ylabel("Amplitude [strain]", y=-0.2) ax2.set_ylabel("") ax2.text( 1.0, 1.01, r"50-250\,Hz bandpass, notches at 60, 120, 180 Hz", transform=ax2.transAxes, ha="right", ) plot.show() .. image-sg:: /examples/signal/images/sphx_glr_gw150914_001.png :alt: LIGO-Hanford strain data around GW150914 :srcset: /examples/signal/images/sphx_glr_gw150914_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 133-135 We see now a spike around 16 seconds into the data, so let's zoom into that time (and prettify): .. GENERATED FROM PYTHON SOURCE LINES 135-144 .. code-block:: Python plot = hfilt.plot(color="gwpy:ligo-hanford") ax = plot.gca() ax.set_title("LIGO-Hanford strain data around GW150914") ax.set_ylabel("Amplitude [strain]") ax.set_xlim(1126259462, 1126259462.6) ax.set_xscale("seconds", epoch=1126259462) plot.show() .. image-sg:: /examples/signal/images/sphx_glr_gw150914_002.png :alt: LIGO-Hanford strain data around GW150914 :srcset: /examples/signal/images/sphx_glr_gw150914_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 145-154 Congratulations, you have succesfully filtered LIGO data to uncover the first ever directly-detected gravitational wave signal, GW150914! Adding the second LIGO detector ------------------------------- But wait, what about LIGO-Livingston? We can easily add that to our figure by following the same procedure. First, we load the new data .. GENERATED FROM PYTHON SOURCE LINES 154-158 .. code-block:: Python ldata = TimeSeries.fetch_open_data("L1", 1126259446, 1126259478) lfilt = ldata.filter(zpk, filtfilt=True) .. GENERATED FROM PYTHON SOURCE LINES 159-163 The article announcing the detector told us that the signals were separated by 6.9 ms between detectors, and that the relative orientation of those detectors means that we need to invert the data from one before comparing them, so we apply those corrections: .. GENERATED FROM PYTHON SOURCE LINES 163-167 .. code-block:: Python lfilt.shift("6.9ms") lfilt *= -1 .. GENERATED FROM PYTHON SOURCE LINES 168-169 and finally make a new plot with both detectors: .. GENERATED FROM PYTHON SOURCE LINES 169-182 .. code-block:: Python plot = Plot(figsize=[12, 4]) ax = plot.gca() ax.plot(hfilt, label="LIGO-Hanford", color="gwpy:ligo-hanford") ax.plot(lfilt, label="LIGO-Livingston", color="gwpy:ligo-livingston") ax.set_title("LIGO strain data around GW150914") ax.set_xlim(1126259462, 1126259462.6) ax.set_xscale("seconds", epoch=1126259462) ax.set_ylabel("Amplitude [strain]") ax.set_ylim(-1e-21, 1e-21) ax.legend() plot.show() .. image-sg:: /examples/signal/images/sphx_glr_gw150914_003.png :alt: LIGO strain data around GW150914 :srcset: /examples/signal/images/sphx_glr_gw150914_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 183-187 The above filtering is (almost) the same as what was applied to LIGO data to produce Figure 1 in `Abbott et al. (2016) `_ (the joint LSC-Virgo publication announcing this detection). .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.187 seconds) .. _sphx_glr_download_examples_signal_gw150914.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: gw150914.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: gw150914.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: gw150914.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_