In [13]:
import sys
from pycbc.catalog import Merger
from pycbc import frame
import pylab
%matplotlib inline
From public data¶
In [4]:
m = Merger('GW150914')
data = {}
for ifo in ['H1', 'L1']:
data[ifo] = m.strain(ifo)
In [5]:
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()
In [6]:
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()
In [7]:
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()
In [8]:
# 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()
In [9]:
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 [10]:
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()
In [11]:
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()
In [12]:
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()
From LDR: proprietary data¶
In [51]:
# for GW190814
merger_time = 1249852257.0
Read in frames directly¶
In [53]:
frame_paths = {
'H1' : '/hdfs/frames/O3/hoft_C01/H1/H-H1_HOFT_C01-12498/H-H1_HOFT_C01-1249849344-4096.gwf',
'L1' : '/hdfs/frames/O3/hoft_C01/L1/L-L1_HOFT_C01-12498/L-L1_HOFT_C01-1249849344-4096.gwf',
'V1' : '/hdfs/frames/O3/V1Online/V-V1Online-12498/V-V1Online-1249850000-2000.gwf,/hdfs/frames/O3/V1Online/V-V1Online-12498/V-V1Online-1249852000-2000.gwf'
}
channel_names = {
'H1' : 'H1:DCS-CALIB_STRAIN_CLEAN_C01',
'L1' : 'L1:DCS-CALIB_STRAIN_CLEAN_C01',
'V1' : 'V1:Hrec_hoft_16384Hz_Gated'
}
frame_types = {'H1': 'H1_HOFT_C01', 'L1': 'L1_HOFT_C01','V1': 'V1Online'}
In [50]:
frame_data = {}
for ifo in frame_paths:
ifo_frame_paths = frame_paths[ifo]
for idx, pth in enumerate(ifo_frame_paths.split(',')):
print("Reading {0} for {1}".format(pth, ifo))
data = frame.read_frame(pth, channel_names[ifo])
print("...read in data with length {0}".format(
data.delta_t * len(data)))
if idx == 0:
all_data = data
if idx > 0:
all_data.append_zeros( len(data) )
all_data[len(all_data) - len(data):len(all_data)] = data
frame_data[ifo] = all_data
Query the LDR server¶
In [55]:
x = frame.query_and_read_frame('H1_HOFT_C01', 'H1:DCS-CALIB_STRAIN_CLEAN_C01',
merger_time - 15, merger_time + 15)
print(len(x) * x.delta_t)
In [56]:
frame_data = {}
for ifo in channel_names:
frame_data[ifo] = frame.query_and_read_frame(frame_types[ifo], channel_names[ifo],
merger_time - 2048, merger_time + 2048)
print("Length of data read in for {0}: {1} seconds".format(ifo,
len(frame_data[ifo]) * frame_data[ifo].delta_t))
In [68]:
# Whiten the data
whitened_data = {}
for ifo in frame_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_data[ifo] = frame_data[ifo].whiten(4, 4)
zoom = whitened_data[ifo].time_slice(merger_time - 1.5, merger_time + 0.5)
pylab.plot(zoom.sample_times, zoom, label=ifo)
pylab.ylabel('Whitened Strain')
pylab.xlabel('Time (s)')
pylab.legend()
pylab.show()
In [72]:
for ifo in whitened_data:
#if 'V1' in ifo: continue
print("Plotting QScan for {0}".format(ifo))
# We'll choose a tighter zoom here.
zoom = whitened_data[ifo].time_slice(merger_time - 5, merger_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=(32, 32),
frange=(20, 512),
)
pylab.figure(figsize=[15, 3])
pylab.pcolormesh(times, freqs, power**0.5)
pylab.xlim(merger_time - 1.5, merger_time + 0.3)
pylab.title(ifo)
pylab.yscale('log')
pylab.show()
In [73]:
for ifo in whitened_data:
#if 'V1' in ifo: continue
print("Plotting QScan for {0}".format(ifo))
# We'll choose a tighter zoom here.
zoom = whitened_data[ifo].time_slice(merger_time - 5, merger_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=(32, 32),
frange=(20, 512),
)
pylab.figure(figsize=[15, 3])
pylab.pcolormesh(times, freqs, power**0.5)
pylab.xlim(merger_time - 1.5, merger_time + 0.3)
pylab.title(ifo)
pylab.yscale('log')
pylab.show()
In [74]:
for ifo in whitened_data:
#if 'V1' in ifo: continue
print("Plotting QScan for {0}".format(ifo))
# We'll choose a tighter zoom here.
zoom = whitened_data[ifo].time_slice(merger_time - 5, merger_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=(64, 64),
frange=(20, 512),
)
pylab.figure(figsize=[15, 3])
pylab.pcolormesh(times, freqs, power**0.5)
pylab.xlim(merger_time - 1.5, merger_time + 0.3)
pylab.title(ifo)
pylab.yscale('log')
pylab.show()
In [ ]: