Detection of Gravitational Waves

Detection of Gravitational Waves#

Overview#

Gravitational waves (GWs) were predicted in 1916 by Einstein’s general theory of relativity, in which mass and energy curve spacetime. By analogy with electromagnetism, a time-varying mass quadrupole moment of a source produces transverse waves of spatial strain that travel at the speed of light.

The first evidence for gravitational waves came from the binary pulsar PSR B1913+16, discovered by Hulse and Taylor in 1974 [1]: a pulsar in orbit around another neutron star. Over the following years its orbit was observed to shrink at exactly the rate predicted by general relativity for energy lost to gravitational radiation [2]. This indirect evidence earned the 1993 Nobel Prize in Physics.

The first direct detection came on September 14, 2015, when the two Advanced LIGO (Laser Interferometer Gravitational-Wave Observatory) detectors recorded the signal of two black holes merging into one. The LIGO Scientific Collaboration and the Virgo Collaboration announced it in February 2016 [3]. This observation, and the many detections since, provide new tests of general relativity in its strong-field regime, and GW observations have become an important new way to learn about the Universe.

In this project you will reproduce key results of Ref. [3] for the event GW150914 with LIGO open data, and then develop a machine-learning method to detect GW signals in simulated data.

Data Sources#

  • LIGO open data from the Gravitational Wave Open Science Center (GWOSC) [4]: https://gwosc.org/data. GWOSC provides public access to the data of LIGO, Virgo, GEO600 and KAGRA. We download the data around GW150914 directly in the notebook with the gwpy package [5].

  • Simulated GW challenge data from the G2Net Gravitational Wave Detection Kaggle competition (2021), used in Question 06. The dataset is large (about 77 GB), and using it requires a Kaggle account and accepting the competition rules. See Question 06 for how to work with it.

Questions#

Answer each question in a markdown cell directly below it (add code cells where needed). For the background questions, aim for one to two well-organized paragraphs, and cite your sources.

Question 01#

What are gravitational waves? How are they generated, and what kinds of astrophysical sources produce waves that we can detect? How fast do they propagate, and how has this been tested experimentally? What could gravitational waves tell us about the early Universe?

Question 02#

How are gravitational waves detected? Describe the principle of a laser interferometer such as LIGO, and the size of the strain it must measure. Where are the two LIGO detectors, and which other detectors form today’s worldwide network? Why is it important to have several detectors far apart?

Filtering a TimeSeries to detect gravitational waves#

Importing and Setup#

We use the gwpy package [5] to download the data and for signal processing. Install it with the cell below. If you are working in Google Colab and an import fails afterwards, restart the runtime and run the notebook again from the top.

Note: do not add %matplotlib inline to this notebook. Inline plotting is already the default in Jupyter and Colab, and with recent versions of matplotlib that line breaks gwpy’s plotting functions.

%pip install --quiet gwpy

import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import numpy as np
from scipy.interpolate import interp1d
import warnings
warnings.filterwarnings('ignore')

The raw “strain” output of the LIGO detectors is recorded as a TimeSeries. It contains contributions from a large number of known and unknown noise sources, as well as possible gravitational-wave signals.

To uncover a signal, we need to filter out the noise that otherwise hides it. We do this with the gwpy.signal module. We design a digital filter that cuts out low- and high-frequency noise, and “notches” out fixed frequencies polluted by known artefacts.

First we download 32 seconds of raw LIGO-Hanford (H1) strain data around the time of GW150914 from the GWOSC public archive. The times are GPS times, in seconds; GW150914 occurred at GPS time 1126259462.4.

from gwpy.timeseries import TimeSeries
hdata = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)

Next we design a zero-pole-gain (ZPK) filter to remove the extraneous noise. First we import the gwpy.signal.filter_design module and create a bandpass() filter to remove both low- and high-frequency content:

from gwpy.signal import filter_design
bp = filter_design.bandpass(50, 250, hdata.sample_rate)

Now we combine the bandpass filter with a series of notch() filters for the first three harmonics of the 60 Hz US mains power:

notches = [filter_design.notch(line, hdata.sample_rate) for line in (60, 120, 180)]

and concatenate the filters into a single ZPK model:

zpk = filter_design.concatenate_zpks(bp, *notches)

Now we apply the combined filter to the data. We use filtfilt=True, which filters both backwards and forwards and so preserves the correct phase at all frequencies:

hfilt = hdata.filter(zpk, filtfilt=True)

Note: the filter_design methods return digital filters by default, so we apply them using TimeSeries.filter. If we had analogue filters (e.g. by passing analog=True to the filter-design method), the easiest way to apply them would be TimeSeries.zpk.

The filter_design methods return infinite-impulse-response filters by default. When applied, these corrupt a small amount of data at the beginning and end of the TimeSeries, so we discard one second at each end with the crop() method (for consistency, we apply this to both data series):

hdata = hdata.crop(*hdata.span.contract(1))
hfilt = hfilt.crop(*hfilt.span.contract(1))

Finally, we plot() the original and filtered data, with some code to prettify the figure:

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()

We now see a spike about 16 seconds into the data, so let’s zoom in on that time:

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()

Question 03#

Congratulations, you have successfully filtered LIGO data to uncover the first directly detected gravitational-wave signal, GW150914!

But what about LIGO-Livingston (L1)? Add the data from LIGO-Livingston to the figure by following the same procedure. You can load the data with:

ldata = TimeSeries.fetch_open_data('L1', 1126259446, 1126259478)

Hint: Ref. [3] reports that the signal arrived at the two detectors at slightly different times, and that because of their relative orientations, one data stream must be inverted before comparing them. If you do it right, you should be able to make a comparison plot resembling the upper-right panel of Fig. 1 of Ref. [3]. (That figure uses a 35–350 Hz bandpass filter instead of 50–250 Hz, so the curves will differ slightly.)

What does the time difference tell you about the direction to the source?

Data in the Frequency Domain#

Plotting the data in the frequency domain shows their frequency content. One way to visualize it is the amplitude spectral density (ASD), the square root of the power spectral density (PSD). The PSD is an average of the squared magnitudes of the fast Fourier transforms (FFTs) of short stretches of the data.

The ASD is an estimate of the “strain-equivalent noise” of the detectors as a function of frequency, which limits their ability to detect GW signals. It has units of strain/\(\sqrt{\text{Hz}}\). To find the root-mean-square (rms) strain noise in a frequency band, integrate the square of the ASD over that band and take the square root.

You will see that the unfiltered LIGO data are dominated by low-frequency noise. Without the signal processing above, there would be no way to see the signal.

fs = int(hdata.sample_rate.value)   # sampling rate: 4096 Hz
NFFT = 1 * fs                       # number of samples per FFT (1 second)
fmin, fmax = 10, 2000

Pxx_H1data, freqs = mlab.psd(hdata.value, Fs=fs, NFFT=NFFT)
Pxx_H1filt, freqs = mlab.psd(hfilt.value, Fs=fs, NFFT=NFFT)

# Interpolations of the PSDs (useful e.g. for whitening)
psd_H1data = interp1d(freqs, Pxx_H1data)
psd_H1filt = interp1d(freqs, Pxx_H1filt)

plt.figure()
plt.loglog(freqs, np.sqrt(Pxx_H1data), 'r', label='H1 strain')
plt.loglog(freqs, np.sqrt(Pxx_H1filt), 'g', label='H1 strain (filtered)')
plt.axis([fmin, fmax, 1e-24, 1e-19])
plt.grid(True)
plt.ylabel(r'ASD (strain/$\sqrt{\mathrm{Hz}}$)')
plt.xlabel('Frequency (Hz)')
plt.legend(loc='upper center')
plt.title('Advanced LIGO strain data near GW150914')
plt.show()

We only plot the data between fmin = 10 Hz and fmax = 2000 Hz.

  • Below fmin, the data are not properly calibrated. That is fine, because the noise there is so high that LIGO cannot detect gravitational waves from astrophysical sources in that band.

  • The sampling rate is fs = 4096 Hz (\(2^{12}\) Hz), so the data cannot capture frequency content above the Nyquist frequency, fs/2 = 2048 Hz. That is also fine, because GW150914 only has detectable frequency content between about 20 and 300 Hz.

You can see strong spectral lines in the data, all of instrumental origin. Some are engineered into the detectors, such as the mirror-suspension resonances near 500 Hz and their harmonics, calibration lines, and control dither lines. Others, such as 60 Hz and its harmonics, are unwanted.

Question 04#

Make the analogous ASD plot for the LIGO-Livingston (L1) data, and a plot that directly compares the H1 and L1 noise spectra, analogous to the strain-noise spectra shown in Fig. 3 of Ref. [3]. In which frequency range are the detectors most sensitive? How does that compare with the frequency content of GW150914?

Evolution of the Frequency over Time: The Spectrogram#

One of the most useful ways to visualize gravitational-wave data is a spectrogram, which shows the frequency content of the data in successive short time intervals. The cell below calculates a spectrogram of the filtered H1 data with TimeSeries.spectrogram2(), using FFTs of 1/16 s with an overlap of 15/256 s between consecutive FFTs. spectrogram2() returns a power spectral density (PSD) spectrogram, so we take the square root (** (1/2.)) to convert it into an amplitude spectral density.

specgram = hfilt.spectrogram2(fftlength=1/16., overlap=15/256.) ** (1/2.)
specgram = specgram.crop_frequencies(20)  # drop everything below the high-pass frequency

plot = specgram.plot(norm='log', cmap='viridis', yscale='log')
ax = plot.gca()
ax.set_title('LIGO-Hanford strain data around GW150914')
ax.set_xlim(1126259462, 1126259463)
ax.colorbar(label=r'Strain ASD [1/$\sqrt{\mathrm{Hz}}$]')
plot.show()

There is a hint of the famous “chirp” signal (around \(t = 0.4\) s in the plot) from the inspiral and merger of the two black holes. The frequency rises as the black holes spiral together. With the additional processing of the next question, it becomes much more visible.

Q-transform#

One of the most useful tools for filtering and visualizing short-duration features in a TimeSeries is the Q-transform. The LIGO Scientific Collaboration and the Virgo Collaboration use it regularly to make high-resolution time–frequency maps of transient noise (“glitches”) and of potential gravitational-wave signals. It was also used to visualize GW150914 in the bottom row of Fig. 1 of Ref. [3].

Question 05#

Following the gwpy Q-transform example, make Q-transform spectrograms of the LIGO-Hanford and LIGO-Livingston data, using TimeSeries.q_transform on the unfiltered data (the Q-transform whitens the data itself).

  • Compute the Q-transform for a short interval around GW150914 (GPS time 1126259462.4), e.g. outseg=(t0 - 0.2, t0 + 0.1). Show the time axis from 0.17 s before to 0.03 s after \(t_0\), and the frequency axis from 20 to 500 Hz on a logarithmic scale.

  • Compare your plots with the bottom row of Fig. 1 of Ref. [3].

  • Explain why the frequency of the signal increases with time (“chirp”), and why the signal ends abruptly.

Machine Learning Approaches to Gravitational Wave Detection#

GW signals are buried in detector noise, and machine-learning methods can be powerful for finding them. You will now develop an ML method to detect GW signals from binary black hole mergers, using simulated data from the G2Net Gravitational Wave Detection challenge.

In this dataset, each sample contains 2 seconds of simulated data, sampled at 2048 Hz, from three detectors: LIGO Hanford, LIGO Livingston, and Virgo. Check the exact format yourself. Each sample has a label: 1 if a simulated signal is present, 0 if it is noise only. Submissions to the challenge were ranked by the area under the ROC curve (AUC). The best teams reached an AUC of about 0.885 on the competition’s test data.

Working with the data. The full dataset (about 77 GB) is too large to download for this project. We recommend one of the following:

  • Work in a Kaggle Notebook attached to the competition (Code → New Notebook on the competition page). The data are then available without downloading, and Kaggle provides a free GPU quota. You need a Kaggle account, and you must accept the competition rules.

  • Download a subset of the training files (e.g. 20,000–50,000 samples) with the Kaggle API.

Question 06#

Develop, train and evaluate an ML model that detects GW signals from binary black hole mergers in the G2Net data.

  • Preprocessing: the signals are far below the noise level. Explain and apply standard GW preprocessing: a window function (e.g. Tukey), a bandpass filter (e.g. 20–500 Hz), and whitening, i.e. dividing the data in the frequency domain by an estimate of the noise ASD. A good noise ASD can be estimated by averaging the spectra of many noise-only (label 0) samples. Plot a sample before and after preprocessing.

  • Model: a 1D CNN acting on the three whitened time series is a good start. Alternatively, convert each sample into a time–frequency image (e.g. a constant-Q transform, as in Question 05) and use a 2D CNN.

  • Evaluation: keep a validation set that is not used for training. Report the ROC curve and AUC, and compare with the top challenge scores. Show the distribution of the network output for signal and noise samples.

It is more important to get one ML model working well than to “play” with many models superficially. You may take inspiration and code from publicly available challenge solutions (see the competition’s Code tab) or from published work, e.g. Ref. [6], which compares a Conv1D network with a generative pre-trained Transformer. You must cite all your sources.

Optional extension: apply your trained model to the real LIGO data around GW150914 from the first part of this project. You will need to resample to 2048 Hz, take 2-second windows, and decide what to do about Virgo, which was not observing at the time. Does your model detect the event? What does this tell you about the difference between simulated and real detector noise?

References#

[1] R. A. Hulse, J. H. Taylor, “Discovery of a pulsar in a binary system,” Astrophys. J. 195, L51 (1975). doi:10.1086/181708, ADS

[2] J. H. Taylor, J. M. Weisberg, “A new test of general relativity: gravitational radiation and the binary pulsar PSR 1913+16,” Astrophys. J. 253, 908 (1982). doi:10.1086/159690

[3] B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116, 061102 (2016). doi:10.1103/PhysRevLett.116.061102, arXiv:1602.03837

[4] R. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), “Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo,” SoftwareX 13, 100658 (2021). doi:10.1016/j.softx.2021.100658

[5] D. M. Macleod et al., “GWpy: A Python package for gravitational-wave astrophysics,” SoftwareX 13, 100657 (2021). doi:10.1016/j.softx.2021.100657. Documentation: https://gwpy.readthedocs.io

[6] A. Dana Ruiz, Detection of Gravitational Wave signals using Machine Learning methods and Generative Pre-trained Transformers, thesis, Universitat de Barcelona (2023). https://diposit.ub.edu/handle/2445/201012

[7] G2Net Gravitational Wave Detection, Kaggle competition (2021). https://www.kaggle.com/competitions/g2net-gravitational-wave-detection

Acknowledgements#

  • Initial version: Mark Neubauer

  • Revised: 2026

  • The filtering, spectrogram and Q-transform examples are adapted from the gwpy documentation and the GWOSC tutorials.

© Copyright 2026