Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added Time/Frequency Series SNR functionality #93

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -6,3 +6,8 @@ dist/
.hypothesis/*
docs/source/detectors_autogen.inc
docs/source/figures/*
*_gwstrain_trim.dat
Merger_data*.ipynb
*.hdf5
*eatmap*

6 changes: 3 additions & 3 deletions GWFish/detectors.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ ET:
lon: (9 + 25. / 60) * np.pi / 180.
opening_angle: np.pi / 3.
azimuth: 70.5674 * np.pi / 180.
psd_data: ET_psd.txt
psd_data: ET_full.txt
duty_factor: 0.85
detector_class: earthDelta
plotrange: 3, 1000, 1e-25, 1e-20
Expand Down Expand Up @@ -33,7 +33,7 @@ CE1:
lon: -119.4 * np.pi / 180.
opening_angle: np.pi / 2.
azimuth: 126. * np.pi / 180.
psd_data: CE1_psd.txt
psd_data: CE20km_full.txt
duty_factor: 0.85
detector_class: earthL
plotrange: 10, 1000, 1e-25, 1e-20
Expand All @@ -48,7 +48,7 @@ CE2:
lon: -119.4 * np.pi / 180.
opening_angle: np.pi / 2.
azimuth: 126. * np.pi / 180.
psd_data: CE2_psd.txt
psd_data: CE40km_full.txt
duty_factor: 0.85
detector_class: earthL
plotrange: 10, 1000, 1e-25, 1e-20
Expand Down
140 changes: 139 additions & 1 deletion GWFish/modules/utilities.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import pandas as pd
import yaml
from pathlib import Path
from lal import CreateREAL8TimeSeries, CreateREAL8Vector, DimensionlessUnit

DEFAULT_CONFIG = Path(__file__).parent.parent / 'detectors.yaml'
PSD_PATH = Path(__file__).parent.parent / 'detector_psd'
Expand Down Expand Up @@ -165,4 +166,141 @@ def get_snr(parameters, network, waveform_model):

return pd.DataFrame.from_dict(snrs, orient='index')


def make_fft_from_time_series(time_series_input, df, dt, title="Ines_Ludo"):
'''
Returns the FFT done through the lal library given a time series. Also returns the frequency array.

Parameters
----------
time_series_input : array
Time series data
df : float
Frequency step
dt : float
Time step
title : str, optional
Title of the time series

Returns
-------
tuple
FFT of the time series and the frequency array
'''
dims = len(time_series_input)
time_series = CreateREAL8Vector(dims)
time_series.data = time_series_input
ts = CreateREAL8TimeSeries(title, 1, 0, dt, DimensionlessUnit, dims)
ts.data = time_series
fft_dat = gw.fft.fft_lal_timeseries(ts, df).data.data
freq_range = np.linspace( 0, df * len(fft_dat), len(fft_dat) )

return fft_dat, freq_range

def _fd_phase_correction_and_output_format_from_stain_series(f_, hp, hc, geo_time = 1395964818):
'''
Prepares the polarizations for GWFish projection function. Combining
the functions "_fd_phase_correction_geocent_time", "_fd_gwfish_output_format" as in LALFD_Waveform class from waveforms.py module.

Parameters
----------
f_ : array
Frequency array
hp : array
Plus polarization
hc : array
Cross polarization
geo_time : int, optional
Geocentric time

Returns
-------
array
Polarizations in form (hp, hc)
'''
phi_in = np.exp( 1.j * (2 * f_ * np.pi * geo_time) ).T[0]
fft_dat_plus = phi_in*np.conjugate( hp )
fft_dat_cross = phi_in*np.conjugate( hc )

# GW Fish format for hfp and hfc
hfp = fft_dat_plus[:, np.newaxis]
hfc = fft_dat_cross[:, np.newaxis]
polarizations = np.hstack((hfp, hfc))

return polarizations

def get_snr(parameters, network, waveform_model = None):

waveform_class = gw.waveforms.LALFD_Waveform

nsignals = len(parameters)

# The SNR is then computed by taking the norm of the signal projected onto the detector
# and dividing by the noise of the detector
snrs = {}
for i in range(nsignals):
snr = {}
for detector in network.detectors:
data_params = {
'frequencyvector': detector.frequencyvector,
'f_ref': 50.
}
waveform_obj = waveform_class(waveform_model, parameters.iloc[i], data_params)
wave = waveform_obj()
t_of_f = waveform_obj.t_of_f
signal = gw.detection.projection(parameters.iloc[i], detector, wave, t_of_f)

snr[detector.name] = np.sqrt(np.sum(gw.detection.SNR(detector, signal)**2))

snr['network'] = np.sqrt(np.sum([snr[detector.name]**2 for detector in network.detectors]))
snrs['event_' + str(i)] = snr

return pd.DataFrame.from_dict(snrs, orient='index')

def get_SNR_from_series(f_in, hp, hc, network, parameters, geo_time = 1395964818, long_wavelength_approx = True):

'''
Given a set of parameters, polarizations, network, timevector and frequency array, returns the SNR associated to the signal

Parameters
----------
f_in : array
Frequency array on which to evaluate the signal
hp : array
Plus polarization without geocentric time phase corrections
hc : array
Cross polarization without geocentric time phase corrections
network : gw.detection.DetectorNetwork
Detector Network object
params : dict
Parameters of the event, needs to include ra, dec, psi
geo_time : int, optional
Geocentric time
long_wavelength_approx : bool, optional
Whether to use the long wavelength approximation or not

Returns
-------
float
Total signal-to-Noise Ratio
'''

polarizations = _fd_phase_correction_and_output_format_from_stain_series(f_in, hp, hc)
timevector = np.ones( len(f_in) ) * geo_time

series_data = (polarizations, timevector, f_in)

# SNR = get_snr(params, polarizations, detector, timevector, f_in, long_wavelength_approx)

polarizations, timevector, f_new = series_data

snrs_series = {}
for detector in network.detectors:
detector.frequencyvector = f_new
args = (parameters, detector, polarizations, timevector)
signal = gw.detection.projection(*args, long_wavelength_approx = long_wavelength_approx)
component_SNRs = gw.detection.SNR(detector, signal, frequencyvector=np.squeeze(f_new))
out_SNR = np.sqrt(np.sum(component_SNRs**2))
snrs_series[detector.name] = out_SNR

out_SNR = np.sqrt(np.sum([snrs_series[detector.name]**2 for detector in network.detectors]))
return out_SNR
66 changes: 66 additions & 0 deletions GWFish/modules/waveforms.py
Original file line number Diff line number Diff line change
Expand Up @@ -1016,3 +1016,69 @@ def plot(self, output_folder='./'):
plt.savefig(output_folder + 'psi_phenomD_zoomed.png')
plt.close()


class Ludo_Waveform(LALFD_Waveform):
"""
"""

def __init__(self, f, hp_func, hc_func, parameters):
self.hf_plus_func = hp_func
self.hf_cross_func = hc_func

self.gw_params = parameters
self.frequencyvector = f

def _update_frequency_range_indices(self):
self.idx_low = int(self.f_min / self.delta_f)
self.idx_high = int(self.f_max / self.delta_f)

def _lal_fd_strain_adjust_frequency_range(self):
""" Frequency array starts from zero, so we need to mask some frequencies """
self._update_frequency_range_indices()
self.hf_cross_out = self._lal_hf_cross.data.data[self.idx_low:self.idx_high+1]
self.hf_plus_out = self._lal_hf_plus.data.data[self.idx_low:self.idx_high+1]

def _lal_fd_phase_correction_by_epoch_and_df(self):
""" This correction is also done in Bilby after calling SimInspiralFD """
# BORIS: weird Bilby correction
dt = 1. / self.delta_f + (self._lal_hf_plus.epoch.gpsSeconds +
self._lal_hf_plus.epoch.gpsNanoSeconds * 1e-9)
self.hf_plus_out *= np.exp(
-1j * 2 * np.pi * dt * self.frequencyvector)
self.hf_cross_out *= np.exp(
-1j * 2 * np.pi * dt * self.frequencyvector)

def _hf_postproccessing_SimInspiralFD(self):
self._lal_fd_strain_adjust_frequency_range()
self._lal_fd_phase_correction_by_epoch_and_df()

def _hf_postproccessing_SimInspiralCFDWS(self):
self.hf_plus_out, self.hf_cross_out = self._lal_hf_plus.data.data, self._lal_hf_cross.data.data

def _fd_phase_correction_geocent_time(self):
""" Add initial 2pi*f*tc - phic - pi/4 to phase """
phi_in = np.exp(1.j*(2*self.frequencyvector*np.pi*self.gw_params['geocent_time']))

hfp = phi_in * np.conjugate(self.hf_plus_out) # it's already multiplied by the phase
hfc = phi_in * np.conjugate(self.hf_cross_out)

return hfp, hfc

def _fd_gwfish_output_format(self, hfp, hfc):

hfp = hfp[:, np.newaxis]
hfc = hfc[:, np.newaxis]

polarizations = np.hstack((hfp, hfc))

return polarizations

def calculate_frequency_domain_strain(self):

self.hf_plus_out = self.hf_plus_func(self.frequencyvector)
self.hf_cross_out = self.hf_cross_func(self.frequencyvector)

hfp, hfc = self._fd_phase_correction_geocent_time()
polarizations = self._fd_gwfish_output_format(hfp, hfc)

self._frequency_domain_strain = polarizations
541 changes: 541 additions & 0 deletions series_snr_tutorial.ipynb

Large diffs are not rendered by default.