This tutorial is intended for educational purposes. The code shown here is not used to produce results papers published by the LIGO Scientific Collaboration, which instead rely on special purpose analysis software packages.¶

For publicly available, gravitational-wave software analysis packages that are used to produce LSC and Virgo Collaboration results papers, see https://losc.ligo.org/software/.¶

This tutorial also assumes that you know a bit about signal processing of digital time series data (or want to learn!). This includes power spectral densities, spectrograms, digital filtering, whitening, audio manipulation. This is a vast and complex set of topics, but we will cover many of the basics in this tutorial.

# Standard python numerical analysis imports:importnumpyasnpfromscipyimportsignalfromscipy.interpolateimportinterp1dfromscipy.signalimportbutter,filtfilt,iirdesign,zpk2tf,freqzimporth5pyimportjson# the IPython magic below must be commented out in the .py file, since it doesn't work there.%matplotlib inline
%config InlineBackend.figure_format = 'retina'
importmatplotlib.pyplotaspltimportmatplotlib.mlabasmlab# LIGO-specific readligo.py importreadligoasrl# you might get a matplotlib warning here; you can ignore it.

Read the event properties from a local json file (download in advance):¶

In [49]:

# Read the event properties from a local json filefnjson="BBH_events_v3.json"try:events=json.load(open(fnjson,"r"))exceptIOError:print("Cannot find resource file "+fnjson)print("You can download it from https://losc.ligo.org/s/events/"+fnjson)print("Quitting.")quit()# did the user select the eventname ?try:events[eventname]except:print('You must select an eventname that is in '+fnjson+'! Quitting.')quit()

In [50]:

# Extract the parameters for the desired event:event=events[eventname]fn_H1=event['fn_H1']# File name for H1 datafn_L1=event['fn_L1']# File name for L1 datafn_template=event['fn_template']# File name for template waveformfs=event['fs']# Set sampling ratetevent=event['tevent']# Set approximate event GPS timefband=event['fband']# frequency band for bandpassing signalprint("Reading in parameters for event "+event["name"])print(event)

#----------------------------------------------------------------# Load LIGO data from a single file.# FIRST, define the filenames fn_H1 and fn_L1, above.#----------------------------------------------------------------try:# read in data from H1 and L1, if available:strain_H1,time_H1,chan_dict_H1=rl.loaddata(fn_H1,'H1')strain_L1,time_L1,chan_dict_L1=rl.loaddata(fn_L1,'L1')except:print("Cannot find data files!")print("You can download them from https://losc.ligo.org/s/events/"+eventname)print("Quitting.")quit()

NOTE that in general, LIGO strain time series data has gaps (filled with NaNs) when the detectors are not taking valid ("science quality") data. Analyzing these data requires the user to
loop over "segments" of valid data stretches.

In this tutorial, for simplicity, we assume there are no data gaps - this will not work for all times! See the
notes on segments for details.

# plot +- deltat seconds around the event:# index into the strain time series for this time interval:deltat=5indxt=np.where((time>=tevent-deltat)&(time<tevent+deltat))print(tevent)ifmake_plots:plt.figure()plt.plot(time[indxt]-tevent,strain_H1[indxt],'r',label='H1 strain')plt.plot(time[indxt]-tevent,strain_L1[indxt],'g',label='L1 strain')plt.xlabel('time (s) since '+str(tevent))plt.ylabel('strain')plt.legend(loc='lower right')plt.title('Advanced LIGO strain data near '+eventname)plt.savefig(eventname+'_strain.'+plottype)

1167559936.6

The data are dominated by low frequency noise; there is no way to see a signal here, without some signal processing.

Plotting these data in the Fourier domain gives us an idea of the frequency content of the data. A way to visualize the frequency content of the data is to plot the amplitude spectral density, ASD.

The ASDs are the square root of the power spectral densities (PSDs), which are averages of the square of the fast fourier transforms (FFTs) of the data.

They are an estimate of the "strain-equivalent noise" of the detectors versus frequency,
which limit the ability of the detectors to identify GW signals.

They are in units of strain/rt(Hz).
So, if you want to know the root-mean-square (rms) strain noise in a frequency band,
integrate (sum) the squares of the ASD over that band, then take the square-root.

There's a signal in these data!
For the moment, let's ignore that, and assume it's all noise.

In [54]:

make_psds=1ifmake_psds:# number of sample for the fast fourier transform:NFFT=4*fsPxx_H1,freqs=mlab.psd(strain_H1,Fs=fs,NFFT=NFFT)Pxx_L1,freqs=mlab.psd(strain_L1,Fs=fs,NFFT=NFFT)# We will use interpolations of the ASDs computed above for whitening:psd_H1=interp1d(freqs,Pxx_H1)psd_L1=interp1d(freqs,Pxx_L1)# Here is an approximate, smoothed PSD for H1 during O1, with no lines. We'll use it later. Pxx=(1.e-22*(18./(0.1+freqs))**2)**2+0.7e-23**2+((freqs/2000.)*4.e-23)**2psd_smooth=interp1d(freqs,Pxx)ifmake_plots:# plot the ASDs, with the template overlaid:f_min=20.f_max=2000.plt.figure(figsize=(10,8))plt.loglog(freqs,np.sqrt(Pxx_L1),'g',label='L1 strain')plt.loglog(freqs,np.sqrt(Pxx_H1),'r',label='H1 strain')plt.loglog(freqs,np.sqrt(Pxx),'k',label='H1 strain, O1 smooth model')plt.axis([f_min,f_max,1e-24,1e-19])plt.grid('on')plt.ylabel('ASD (strain/rtHz)')plt.xlabel('Freq (Hz)')plt.legend(loc='upper center')plt.title('Advanced LIGO strain data near '+eventname)plt.savefig(eventname+'_ASDs.'+plottype)

NOTE that we only plot the data between f_min = 20 Hz and f_max = 2000 Hz.

Below f_min, the data are not properly calibrated. That's OK, because the noise is so high below f_min that LIGO cannot sense gravitational wave strain from astrophysical sources in that band.

The sample rate is fs = 4096 Hz (2^12 Hz), so the data cannot capture frequency content above the Nyquist frequency = fs/2 = 2048 Hz. That's OK, because our events only have detectable frequency content in the range given by fband, defined above; the upper end will (almost) always be below the Nyquist frequency. We set f_max = 2000, a bit below Nyquist.

You can see strong spectral lines in the data; they are all of instrumental origin. Some are engineered into the detectors (mirror suspension resonances at ~500 Hz and harmonics, calibration lines, control dither lines, etc) and some (60 Hz and harmonics) are unwanted. We'll return to these, later.

You can't see the signal in this plot, since it is relatively weak and less than a second long, while this plot averages over 32 seconds of data. So this plot is entirely dominated by instrumental noise.

The smooth model is hard-coded and tuned by eye; it won't be right for arbitrary times. We will only use it below for things that don't require much accuracy.

A standard metric that LIGO uses to evaluate the sensitivity of our detectors, based on the detector noise ASD, is the BNS range.

This is defined as the distance to which a LIGO detector can register a BNS signal with a single detector signal-to-noise ratio (SNR) of 8, averaged over source direction and orientation. Here, SNR 8 is used as a nominal detection threshold, similar to typical CBC detection thresholds of SNR 6-8.

We take each neutron star in the BNS system to have a mass of 1.4 times the mass of the sun, and negligible spin.

GWs from BNS mergers are like "standard sirens"; we know their amplitude at the source from theoretical calculations. The amplitude falls off like 1/r, so their amplitude at the detectors on Earth tells us how far away they are. This is great, because it is hard, in general, to know the distance to astronomical sources.

The amplitude at the source is computed in the post-Newtonian "quadrupole approximation". This is valid for the inspiral phase only, and is approximate at best; there is no simple expression for the post-inspiral (merger and ringdown) phase. So this won't work for high-mass binary black holes like GW150914, which have a lot of signal strength in the post-inspiral phase.

But, in order to use them as standard sirens, we need to know the source direction and orientation relative to the detector and its "quadrupole antenna pattern" response to such signals. It is a standard (if non-trivial) computation to average over all source directions and orientations; the average amplitude is 1./2.2648 times the maximum value.

NOTE that, since mass is the source of gravity and thus also of gravitational waves, systems with higher masses (such as the binary black hole merger GW150914) are much "louder" and can be detected to much higher distances than the BNS range. We'll compute the BBH range, using a template with specific masses, below.

From the ASD above, we can see that the data are very strongly "colored" - noise fluctuations are much larger at low and high frequencies and near spectral lines, reaching a roughly flat ("white") minimum in the band around 80 to 300 Hz.

We can "whiten" the data (dividing it by the noise amplitude spectrum, in the fourier domain), suppressing the extra noise at low frequencies and at the spectral lines, to better see the weak signals in the most sensitive band.

Whitening is always one of the first steps in astrophysical data analysis (searches, parameter estimation).
Whitening requires no prior knowledge of spectral lines, etc; only the data are needed.

To get rid of remaining high frequency noise, we will also bandpass the data.

The resulting time series is no longer in units of strain; now in units of "sigmas" away from the mean.

We will plot the whitened strain data, along with the signal template, after the matched filtering section, below.

In [56]:

# function to whiten datadefwhiten(strain,interp_psd,dt):Nt=len(strain)freqs=np.fft.rfftfreq(Nt,dt)freqs1=np.linspace(0,2048.,Nt/2+1)# whitening: transform to freq domain, divide by asd, then transform back, # taking care to get normalization right.hf=np.fft.rfft(strain)norm=1./np.sqrt(1./(dt*2))white_hf=hf/np.sqrt(interp_psd(freqs))*normwhite_ht=np.fft.irfft(white_hf,n=Nt)returnwhite_htwhiten_data=1ifwhiten_data:# now whiten the data from H1 and L1, and the template (use H1 PSD):strain_H1_whiten=whiten(strain_H1,psd_H1,dt)strain_L1_whiten=whiten(strain_L1,psd_L1,dt)# We need to suppress the high frequency noise (no signal!) with some bandpassing:bb,ab=butter(4,[fband[0]*2./fs,fband[1]*2./fs],btype='band')normalization=np.sqrt((fband[1]-fband[0])/(fs/2))strain_H1_whitenbp=filtfilt(bb,ab,strain_H1_whiten)/normalizationstrain_L1_whitenbp=filtfilt(bb,ab,strain_L1_whiten)/normalization

ifmake_plots:# index into the strain time series for this time interval:indxt=np.where((time>=tevent-deltat)&(time<tevent+deltat))# pick a shorter FTT time interval, like 1/8 of a second:NFFT=int(fs/8)# and with a lot of overlap, to resolve short-time features:NOVL=int(NFFT*15./16)# and choose a window that minimizes "spectral leakage" # (https://en.wikipedia.org/wiki/Spectral_leakage)window=np.blackman(NFFT)# the right colormap is all-important! See:# http://matplotlib.org/examples/color/colormaps_reference.html# viridis seems to be the best for our purposes, but it's new; if you don't have it, you can settle for ocean.#spec_cmap='viridis'spec_cmap='ocean'# Plot the H1 spectrogram:plt.figure(figsize=(10,6))spec_H1,freqs,bins,im=plt.specgram(strain_H1[indxt],NFFT=NFFT,Fs=fs,window=window,noverlap=NOVL,cmap=spec_cmap,xextent=[-deltat,deltat])plt.xlabel('time (s) since '+str(tevent))plt.ylabel('Frequency (Hz)')plt.colorbar()plt.axis([-deltat,deltat,0,2000])plt.title('aLIGO H1 strain data near '+eventname)plt.savefig(eventname+'_H1_spectrogram.'+plottype)# Plot the L1 spectrogram:plt.figure(figsize=(10,6))spec_H1,freqs,bins,im=plt.specgram(strain_L1[indxt],NFFT=NFFT,Fs=fs,window=window,noverlap=NOVL,cmap=spec_cmap,xextent=[-deltat,deltat])plt.xlabel('time (s) since '+str(tevent))plt.ylabel('Frequency (Hz)')plt.colorbar()plt.axis([-deltat,deltat,0,2000])plt.title('aLIGO L1 strain data near '+eventname)plt.savefig(eventname+'_L1_spectrogram.'+plottype)