We will be using the PyCBC library, which is used to study gravitational-wave data, find astrophysical sources due to compact binary mergers, and study their parameters. These are some of the same tools that the LIGO and Virgo collaborations use to find gravitational waves in LIGO/Virgo data
In this tutorial we will walk through how to visualize LIGO/Virgo data and how to perform some basic signal processing on it, including high/low passing, psd estimation, and whitening.
import sys
!{sys.executable} -m pip install pycbc ligo-common --no-cache-dir
Requirement already satisfied: pycbc in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (1.14.4) Requirement already satisfied: lalsuite in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (6.62) Requirement already satisfied: ligo-common in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (1.0.3) Requirement already satisfied: decorator>=3.4.2 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (4.3.0) Requirement already satisfied: Mako>=1.0.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.1.0) Requirement already satisfied: emcee==2.2.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.2.1) Requirement already satisfied: matplotlib>=1.5.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (3.0.0) Requirement already satisfied: lscsoft-glue>=1.59.3 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.0.0) Requirement already satisfied: six>=1.10.0 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.11.0) Requirement already satisfied: requests>=1.2.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.22.0) Requirement already satisfied: ligo-segments in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.2.0) Requirement already satisfied: tqdm in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (4.26.0) Requirement already satisfied: mpld3>=0.3 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (0.3) Requirement already satisfied: astropy>=2.0.3; python_version > "3.0" in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (3.0.4) Requirement already satisfied: cython in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (0.29.14) Requirement already satisfied: scipy>=0.16.0; python_version >= "3.5" in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.1.0) Requirement already satisfied: beautifulsoup4>=4.6.0 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (4.6.3) Requirement already satisfied: pillow in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (5.3.0) Requirement already satisfied: h5py>=2.5 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.8.0) Requirement already satisfied: numpy>=1.13.0; python_version > "3.0" in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (1.16.2) Requirement already satisfied: jinja2 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pycbc) (2.10) Requirement already satisfied: python-dateutil in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from lalsuite) (2.8.1) Requirement already satisfied: MarkupSafe>=0.9.2 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from Mako>=1.0.1->pycbc) (1.1.0) Requirement already satisfied: cycler>=0.10 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from matplotlib>=1.5.1->pycbc) (0.10.0) Requirement already satisfied: kiwisolver>=1.0.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from matplotlib>=1.5.1->pycbc) (1.0.1) Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from matplotlib>=1.5.1->pycbc) (2.3.0) Requirement already satisfied: pyOpenSSL in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from lscsoft-glue>=1.59.3->pycbc) (18.0.0) Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (3.0.4) Requirement already satisfied: certifi>=2017.4.17 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (2018.10.15) Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (1.23) Requirement already satisfied: idna<2.9,>=2.5 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from requests>=1.2.1->pycbc) (2.7) Requirement already satisfied: setuptools in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from kiwisolver>=1.0.1->matplotlib>=1.5.1->pycbc) (41.6.0) Requirement already satisfied: cryptography>=2.2.1 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.3.1) Requirement already satisfied: asn1crypto>=0.21.0 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (0.24.0) Requirement already satisfied: cffi!=1.11.3,>=1.7 in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (1.11.5) Requirement already satisfied: pycparser in /home/nbuser/anaconda3_501/lib/python3.6/site-packages (from cffi!=1.11.3,>=1.7->cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.19)
Below we will view the raw ligo data. You should notice that there is signifant low frequency content (indicated by the large oscilations) and even some DC offset in the data.
%matplotlib inline
# Read in the data around GW150914
from pycbc.catalog import Merger
import pylab
m = Merger('GW150914')
data = {}
for ifo in ['H1', 'L1']:
data[ifo] = m.strain(ifo)
for ifo in data:
pylab.plot(data[ifo].sample_times, data[ifo], label=ifo)
pylab.ylabel('Strain')
pylab.xlabel('GPS Time (s)')
pylab.legend()
pylab.show()
When you just zoom in to one second around GW150914, all you can see is the low frequency behavior of the noise, since it is much louder than the higher frequency noise (and signal).
for ifo in data:
# The time slice method can give just a portion of the time
# series using GPS start and stop times
zoom = data[ifo].time_slice(m.time - 0.5, m.time + 0.5)
pylab.plot(zoom.sample_times, zoom, label=ifo)
pylab.legend()
pylab.show()
Below we apply a highpass filter to the data to suppress the low frequency noise of the instrument. We can see that this brings the dynamic range of the data largely into the same range. However, there is clearly still some dominant frerquencies. To equalize this, we would need to apply a whitening filter.
for ifo in data:
# Apply a highpass filter to the data. This suppresses the low
# frequency content of the data. We choose here a finite-impulse-response (FIR).
# Options
# 1) highpass frequency
# 2) half sample length of highpass filter
#(higher value will give less ripple in passband)
high_data = data[ifo].highpass_fir(15, 512) # Highpass point is 15 Hz
# The time slice method can give just a portion of the time
# series using GPS start and stop times
zoom = high_data.time_slice(m.time - 0.5, m.time + 0.5)
pylab.plot(zoom.sample_times, zoom, label=ifo)
pylab.legend()
pylab.show()
Understanding how the noise power varies over frequency is important for LIGO data analysis. In this section we use a version of Welch's method to estimate the power spectral density of the data.
Note that there is a significant amount of noise at low frequencies (orders of magnitude). Note also that there is a large amount of power at a few specific frequencies. The causes for these include 60 Hz power line noise, violin modes of the hanging mirrors, and various other instrumental resonances. The downturn in power before 2 KHz is due to a low pass filter that was applied prior to resampling the data to 4096 Hz.
for ifo in data:
# This estimates the PSD by sub-dividing the data into overlapping
# 4s long segments. (See Welch's method)
psd = data[ifo].psd(4)
# Note that the psd is a FrequencySeries!
pylab.loglog(psd.sample_frequencies, psd)
pylab.ylabel('$Strain^2 / Hz$')
pylab.xlabel('Frequency (Hz)')
pylab.grid()
pylab.xlim(10, 2048)
pylab.show()
To visualize deviations from the noise, it is useful to "whiten" the data within some frequency range. In this way excesses in the data are visible as deviations from zero. Whitening takes the data and attempts to make the power spectral density flat, so that all frequencies contribute equally.
Below we will whiten the data, and then bandpass the result to focus on a specific frequency range.
# Whiten the data
whitened = {}
for ifo in data:
# This produces a whitened set.
# This works by estimating the power spectral density from the
# data and then flattening the frequency response.
# (1) The first option sets the duration in seconds of each
# sample of the data used as part of the PSD estimate.
# (2) The second option sets the duration of the filter to apply
whitened[ifo] = data[ifo].whiten(4, 4)
zoom = whitened[ifo].time_slice(m.time - 0.5, m.time + 0.5)
pylab.plot(zoom.sample_times, zoom, label=ifo)
pylab.ylabel('Whitened Strain')
pylab.xlabel('Time (s)')
pylab.legend()
pylab.show()
We will now bandpass the data around GW150914 between 30 - 250 Hz. This will remove frequency ranges which won't contribute to this kind of signal and make it possible to see the signal in question.
for ifo in whitened:
# Apply a highpass filter (at 30 Hz) followed by an lowpass filter (at 250 Hz)
bpsd = whitened[ifo].highpass_fir(30, 512).lowpass_fir(250, 512)
zoom = bpsd.time_slice(m.time - 0.5, m.time + 0.5)
pylab.plot(zoom.sample_times, zoom, label=ifo)
pylab.grid()
pylab.legend()
pylab.show()
In the above plot we can see that there is some excess signal that sticks above the noise. Let's zoom around this time now, and align the two time series.
pylab.figure(figsize=[15, 3])
for ifo in whitened:
# Apply a highpass filter (at 30 Hz) followed by an lowpass filter (at 250 Hz)
bpsd = whitened[ifo].highpass_fir(30, 512).lowpass_fir(250, 512)
# We'll choose a tighter zoom here.
zoom = bpsd.time_slice(m.time - 0.2, m.time + .1)
pylab.plot(zoom.sample_times, zoom, label=ifo)
pylab.grid()
pylab.legend()
pylab.show()
pylab.figure(figsize=[15, 3])
for ifo in whitened:
# Apply a highpass filter (at 30 Hz) followed by an lowpass filter (at 250 Hz)
bpsd = whitened[ifo].highpass_fir(30, 512).lowpass_fir(250, 512)
# Now we'll specially align the L1 data. Where does this come from?
# (1) We already knew that the signal was ~ 7 ms separated between detectors.
# (2) The two LIGO interferometers are roughly aligned so that the output of
# one is a sign change of the other for *many* sky locations. This is an
# approximation and doesn't hold for all possible source sky locations.
# A later tutorial will show how to estimate this alignment more precisely.
if ifo == 'L1':
bpsd.roll(int(bpsd.sample_rate * .007))
bpsd *= -1
# We'll choose a tighter zoom here.
zoom = bpsd.time_slice(m.time - 0.2, m.time + .1)
pylab.plot(zoom.sample_times, zoom, label=ifo)
pylab.grid()
pylab.legend()
pylab.show()
We can now see that there is a coherent signal that matches in phase for a few cycles, which ends at about 0.44 on the plot above.
A common way to visualize gravitational-wave data is with a time-frequency representation known as the constant-Q transform. It is similar to a standard spectrogram made with short-time Fourier transforms with the advantage that frequency bins are more sparsely spaced at high freqeuncys.
for ifo in whitened:
# We'll choose a tighter zoom here.
zoom = whitened[ifo].time_slice(m.time - 5, m.time + 5)
# The qtransform method returns a vector of the sample times, frequencies, and a 2-d vector of the
# power in each time-frequency bin. The free parameter is the choice of the Q-value. Larger Q-values
# are generally more appropriate for viewing long duration features of the data and vice versa.
# The options here:
# (1) The time spacing for the output image (i.e. 1 ms in this case)
# (2) The number of frequency bins in the output, logarithmically spaced
# (3) The qrange to maximize over. We'll pick a constant at 8 here
# Typically higher values will be more appropriate for longer duration
# signals
# (4) The frequency range to output
times, freqs, power = zoom.qtransform(.001, logfsteps=100,
qrange=(8, 8),
frange=(20, 512),
)
pylab.figure(figsize=[15, 3])
pylab.pcolormesh(times, freqs, power**0.5)
pylab.xlim(m.time - 0.5, m.time + 0.3)
pylab.title(ifo)
pylab.yscale('log')
pylab.show()
The GW150914 signal is relatively clear in the qtransform output, which is why the Q-transform can be a powerful diagnostic. However, note that quieter signals, especially those with lower masses than GW150914 had, will be harder to spot visually. In practice, we use a technique called matched filtering to find signals in our data.
Use the methods demonstrated above to see if you can visually spot a chirping signal in these data sets. Which of them contain a signal? Which contains just Gaussian noise?
Information that may be useful:
# Download the challenge set files
from six.moves.urllib import request
def get_file(fname):
url = "https://github.com/ahnitz/odw-storage/raw/master/{}"
url = url.format(fname)
request.urlretrieve(url, fname)
print('Getting : {}'.format(url))
files = ['PyCBC_T2_0.gwf', 'PyCBC_T2_1.gwf', 'PyCBC_T2_2.gwf',
'PyCBC_T2_3.gwf', 'PyCBC_T2_4.gwf']
for fname in files:
get_file(fname)
Getting : https://github.com/ahnitz/odw-storage/raw/master/PyCBC_T2_0.gwf Getting : https://github.com/ahnitz/odw-storage/raw/master/PyCBC_T2_1.gwf Getting : https://github.com/ahnitz/odw-storage/raw/master/PyCBC_T2_2.gwf Getting : https://github.com/ahnitz/odw-storage/raw/master/PyCBC_T2_3.gwf Getting : https://github.com/ahnitz/odw-storage/raw/master/PyCBC_T2_4.gwf
import pycbc.frame
# example of reading the strain in (note this is the same as in PyCBC tutorial 1)
ts = pycbc.frame.read_frame('PyCBC_T2_0.gwf', 'H1:TEST-STRAIN', 0, 128)
pylab.plot(ts.sample_times, ts)
pylab.show()