<jupyter_start><jupyter_text>Define functions we will need later on<jupyter_code>def apply_psd(signal_t, psd, sampling_rate=4096, apply_butter=False):
    """
    Take a signal in the time domain, and a precalculated Power Spectral
    Density, and color the signal according to the given PSD.

    Args:
        signal_t: A signal in time domain (i.e. a 1D numpy array)
        psd: A Power Spectral Density, e.g. calculated from the detector noise.
            Should be a function: psd(frequency)
        sampling_rate: Sampling rate of signal_t
        apply_butter: Whether or not to apply a Butterworth filter to the data.

    Returns: color_signal_t, the colored signal in the time domain.
    """

    # First set some parameters for computing power spectra
    signal_size = len(signal_t)
    delta_t = 1 / sampling_rate

    # Go into Fourier (frequency) space: signal_t -> signal_f
    frequencies = np.fft.rfftfreq(signal_size, delta_t)
    signal_f = np.fft.rfft(signal_t)

    # Divide by the given Power Spectral Density (PSD)
    # This is the 'whitening' = actually adding color
    color_signal_f = signal_f / (np.sqrt(psd(frequencies) / delta_t / 2))

    # Go back into time space: color_signal_f -> color_signal_t
    color_signal_t = np.fft.irfft(color_signal_f, n=signal_size)

    # In case we want to use a Butterworth-filter, here's how to do it:
    if apply_butter:

        # Define cut-off frequencies for the filter
        f_low = 42
        f_high = 800

        # Calculate Butterworth-filter and normalization
        numerator, denominator = butter(4, [f_low*2/4096, f_high*2/4096],
                                        btype="bandpass")
        normalization = np.sqrt((f_high - f_low) / (sampling_rate / 2))

        # Apply filter and normalize
        color_signal_t = filtfilt(numerator, denominator, color_signal_t)
        color_signal_t = color_signal_t / normalization

    return color_signal_t


# -----------------------------------------------------------------------------


def get_psd(real_strain, sampling_rate=4096):
    """
    Take a detector recording and calculate the Power Spectral Density (PSD).

    Args:
        real_strain: The detector recording to be used.
        sampling_rate: The sampling rate (in Hz) of the recording

    Returns:
        psd: The Power Spectral Density of the detector recordings
    """

    # Define some constants
    nfft = 2 * sampling_rate  # Bigger values yield better resolution?

    # Use matplotlib.mlab to calculate the PSD from the real strain
    power_spectrum, frequencies = mlab.psd(real_strain,
                                           NFFT=nfft,
                                           Fs=sampling_rate)

    # Interpolate it linearly, so we can re-sample the spectrum arbitrarily
    psd = interp1d(frequencies, power_spectrum)

    return psd<jupyter_output><empty_output><jupyter_text>Read in the strain files for all events<jupyter_code># Path to the directory where all data is stored
data_path = '../data'

strains = dict()

for event in ['GW150914', 'GW151226', 'GW170104']:
    
    strains[event] = dict()
    
    with h5py.File(os.path.join(data_path, 'strain', '{}_H1_STRAIN_4096.h5'.format(event)), 'r') as file:
        strains[event]['H1'] = np.array(file['strain/Strain'])
    with h5py.File(os.path.join(data_path, 'strain', '{}_L1_STRAIN_4096.h5'.format(event)), 'r') as file:
        strains[event]['L1'] = np.array(file['strain/Strain'])
psds = dict()

for event in ['GW150914', 'GW151226', 'GW170104']:
    
    psds[event] = dict()
    psds[event]['H1'] = get_psd(strains[event]['H1'])
    psds[event]['L1'] = get_psd(strains[event]['L1'])<jupyter_output><empty_output><jupyter_text>Read in and whiten the waveform files<jupyter_code>waveforms = dict()
for dist in ['0100_0300', '0250_0500', '0400_0800', '0700_1200', '1000_1700']:
    with h5py.File(os.path.join(data_path, 'waveforms', 'waveforms_3s_{}_testing.h5'.format(dist)), 'r') as file:
        waveforms_ = np.array(file['waveforms'])
        failed_ = np.array(file['failed'])
        waveforms[dist] = [_ for i, _ in enumerate(waveforms_) if i not in failed_]
white_waveforms = dict()
for event in ['GW150914', 'GW151226', 'GW170104']:
    white_waveforms[event] = dict()
    for dist in ['0100_0300', '0250_0500', '0400_0800', '0700_1200', '1000_1700']:
    
        white_waveforms[event][dist] = dict()
        white_waveforms[event][dist]['H1'] = []
        white_waveforms[event][dist]['L1'] = []
    
        for waveform in waveforms[dist]:
            white_waveforms[event][dist]['H1'].append(apply_psd(waveform, psds[event]['H1']))
            white_waveforms[event][dist]['L1'].append(apply_psd(waveform, psds[event]['L1']))<jupyter_output><empty_output><jupyter_text>Calculate the SNRs<jupyter_code>stds = {'GW150914': {'H1': 0.98780110486117378, 'L1': 0.98664230166344458},
        'GW151226': {'H1': 0.99035140654295228, 'L1': 0.98776116540338088},
        'GW170104': {'H1': 0.99212202082054124, 'L1': 0.98727435617370363}}
snrs = dict()
for event in ['GW150914', 'GW151226', 'GW170104']:
    snrs[event] = dict()
    for dist in ['0100_0300', '0250_0500', '0400_0800', '0700_1200', '1000_1700']:
        
        snrs[event][dist] = dict()
        snrs[event][dist]['H1'] = []
        snrs[event][dist]['L1'] = []
        
        for i, white_waveform in enumerate(white_waveforms[event][dist]['H1']):
            maximum = np.max(np.abs(white_waveform))
            snrs[event][dist]['H1'].append(maximum / stds[event]['H1'])
            if maximum == 0:
                print(event, dist, i)
        for i, white_waveform in enumerate(white_waveforms[event][dist]['L1']):
            maximum = np.max(np.abs(white_waveform))
            snrs[event][dist]['L1'].append(maximum / stds[event]['L1'])
            if maximum == 0:
                print(event, dist, i)
for event in ['GW150914', 'GW151226', 'GW170104']:
    display(HTML('<h3>SNRs for {}</h3>'.format(event)))
    rows = []
    for dist in ['0100_0300', '0250_0500', '0400_0800', '0700_1200', '1000_1700']:
        median_H1 = '{:.2f}'.format(np.nanmedian(snrs[event][dist]['H1']))
        median_L1 = '{:.2f}'.format(np.nanmedian(snrs[event][dist]['L1']))
        min_H1    = '{:.2f}'.format(np.nanmin(snrs[event][dist]['H1']))
        min_L1    = '{:.2f}'.format(np.nanmin(snrs[event][dist]['L1']))
        max_H1    = '{:.2f}'.format(np.nanmax(snrs[event][dist]['H1']))
        max_L1    = '{:.2f}'.format(np.nanmax(snrs[event][dist]['L1']))
        rows.append([dist, median_H1, min_H1, max_H1, median_L1, min_L1, max_L1])
    display(HTML(tabulate(rows, tablefmt='html', 
                          headers=['Distances', 'Median H1', 'Minimum H1', 'Maximum H1', 'Median L1', 'Minimum L1', 'Maximum L1'])))
    # print(tabulate(rows, tablefmt='latex', 
    #                      headers=['Distances', 'Median H1', 'Minimum H1', 'Maximum H1', 'Median L1', 'Minimum L1', 'Maximum L1']))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15297814597914583, 'fraction_numerical': 0.06984716469075847, 'mean_word_length': 3.0995316159250588, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '554405', 'license': 'Apache-2.0', 'n_tokens_mistral': 2716, 'n_tokens_neox': 2385, 'n_words': 539, 'path': 'quantum gravity/convwave/src/jupyter/get_snr_distributions_from_testing.ipynb', 'repo_name': 'HiteshDhola/Artificial-Intelligence-Deep-Learning-Machine-Learning-Tutorials', 'size': '7001'}