In-Class Assignment 19#

Credit: Gravitational Wave Open Science Center

Learning Objectives#

  • Use GWpy to analyze open-source LIGO data

  • Practice applying timeseries analysis techniques: FFT, whitening, bandpasses

  • Convert whitened data to audio WAV file

To Begin#

Load some imports including GWpy - a python package for gravitational-wave astrophysics.

We will be utilizing the TimeSeries tools described in Documentation.

import requests, os
import matplotlib.pyplot as plt
import numpy as np
%config InlineBackend.figure_format = 'retina'

try:
    from gwpy.timeseries import TimeSeries
except:
    ! pip install -q "gwpy==3.0.9"
    ! pip install -q "matplotlib==3.7.3"
    ! pip install -q "astropy==6.1.4"
    from gwpy.timeseries import TimeSeries

a. - Set parameters for the signal we want to explore#

The two first parameters we need to set are:

  • GPS time of the data of interest (Learn more)

  • Detector could be H1 (LIGO Hanford), L1 (LIGO Livingston), or V1 (Virgo)

Example GPS times for different events are:

You might try some of these examples times in the H1 detector:

t0 = 1126259462.4    # -- GW150914
t0 = 1187008882.4    # -- GW170817
t0 = 933200215       # -- Loud hardware injection
t0 = 1132401286.33   # -- Koi Fish Glitch
  1. Choose and set a GPS time

  2. Choose a detector

  3. Round the GPS time to integer then query for the strain using the integer GPS time, and the detector using TimeSeries.fetch_open_data.

  4. Plot the resulting strain as a function of time.

# a result here

b. - Whiten and band-pass the data#

  • Whitening is a process that re-weights a signal, so that all frequency bins have a nearly equal amount of noise.

  • A band-pass filter uses both a low frequency cutoff and a high frequency cutoff, and only passes signals in the frequency band between these values.

See also:

  1. Createn a new array of the whitened strain data

  2. Apply a bandpass to the whitened data from 30 to 400 Hz.

  3. Plot the whitened and bandpassed data from (our t0 GPS time - 0.2 seconds, to our t0 GPS time + 0.1 seconds.)

Can we identify the inspiral, merger, and or ringdown of the post-processed data?

# b result here

c. - Plot a q-transform of the data#

A Q-transform plot shows how a signal’s frequency changes with time.

  • The x-axis shows time

  • The y-axis shows frequency

  • The color scale shows the amount of “energy” or “signal power” in each time-frequency pixel.

A parameter called “Q” refers to the quality factor. A higher quality factor corresponds to a larger number of cycles in each time-frequency pixel.

  1. Using a dt=1, use GWpy to apply a [q-transform] with an out segment of \(t_{0}\pm dt\).

  2. Plot the resulting transformation using the GWpy native plotting tool. Make sure y-scale is set to log.

# c results here

Some helper functions for d#

# function to keep the data within integer limits, and write to wavfile:
def write_wavfile(filename,fs,data):
    d = np.int16(data/np.max(np.abs(data)) * 32767 * 0.9)
    wavfile.write(filename,int(fs), d)

deltat_sound = 2.                     # seconds around the event

# function that shifts frequency of a band-passed signal
def reqshift(data,fshift=100,sample_rate=4096):
    """Frequency shift the signal by constant
    """
    x = np.fft.rfft(data)
    T = len(data)/float(sample_rate)
    df = 1.0/T
    nbins = int(fshift/df)
    # print T,df,nbins,x.real.shape
    y = np.roll(x.real,nbins) + 1j*np.roll(x.imag,nbins)
    y[0:nbins]=0.
    z = np.fft.irfft(y)
    return z

d. - Convert the data to audio#

  1. Define the time array of our final bandpassed and whitened strain data. We only want the data not the units (.value). Use tevent = t0. Compute the indices of interested and write the wav file for our data.

  2. Listen to the data wav file.

  3. Finally, shift the data using the function line, write the wav file of the shifted data.

  4. Listen to the result.

Can you hear the chirp heard around the world?

# make wav (sound) files from the whitened+BP data, +-2s around the event.
import scipy
from scipy.io import wavfile
from IPython.display import Audio

Listening to the Whitened-Bandpassed Strain Data in the original Frequency Range#

#time = bp_data.times
#time = time.value
#tevent = 
#fs = 4096 # sampling freq

# index into the strain time series for this time interval:
#indxd = np.where((time >= tevent-deltat_sound) & (time < tevent+deltat_sound))

# write the files:
#write_wavfile(detector+"_whitenbp.wav",int(fs), strain[indxd])
# listen to the wav file here
#fna = detector+"_whitenbp.wav"
#print(fna)
#Audio(fna)

Listening to the frequency shifted Whitened-Bandpassed Strain Data#

# using our helper function, shift the data and make the new wav file
#fshift = 400.
#strain_shifted = reqshift(bp_data,fshift=fshift,sample_rate=fs)
#write_wavfile(detector+"_whitenbp_shifted.wav",int(fs), strain_shifted[indxd])
# listen to the wav file here
#fna = detector+"_whitenbp_shifted.wav"
#print(fna)
#Audio(fna)