<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Matched Filtering Tutorial
This tutorial introduces matched filter in the context of detecting gravitational waves with LIGO. 
It borrows from a problem set developed by Christopher Berry.
Author
Step1:  Before we begin - download the LIGO data necessary for this notebook. This download may take a bit of time.
Problem 1) Detector Noise Curves
Read in the Hanford (H1) and Livingston (L1) amplitude spectral densities (ASD) data files (H1-ASD.txt, L1-ASD.txt).
Step2:  Problem 1a
From the ASDs, calculate the power spectral densities (PSDs) $S_n(f)$.
Step3:  Problem 1b
Plot the PSDs for Hanford and Livingston.
Step4:  Problem 1c
(no code) Looking at your plot of the PSDs, at approximately what frequencies are the detectors most sensitive?
Answer
Step5:  Answer
Step6:  You have also been given template waveforms $T(t)$ for a binary black hole, neutron star black hole, and binary neutron star coalescence
Step7:  Problem 2a
Plot the template waveforms $T(t)$. Try zooming into the last one or two seconds before merger (time = 0).
Step8:  Problem 2b
Plot the datastreams $h(t)$.
Step9:  Can you identify the presence of signals by eye?
Answer
Step10:  Problem 3b
Compute the Fourier transform the templates and plot the Fourier-transformed templates $\tilde{T}(f)$ as a function of frequency. Since the templates are all real, you only need to plot the positive Fourier frequencies and the magnitude of the Fourier transform.
Notice that the sampling rate is the same for all the templates (and all the timestreams), so the sampled frequencies are the same. (In other words, since the time bins are the same, the frequency bins are the same.)
Step12:  Problem 3c
Compute the matched filter signal-to-noise ratio timeseries between the templates and the datastreams. Start by computing the matched filter SNR timeseries between the BBH template and the "(a) datastream" in the Hanford detector (H1-a).
Note
Step13:  Problem 3d
Plot the SNR timeseries for each template in each of the (a), (b) and (c) datastreams, in each detector (18 combinations total). Do you see spikes in the SNR value as a function of the time shift? Which waveform signals are embedded within which datastreams? 
Note
Step14:  Which waveform signals are embedded within which datastreams? 
Answer
Step15:  Problem 4 (challenge)
Step16:  Problem 4b
Experiment with bandpassing the whitened data (e.g. using a Butterworth filter) to remove the noisiest frequencies, and remake some plots from Problem 3e.
Step17:  Problem 4c
Experiment with making time-frequency spectograms to visualize the data. Check out the Q-transform from gwpy.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
Explanation: Matched Filtering Tutorial
This tutorial introduces matched filter in the context of detecting gravitational waves with LIGO. 
It borrows from a problem set developed by Christopher Berry.
Author: Maya Fishbach
End of explanation
# Load the ASD data
dir = './Matched_Filter_Data/' #point to your local directory containing the data files
H1 = pd.read_csv(dir+'H1-ASD.txt', delim_whitespace=True)
L1 = pd.read_csv(dir+'L1-ASD.txt', delim_whitespace=True)
Explanation: Before we begin - download the LIGO data necessary for this notebook. This download may take a bit of time.
Problem 1) Detector Noise Curves
Read in the Hanford (H1) and Livingston (L1) amplitude spectral densities (ASD) data files (H1-ASD.txt, L1-ASD.txt).
End of explanation
H1['PSD'] = #complete
L1['PSD'] = #complete
Explanation: Problem 1a
From the ASDs, calculate the power spectral densities (PSDs) $S_n(f)$.
End of explanation
plt.loglog(#complete
plt.loglog(#complete
plt.ylabel(#complete
plt.xlabel(#complete
    
plt.legend(loc = 'best')
plt.xlim(10, 2000)
plt.ylim(10**-47,10**-37)
plt.show()
Explanation: Problem 1b
Plot the PSDs for Hanford and Livingston.
End of explanation
min_sampling_rate = #complete 
print('The minimum sampling rate is %.f'%(min_sampling_rate))
Explanation: Problem 1c
(no code) Looking at your plot of the PSDs, at approximately what frequencies are the detectors most sensitive?
Answer: 
Problem 1d
These ASDs have been calculated from data sampled in the time domain. What is the maximum frequency of the ASDs? What does this imply for the minimum sampling rate for the time domain data? Explain why.
Hint: You can find the actual sampling rate in the GWOSC Open Data paper.
End of explanation
# load the data
H1a = pd.read_csv(dir+'H1-a.txt', delim_whitespace=True)
H1b = pd.read_csv(dir+'H1-b.txt', delim_whitespace=True)
H1c = pd.read_csv(dir+'H1-c.txt', delim_whitespace=True)
L1a = pd.read_csv(dir+'L1-a.txt', delim_whitespace=True)
L1b = pd.read_csv(dir+'L1-b.txt', delim_whitespace=True)
L1c = pd.read_csv(dir+'L1-c.txt', delim_whitespace=True)
Explanation: Answer: 
Problem 2) datastreams and template waveforms
You have been given three noisy datastreams $h(t)$, each 600 seconds long, for the Hanford and Livingston detectors: *1-a.txt, *1-b.txt and *1-c.txt, where * is H for Hanford and L for Livingston.
These data streams all contain Gaussian noise drawn from the PSDs in Problem 1. They may additionally contain signals from compact binary inspirals.
End of explanation
#load the waveform templates
bbh = pd.read_csv(dir+'bbh_waveform.txt', delim_whitespace=True)
nsbh = pd.read_csv(dir+'nsbh_waveform.txt', delim_whitespace=True)
bns = pd.read_csv(dir+'bns_waveform.txt', delim_whitespace=True)
Explanation: You have also been given template waveforms $T(t)$ for a binary black hole, neutron star black hole, and binary neutron star coalescence: bbh-waveform.txt, nsbh-waveform.txt and bns-waveform.txt.
End of explanation
templates = [bbh, nsbh, bns]
for i, T in enumerate(templates):
    plt.plot(#complete
plt.legend(loc = 'best')
    
plt.xlabel(#complete
plt.ylabel(#complete
#plt.xlim(-0.5, 0) #zoom into the last half second before merger
plt.show()
Explanation: Problem 2a
Plot the template waveforms $T(t)$. Try zooming into the last one or two seconds before merger (time = 0).
End of explanation
H1_data = [H1a,H1b,H1c]
L1_data = [L1a,L1b,L1c]
labels = ['a','b','c']
# plot the strain data
f, axs = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=(16,12))
# label plots
axs[0,0].set_title('Hanford', fontsize=20)
axs[0,1].set_title('Livingston', fontsize=20)
axs[2,0].set_xlabel('time (s)', fontsize=15)
axs[2,1].set_xlabel('time (s)', fontsize=15)
axs[0,0].set_ylabel('strain', fontsize=15)
axs[1,0].set_ylabel('strain', fontsize=15)
axs[2,0].set_ylabel('strain', fontsize=15)
###complete
Explanation: Problem 2b
Plot the datastreams $h(t)$.
End of explanation
#Here we pad the template by placing zeros at the end
def pad_template(template):
    #complete
    
padded_bbh = pad_template(bbh)
padded_bns = pad_template(bns)
padded_nsbh = pad_template(nsbh)
Explanation: Can you identify the presence of signals by eye?
Answer:
Problem 3) Matched filter signal-to-noise ratio (SNR) time series
Using template waveforms $T(t)$, we will calculate and plot the matched filter SNR  as a function of the time shift $\tau$ between the template $T(t)$ and the data $h(t)$. We will repeat this for each of the three templates in each of the three datastreams, looking for peaks in the SNR time series that indicate the presence of gravitational-wave signals.
NOTE 1:
We will first need to Fourier transform both the template $T(t)$ and the data $h(t)$. Do not worry about normalizing your Fourier transforms, as we are only interested in qualitative results here.
NOTE 2:
To do the SNR calculation as a function of the time shift $\tau$, we will need the template to be the same length of the time as the data. There are two ways to achieve this: 
OPTION A: Artificially increase the length of the template waveform by zero-padding the template, such that the template is the same length of time as the data.
OPTION B: Look at a small chunk of data at a time, taking chunks of the data with time durations equal to the template. Perform the SNR calculation by "sliding" the template across the shortened data. 
Here, we will go with the method in which we zero-pad the template
Problem 3a
Zero-pad the BBH, NSBH, and BNS templates so that the templates are the same length as the H1-(a,b,c) & L1-(a,b,c) datastreams. The data from H1-(a,b,c) & L1-(a,b,c) are the same size, so the padding on each template is the same regardless of the datastream.
End of explanation
#first compute the corresponding frequency bins for the Fourier transform
sample_spacing = #complete
window_length = #complete
freq = np.fft.fftfreq(window_length, d=sample_spacing)
#now take the Fourier transforms of the templates
Tf_bbh = np.fft.fft(padded_bbh) #Fourier transform of the BBH padded template
Tf_nsbh = #complete #Fourier transform of the NSBH padded template
Tf_bns = #complete #Fourier transform of the BNS padded template
#select the positive frequencies and plot in log-log space
sel = freq > 0
labels = ['BBH', 'NSBH', 'BNS']
for i, Tf in enumerate([Tf_bbh, Tf_nsbh, Tf_bns]):
    plt.loglog(freq[sel], np.abs(Tf[sel]), label = labels[i])
plt.legend(loc = 'best')
plt.ylabel(#complete
plt.xlabel(#complete
plt.show()
Explanation: Problem 3b
Compute the Fourier transform the templates and plot the Fourier-transformed templates $\tilde{T}(f)$ as a function of frequency. Since the templates are all real, you only need to plot the positive Fourier frequencies and the magnitude of the Fourier transform.
Notice that the sampling rate is the same for all the templates (and all the timestreams), so the sampled frequencies are the same. (In other words, since the time bins are the same, the frequency bins are the same.)
End of explanation
def matched_filter_snr(template_Tf, freq_sampled, data_ht, psd):
     Compute the matched filter SNR time series
    template_Tf: waveform template in the frequency domain (Tf_bbh, Tf_nsbh, or Tf_bns from Problem 2)
    freq_sampled: frequency bins corresponding to template_Tf, (freq from Problem 2)
    data_ht: the datastream in the time domain, e.g. H1-a, L1-b
    psd: the dataset containing the PSD data from problem 1, either H1 or L1.
    returns: time, SNR values
    
    
    #evaluate the PSD at the sampled frequencies: 
    #because we are working only with real timeseries, the negative frequencies are the same as the positive frequencies
    
    psd_interp = np.interp(np.abs(freq_sampled), #complete
    
    #Fourier transform the strain data to the frequency domain
    data_hf = #complete
    
    #compute (unnormalized) matched filter SNR from template_Tf, data_hf, and psd_interp, in the frequency domain
    snr_freq = #complete
    
    #inverse Fourier transform snr_freq to the time domain
    snr_time = #complete
    
    return data_ht['time(s)'], np.abs(snr_time)
#compute the matched-filter SNR timeseries between the BBH template and H1-a
time, snr_time = matched_filter_snr(#complete
#plot the timeseries
plt.plot(time, snr_time)
plt.xlabel(#complete
plt.ylabel(#complete
plt.show()
Explanation: Problem 3c
Compute the matched filter signal-to-noise ratio timeseries between the templates and the datastreams. Start by computing the matched filter SNR timeseries between the BBH template and the "(a) datastream" in the Hanford detector (H1-a).
Note: we need to interpolate the PSDs so that we can evaluate them at the same frequencies as the templates and datastreams.
Recall that one can calculate the matched filter signal-to-noise ratio timeseries by taking the inverse Fourier transform of 
$$
\frac{\tilde{T}^*(f) \tilde{h}(f)}{S_n(f)}
$$
where $\tilde{T}^*(f)$ is the complex conjugate of the template in the frequency domain, $\tilde{h}(f)$ is the datastream in the frequency domain, and $S_n(f)$ is the PSD.
End of explanation
#complete
Explanation: Problem 3d
Plot the SNR timeseries for each template in each of the (a), (b) and (c) datastreams, in each detector (18 combinations total). Do you see spikes in the SNR value as a function of the time shift? Which waveform signals are embedded within which datastreams? 
Note: each datastream may contain multiple signals.
End of explanation
# first we need to calculate the time at which the SNR timeseries peaks. We pick the loudest peak.
def calculate_tau_maxSNR(#complete
    
    return tau
    
tau_a = calculate_tau_maxSNR(#complete #time at which the correct template peaks in the H1a datastream
tau_b = calculate_tau_maxSNR(#complete #time at which the correct template peaks in the H1b datastream
tau_c = calculate_tau_maxSNR(#complete #time at which the correct template peaks in the H1c datastream
times = H1a['time(s)']
#make the plot for H1-a
plt.title("Datastream a")
plt.plot(times, H1a['strain']) #plot the datastream
plt.plot(#complete #overplot the correct template with the correct timelag 
plt.xlabel(#complete
plt.ylabel(#complete
plt.xlim(0, 600)
plt.show()
#make the plot for H1-b
plt.title("Datastream b")
plt.plot(times, H1b['strain'])
plt.plot(#complete #overplot the correct template with the correct timelag 
plt.xlabel(#complete
plt.ylabel(#complete
plt.xlim(0, 600)
plt.show()
#make the plot for H1-c
plt.title("Datastream c")
plt.plot(times, H1c['strain'])
plt.plot(#complete #overplot the correct template with the correct timelag 
plt.xlabel(#complete
plt.ylabel(#complete
plt.xlim(0, 600)
plt.show()
Explanation: Which waveform signals are embedded within which datastreams? 
Answer:
Problem 3e
Using the results from part (d), overplot the signal time-domain waveforms over the datastreams that contain them. You can do this just for H1 datastreams, and just for one signal in each datastream. Your results should illustrate the power of matched filtering.  
Note: think about how the zero time $\tau = 0$ is defined.
End of explanation
#complete
Explanation: Problem 4 (challenge): visualize data
Experiment with whitening and bandpassing the data. This can make signals easier to see (for visualization purposes, as our eyes are not as good as matched filtering). As an example, have a look at Figure 3 of the LIGO & Virgo Data Analysis Guide.
Problem 4a
Remake the plots from Problem 3e using whitened data.
End of explanation
#complete
Explanation: Problem 4b
Experiment with bandpassing the whitened data (e.g. using a Butterworth filter) to remove the noisiest frequencies, and remake some plots from Problem 3e.
End of explanation
#complete
Explanation: Problem 4c
Experiment with making time-frequency spectograms to visualize the data. Check out the Q-transform from gwpy.
End of explanation
<END_TASK>