home/comp_snr_ts.md

@@ -129,3 +129,58 @@ print(optimal_snr)
 ```
 
 ![Screen_Shot_2021-03-11_at_4.42.14_PM](uploads/270beaa00acd95de52983771f318d7f4/Screen_Shot_2021-03-11_at_4.42.14_PM.png)
+
+
+- find the snr time series in the whitened scale in *gwpy*
+```
+psd_whitened = pycbc_interpolate(psd_whitened, d_pycbc.delta_f)
+opt_snr_inj = pycbc.filter.matchedfilter.sigma(d_pycbc, psd=psd_whitened, low_frequency_cutoff=10, high_frequency_cutoff=1024)
+print('optimal SNR of the noise: {:.2f}'.format(opt_snr_inj))
+d_pycbc_dummy = d_pycbc.copy()
+peak_time_index = np.argmax(d_pycbc.data)
+peak_time = d_pycbc.sample_times.data[peak_time_index]
+peak_window_start_index = peak_time_index - int(d_pycbc.sample_rate * 1)
+peak_window_end_index = peak_time_index + int(d_pycbc.sample_rate * 1)
+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_whitened, low_frequency_cutoff=10, high_frequency_cutoff=1024)
+plt.plot(snr_timeseries_n.sample_times, np.abs(snr_timeseries_n.data))
+#plt.ylim(0, 10)
+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.52.15_PM](uploads/8f8c2ddbe698769a63b44bbe23aafa6e/Screen_Shot_2021-03-11_at_4.52.15_PM.png)
+
+- find the snr time series in the whitened scale in *pycbc*
+
+```
+whitened_pycbc = d_pycbc.whiten(1, 1)
+d_pycbc = whitened_pycbc
+psd_whitened = whitened_pycbc.psd(1)
+
+psd_whitened = pycbc_interpolate(psd_whitened, d_pycbc.delta_f)
+opt_snr_inj = pycbc.filter.matchedfilter.sigma(d_pycbc, psd=psd_whitened, low_frequency_cutoff=10, high_frequency_cutoff=1024)
+print('optimal SNR of the noise: {:.2f}'.format(opt_snr_inj))
+d_pycbc_dummy = d_pycbc.copy()
+peak_time_index = np.argmax(d_pycbc.data)
+peak_time = d_pycbc.sample_times.data[peak_time_index]
+peak_window_start_index = peak_time_index - int(d_pycbc.sample_rate * 1)
+peak_window_end_index = peak_time_index + int(d_pycbc.sample_rate * 1)
+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_whitened, low_frequency_cutoff=10, high_frequency_cutoff=1024)
+plt.plot(snr_timeseries_n.sample_times, np.abs(snr_timeseries_n.data))
+#plt.ylim(0, 10)
+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.56.55_PM](uploads/03ddfe236b7eb2238360a99f7ab21b77/Screen_Shot_2021-03-11_at_4.56.55_PM.png)
\ No newline at end of file