<jupyter_start><jupyter_text>Generate the estimated gravitational waveform and calculate the signal-to-noise for GW150914 with the Hanford Data<jupyter_code>import sys
!{sys.executable} -m pip install pycbc lalsuite ligo-common --no-cache-dir
#download the data
!curl -O -J https://raw.githubusercontent.com/ligo-cbc/binder/master/H-H1_LOSC_4_V2-1126259446-32.gwf
%matplotlib inline
import pylab
# Import the functions we need for later!
from pycbc.catalog import Merger
from pycbc.filter import highpass_fir, matched_filter
from pycbc.waveform import get_td_waveform
from pycbc.psd import welch, interpolate

# Read the Hanford data and remove low frequency content
h1 = Merger("GW150914").strain("H1")
h1 = highpass_fir(h1, 15, 16)

# Can you spot where the signal is beforehand?
pylab.plot(h1.sample_times, h1)
# Generate a waveform similar to GW150914 
# Change the parameters and see what happens to the waveform
# and the resulting SNR.

# Mass in Solar masses. 
m1 = 35.2
m2 = 34.0

# The intrinsic spin of each black hole
s1z = -0.228
s2z = -0.003

# The frequency to start generating the waveform
fstart = 15.0

hp, hc = get_td_waveform(approximant="SEOBNRv2",
                         mass1=m1, spin1z=s1z,
                         mass2=m2, spin2z=s2z,
                         f_lower=fstart,
                         delta_t=h1.delta_t)
pylab.plot(hp.sample_times, hp)
pylab.xlabel('Time (s)')
pylab.ylabel('Strain')

# Move the waveform so that the merge is about at the end
# This means that an SNR spike later on in the data lines up with this point
hp.resize(len(h1))
hp.roll(int(hp.start_time * hp.sample_rate))
# Estimate the noise spectrum
# Normally we use more data to estimate the psd, but this is illustrative
psd = interpolate(welch(h1), 1.0 / h1.duration)

pylab.loglog(psd.sample_frequencies, psd)
pylab.xlim(20, 1024)
pylab.xlabel('Frequency (Hz)')
# Calculate the complex Signal-to-noise. This is a complex vector
# because the signal could have any phase.
snr = matched_filter(hp, h1, psd=psd, low_frequency_cutoff=30.0)

# Remove regions corrupted by filter wraparound
snr = snr[len(snr) // 4: len(snr) * 3 // 4]

# Now you should be able to spot where the signal is!
pylab.plot(snr.sample_times, abs(snr))
pylab.ylabel('signal-to-noise')
pylab.xlabel('GPS Time (s)')
pylab.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09696710807347288, 'fraction_numerical': 0.03929944468175993, 'mean_word_length': 3.951374207188161, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '665011', 'license': 'MIT', 'n_tokens_mistral': 894, 'n_tokens_neox': 832, 'n_words': 256, 'path': 'PyCBC-Tutorials-master/examples/gw150914_snr.ipynb', 'repo_name': 'basuparth/ICERM_Workshop', 'size': '2341'}