<jupyter_start><jupyter_text>Analysis of LIGO/Virgo strain data GW150914<jupyter_code>import os
import numpy as np
import matplotlib.pyplot as plt
import json
import h5py
from pathlib import Path

from scipy.signal import tukey

from gwpy.timeseries import TimeSeries
import pycbc.psd<jupyter_output><empty_output><jupyter_text>Obtain dataUse gwpy to download the data for GW150914. Consider only the "H1" data for now; "L1" should be treated similarly.<jupyter_code>t_event = 1126259462.391  # GPS time of coalescence

T = 8.0  # number of seconds to analyze in a segment
T_psd = 1024.0  # number of seconds of data for estimating PSD
T_buffer = 2.0  # buffer time after the event to include
h1_event = TimeSeries.fetch_open_data('H1', t_event + T_buffer - T, t_event + T_buffer, cache=True)

# To calculate the PSD, take data before the event segment
h1_psd = TimeSeries.fetch_open_data('H1', t_event + T_buffer - T - T_psd, t_event + T_buffer - T, cache=True)
# Check sampling frequency
f_s = len(h1_event)/T
f_s
plot = h1_event.plot()
plt.show()
plot = h1_psd.plot()
plt.show()<jupyter_output><empty_output><jupyter_text>Convert to pycbc TimeSeries objects.<jupyter_code>h1_event = h1_event.to_pycbc()
h1_psd = h1_psd.to_pycbc()<jupyter_output><empty_output><jupyter_text>PSDWe use the "welch" method. Default settings for this method are to use a Hann window, and in each frequency bin to take the median over all segments. However, we will use a Tukey window, with 0.4 s roll-off, and no overlap, as is described in the LALInference paper (https://arxiv.org/abs/1409.7215).  Window functionThis window will be used both for PSD estimation and for the event segment.<jupyter_code>roll_off = 0.4
alpha = 2 * roll_off / T

w = tukey(int(T * f_s), alpha)
plt.plot(h1_event.sample_times, w)
plt.xlabel('time (s)')
plt.ylabel('w')
plt.show()<jupyter_output><empty_output><jupyter_text>PSD estimation<jupyter_code># psd = pycbc.psd.estimate.welch(h1_psd, seg_len=int(T * f_s), seg_stride=int(T * f_s / 2), window='hann', avg_method='median')
psd = pycbc.psd.estimate.welch(h1_psd, seg_len=int(T * f_s), seg_stride=int(T * f_s), window=w, avg_method='median')
plt.figure(figsize=(15,10))
plt.plot(psd.sample_frequencies, psd)
plt.yscale('log')
plt.xlabel('frequency (Hz)')
plt.ylabel('1/Hz')
plt.title('Estimated PSD for GW150914')
plt.xlim((20,1024))
plt.ylim((1e-47, 1e-39))
plt.show()<jupyter_output><empty_output><jupyter_text>Aside: BayesWave PSDFor comparison to estimated PSD.<jupyter_code># I downloaded this data file from GWOSC
psd_bayeswave = np.loadtxt('../data/psds/GWTC1_GW150914_PSDs.dat.txt')

# Add a final point at 1024 Hz
psd_bayeswave = np.append(psd_bayeswave, np.reshape(psd_bayeswave[-1], (1,3)), axis=0)
psd_bayeswave[-1,0] = 1024.0
psd_bayeswave = pycbc.psd.from_numpy_arrays(psd_bayeswave[:,0], psd_bayeswave[:,1], 8193, 1/8, 20.0)
#plt.figure(figsize=(15,11))
plt.plot(psd_bayeswave.sample_frequencies, psd_bayeswave, label='BayesWave')
plt.plot(psd.sample_frequencies, psd, label='Estimated')
plt.yscale('log')
plt.xlim((950,1024))
plt.ylim((1e-47, 1e-39))
plt.legend()
plt.xlabel('frequency (Hz)')
plt.ylabel('1/Hz')
plt.show()<jupyter_output><empty_output><jupyter_text>The BayesWave PSD has a weird dip just after 1000 Hz. I think this is an effect of using a spline to fit the big jump at slightly lower frequencies. This gives a whitened waveform that has a lot of extra power just beyond 1000 Hz, so I will use the PSD estimated using the welch method.  EventWindow the 8 s segment that contains the event, take the FFT, and scale by the inverse ASD.<jupyter_code>h1_event_FD = (h1_event * w).to_frequencyseries()
h1_event_FD_whitened = h1_event_FD / psd ** 0.5
h1_event_FD_whitened.data[1000]
plt.figure(figsize=(15,10))
plt.plot(h1_event_FD_whitened.sample_frequencies, h1_event_FD_whitened.real())
plt.xlabel('frequency (Hz)')
#plt.xlim((995,1005))
plt.xscale('log')
plt.xlim((10,2048))
plt.show()<jupyter_output><empty_output><jupyter_text>Zero below 20 Hz and truncate above 1024 Hz. This is all that is needed for the likelihood integral.<jupyter_code>f_min = 20.0
f_max = 1024.0  # New f_max. Previously it was 2048 Hz.

h1_event_FD_whitened[:int(f_min / h1_event_FD_whitened.delta_f)] = 0.0
h1_event_FD_whitened = h1_event_FD_whitened[:int(f_max / h1_event_FD_whitened.delta_f)+1]
plt.plot(h1_event_FD_whitened.sample_frequencies, h1_event_FD_whitened.real())
plt.xlabel('frequency (Hz)')
plt.show()<jupyter_output><empty_output><jupyter_text>Apply a time shift to shift the merger to $t=0$. This should make the event visible in the Fourier transform.<jupyter_code>h1_event_FD_whitened = h1_event_FD_whitened.cyclic_time_shift(T_buffer)
plt.plot(h1_event_FD_whitened.sample_frequencies, h1_event_FD_whitened.real())
plt.xlabel('frequency (Hz)')
plt.xscale('log')
plt.xlim((20,1024))
plt.show()<jupyter_output><empty_output><jupyter_text>NormalizationThe noise is not quite normalized to unit variance in the real and imaginary parts.<jupyter_code># Compute the variance of the real part of the whitened FD strain data.

np.var(h1_event_FD_whitened.numpy().real)<jupyter_output><empty_output><jupyter_text>There are two issues here:* We need to account for $\delta f$ for the DFT, the 1-sided PSD, and the real and imaginary parts.* We need to account for the window function that we used. This requires multiplication by a "window factor." Since the signal is much shorter than $T$, the windowing only affects the amount of noise in the data, not the signal.<jupyter_code>delta_f = h1_event_FD_whitened.delta_f

normalization_1 = 1 / np.sqrt(4 * delta_f)
window_factor = np.sum(w ** 2) / (T * f_s)

normalization = normalization_1 * np.sqrt(window_factor)
normalization ** 2<jupyter_output><empty_output><jupyter_text>They match, more or less.<jupyter_code>normalized_FD = h1_event_FD_whitened.numpy() / normalization
print(np.var(normalized_FD.real), np.var(normalized_FD.imag))<jupyter_output>1.017263252150474 0.9981325458433057<jupyter_text>The point of this is that when we train our network, we need to generate noise, which is added to whitened waveforms. Ultimately we plan to evaluate the network on the whitened event strain data, which has been windowed before FFTing, as above. Thus the noise that is added during training should have a standard deviation given by `normalization` in each of the real and imaginary parts.  Save PSD to file<jupyter_code>event_dir = Path('../data/events/GW150914')
event_dir.mkdir(parents=True, exist_ok=True)

psd.save(event_dir / 'PSD_H1.txt')<jupyter_output><empty_output><jupyter_text>L1 detectorDo the same calculation, save to file.<jupyter_code>l1_psd_data = TimeSeries.fetch_open_data('L1', t_event + T_buffer - T - T_psd, t_event + T_buffer - T, cache=True)
l1_psd_data = l1_psd_data.to_pycbc()
l1_psd = pycbc.psd.estimate.welch(l1_psd_data, seg_len=int(T * f_s), seg_stride=int(T * f_s), window=w, avg_method='median')
plt.plot(l1_psd.sample_frequencies, l1_psd)
plt.yscale('log')
plt.xlabel('frequency (Hz)')
plt.ylabel('1/Hz')
plt.title('Estimated PSD for GW150914 L1')
plt.xlim((20,1024))
plt.ylim((1e-47, 1e-39))
plt.show()
l1_event_data = TimeSeries.fetch_open_data('L1', t_event + T_buffer - T, t_event + T_buffer, cache=True)
l1_event_data = l1_event_data.to_pycbc()
l1_event_FD = (l1_event_data * w).to_frequencyseries()
l1_event_FD_whitened = l1_event_FD / l1_psd ** 0.5

l1_event_FD_whitened[:int(f_min / delta_f)] = 0.0
l1_event_FD_whitened = l1_event_FD_whitened[:int(f_max / delta_f)+1]

l1_event_FD_whitened = l1_event_FD_whitened.cyclic_time_shift(T_buffer)
plt.plot(l1_event_FD_whitened.sample_frequencies, l1_event_FD_whitened.real())
plt.xlabel('frequency (Hz)')
plt.xscale('log')
plt.xlim((20,1024))
plt.title('L1 whitened strain')
plt.show()<jupyter_output><empty_output><jupyter_text>Compare H1 and L1 PSDs<jupyter_code>plt.figure(figsize=(15,11))
plt.plot(psd.sample_frequencies, psd, label='H1')
plt.plot(l1_psd.sample_frequencies, l1_psd, label='L1')
plt.yscale('log')
plt.xlabel('frequency (Hz)')
plt.ylabel('1/Hz')
plt.title('Estimated PSDs for GW150914')
plt.xlim((20,1024))
plt.ylim((1e-47, 1e-38))
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>These PSDs have a fair bit more structure than the design sensitivity PSDs. Since we train our reduced basis to represent whitened waveforms, this may mean we need more basis elements to obtain a good representation.<jupyter_code>l1_psd.save(event_dir / 'PSD_L1.txt')<jupyter_output><empty_output><jupyter_text>Save<jupyter_code># Whitened FD strain data

with h5py.File(event_dir / 'strain_FD_whitened.hdf5', 'w') as f:
    f.create_dataset('sample_frequencies', data=np.array(h1_event_FD_whitened.sample_frequencies.numpy()))
    f.create_dataset('H1', data=h1_event_FD_whitened)
    f.create_dataset('L1', data=l1_event_FD_whitened)
# Additional information

with open(event_dir / 'event_info.json', 'w') as f:
    json.dump(dict(
        event='GW150914',
        detectors=['H1', 'L1'],
        t_event=t_event,
        f_min=f_min,
        f_max=f_max,
        T=T,
        T_psd=T_psd,
        T_buffer=T_buffer,
        window_factor=window_factor,
        window='tukey',
        roll_off=roll_off
        ), f, indent=4)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11774245864417775, 'fraction_numerical': 0.040112444588604176, 'mean_word_length': 6.099002302379125, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 68, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '656868', 'license': 'MIT', 'n_tokens_mistral': 3871, 'n_tokens_neox': 3562, 'n_words': 817, 'path': 'notebooks/GW150914_data.ipynb', 'repo_name': 'iphysresearch/lfi-gw', 'size': '9249'}