home/comp_snr_ts.md

@@ -78,3 +78,30 @@ plt.close()
 ```
 ![Screen_Shot_2021-03-11_at_4.25.13_PM](uploads/f238606e0dd7f52ff3bdccf3ed8dfe67/Screen_Shot_2021-03-11_at_4.25.13_PM.png)
 
+-- find the snr time series in pycbc
+```
+# let's do the same thing in pycbc object
+d_pycbc = d_gwpy.to_pycbc()
+psd = d_pycbc.psd(1)
+psd = pycbc_interpolate(psd, d_pycbc.delta_f)
+opt_snr_inj = pycbc.filter.matchedfilter.sigma(d_pycbc, psd=psd, low_frequency_cutoff=10, high_frequency_cutoff=1024)
+print('optimal SNR of the noise: {:.2f}'.format(opt_snr_inj))
+w = d_pycbc.whiten(1, 1)
+d_pycbc = d_pycbc.copy()
+peak_time_index = np.argmax(d_pycbc.data)
+peak_time = d_pycbc.sample_times.data[peak_time_index]
+window = 1
+peak_window_start_index = peak_time_index - int(d_pycbc.sample_rate * window)
+peak_window_end_index = peak_time_index + int(d_pycbc.sample_rate * window)
+d_pycbc_dummy = d_pycbc_dummy.cyclic_time_shift(np.round(d_pycbc.sample_times.data[-1] - peak_time, 3))
+snr_timeseries_n = pycbc.filter.matchedfilter.matched_filter(d_pycbc_dummy, d_pycbc, psd=psd, low_frequency_cutoff=10, high_frequency_cutoff=1024)
+plt.plot(snr_timeseries_n.sample_times, np.abs(snr_timeseries_n.data))
+plt.axvline(peak_time - window, c='black', ls='--')
+plt.axvline(peak_time + window, c='black', ls='--')
+plt.savefig('test.png')
+plt.close()
+optimal_snr = np.abs(snr_timeseries_n.data)[peak_window_start_index: peak_window_end_index].max()
+print(optimal_snr)
+
+```
+![Screen_Shot_2021-03-11_at_4.34.42_PM](uploads/63e5b574a76a6b0ec1a8bf3b7a927684/Screen_Shot_2021-03-11_at_4.34.42_PM.png)