<jupyter_start><jupyter_text>Read in the strain (only H1)<jupyter_code>with h5py.File('../data/strain/H1_2017_4096.hdf5') as file:
    raw_strain = np.array(file['strain/Strain'])<jupyter_output><empty_output><jupyter_text>Read in and crop the waveforms<jupyter_code>with h5py.File('../data/waveforms/waveforms_3s_0100_0300.h5') as file:
    
    # Empty list for the final, cropped waveforms
    waveforms = []
    
    # Read in the waveforms and the indices of the ones where the generation failed
    raw_waveforms = np.array(file['waveforms'])
    failed = np.array(file['failed'])
    
    # Remove failed waveforms
    raw_waveforms = np.array([raw_waveforms[i] for i in range(len(raw_waveforms)) if i not in failed])
    
    # Remove zero-padding from the waveforms
    for i in range(len(raw_waveforms)):
        waveform = raw_waveforms[i]
        for j in range(len(waveform)):
            if waveform[j] != 0:
                start = j
                break
        for j in range(len(waveform))[::-1]:
            if waveform[j] != 0:
                end = j
                break
        waveforms.append(waveform[start:end])<jupyter_output><empty_output><jupyter_text>Calculate SNR in the 'classical' way<jupyter_code>waveform = np.random.choice(waveforms)
start = int(np.random.uniform(0, len(raw_strain)-len(waveform)))
noise = raw_strain[start:start+len(waveform)]

signal = noise # + 0.0*waveform

plt.plot(signal)
plt.plot(noise)
plt.show()
power_signal = np.mean(np.square(np.abs(waveform)))
power_noise = np.mean(np.square(np.abs(noise)))

print('SNR = {}'.format(power_signal / power_noise))<jupyter_output>SNR = 2.08418193856e-06<jupyter_text>Calculate SNR using the 'matched filtering' way<jupyter_code>from pycbc.filter import matched_filter
from pycbc.types import TimeSeries
from pycbc.psd import welch, interpolate

raw_strain_ts = TimeSeries(raw_strain, delta_t=1.0/4096)
signal_ts = TimeSeries(signal, delta_t=1.0/4096)
noise_ts = TimeSeries(noise, delta_t=1.0/4096)
waveform_ts = TimeSeries(waveform, delta_t=1.0/4096)

# Calculate the complex (two-phase SNR)
psd = interpolate(welch(raw_strain_ts), delta_f=waveform_ts.delta_f)
snr = matched_filter(noise_ts, signal_ts, psd=psd)

snr = snr[len(snr) / 4: len(snr) * 3 / 4]

import pylab
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.12728026533996684, 'fraction_numerical': 0.02902155887230514, 'mean_word_length': 4.245652173913044, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '457628', 'license': 'Apache-2.0', 'n_tokens_mistral': 950, 'n_tokens_neox': 881, 'n_words': 173, 'path': 'quantum gravity/convwave/src/jupyter/snr_experiments.ipynb', 'repo_name': 'HiteshDhola/Artificial-Intelligence-Deep-Learning-Machine-Learning-Tutorials', 'size': '2412'}