From cbd4fad534b4680a3a3c4da81368f126ecdc0df0 Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Sun, 1 Mar 2026 23:52:02 +0100 Subject: [PATCH 01/10] Add inject_glitch method Fix spelling mistakes --- bilby/gw/detector/interferometer.py | 75 +++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 94267d30c..9125dab36 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -392,6 +392,81 @@ def check_signal_duration(self, parameters, raise_error=True): else: logger.warning(msg) + def adjust_optimal_snr_and_onset_time_of_glitch(self, glitch, target_glitch_snr, glitch_onset_time): + """ Scale the time-domain strain data of the glitch to achieve a target optimal SNR. + Roll the glitch in such a way that the maxima coincides with the glitch_onset_time. + + Parameters + ========== + glitch: numpy array + Time-domain strain data of the glitch to scale. + target_glitch_snr: float + Desired optimal SNR of the glitch. + glitch_onset_time: float + GPS time of the occurrence of the glitch + + Returns + ======= + scaled_glitch: numpy array + frequency-domain strain of the glitch rescaled so that its optimal SNR matches + `target_glitch_snr` and rolled in such a way thatthe maxima coincides with the glitch_onset_time. + """ + n_full = len(self.strain_data.time_array) + padded_glitch = np.zeros(n_full) + padded_glitch[:len(glitch)] = glitch + peak_idx = np.argmax(np.abs(glitch)) + relative_glitch_onset_time = glitch_onset_time - self.start_time + target_idx = int(relative_glitch_onset_time * self.sampling_frequency) + padded_glitch = np.roll(padded_glitch, target_idx - peak_idx) + glitch_fft, _ = utils.nfft(padded_glitch, sampling_frequency=self.sampling_frequency) + temporary_snr_squared = np.real(self.optimal_snr_squared(signal = glitch_fft)) + scaled_glitch = glitch_fft * target_glitch_snr / np.sqrt(temporary_snr_squared) + return scaled_glitch + + + def check_matched_filter_glitch_snr(self, glitch): + """ Compute filter SNR of a glitch in detector. + + Parameters + ========== + glitch: numpy array + Time-domain strain data of the glitch. + """ + matched_filter_glitch_snr = self.matched_filter_snr(signal = glitch) + logger.info("Matched filter SNR of the glitch: {}".format(matched_filter_glitch_snr)) + + + def inject_glitch(self, glitch, glitch_onset_time, glitch_sampling_frequency, glitch_snr): + """ Inject a glitch into the interferometer data. + + Parameters + ========== + glitch: numpy array + Time-domain strain data of the glitch to inject. + glitch_onset_time: float + GPS time at which to inject the glitch. The maxima of the glitch array will coincide with the glitch onset time. + glitch_sampling_frequency: float + Sampling frequency of the glitch. + glitch_snr: float + The glitch will be rescaled in such a way that it's optimal SNR matches with the glitch_snr. + """ + + glitch_time_array = np.arange(0, len(glitch), 1) / glitch_sampling_frequency + glitch_duration = glitch_time_array[-1] + + if glitch_sampling_frequency != self.sampling_frequency: + _temporary_time_array = np.arange(0, glitch_duration, 1 / self.sampling_frequency) + glitch = np.interp(_temporary_time_array, glitch_time_array, glitch) + glitch_time_array = _temporary_time_array + + scaled_glitch = self.adjust_optimal_snr_and_onset_time_of_glitch(glitch, glitch_snr, glitch_onset_time) + self.strain_data.frequency_domain_strain += scaled_glitch + logger.info("Injected a glitch in: {}".format(self.name)) + logger.info("Optimal SNR of the glitch: {}".format(glitch_snr)) + self.check_matched_filter_glitch_snr(scaled_glitch) + + + def inject_signal(self, parameters, injection_polarizations=None, waveform_generator=None, raise_error=True): """ General signal injection method. From fb3a07f5e18b32986386b5ff39616f86e1f15b7d Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Mon, 2 Mar 2026 00:10:49 +0100 Subject: [PATCH 02/10] Add inject glitch method --- bilby/gw/detector/interferometer.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 9125dab36..4e5bb4d21 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -393,7 +393,7 @@ def check_signal_duration(self, parameters, raise_error=True): logger.warning(msg) def adjust_optimal_snr_and_onset_time_of_glitch(self, glitch, target_glitch_snr, glitch_onset_time): - """ Scale the time-domain strain data of the glitch to achieve a target optimal SNR. + """ Scale the time-domain strain data of the glitch to achieve a target optimal SNR. Roll the glitch in such a way that the maxima coincides with the glitch_onset_time. Parameters @@ -419,23 +419,21 @@ def adjust_optimal_snr_and_onset_time_of_glitch(self, glitch, target_glitch_snr, target_idx = int(relative_glitch_onset_time * self.sampling_frequency) padded_glitch = np.roll(padded_glitch, target_idx - peak_idx) glitch_fft, _ = utils.nfft(padded_glitch, sampling_frequency=self.sampling_frequency) - temporary_snr_squared = np.real(self.optimal_snr_squared(signal = glitch_fft)) + temporary_snr_squared = np.real(self.optimal_snr_squared(signal=glitch_fft)) scaled_glitch = glitch_fft * target_glitch_snr / np.sqrt(temporary_snr_squared) return scaled_glitch - def check_matched_filter_glitch_snr(self, glitch): - """ Compute filter SNR of a glitch in detector. + """ Compute matched filter SNR of a glitch in detector. Parameters ========== glitch: numpy array - Time-domain strain data of the glitch. + frequency-domain strain data of the glitch. """ - matched_filter_glitch_snr = self.matched_filter_snr(signal = glitch) + matched_filter_glitch_snr = self.matched_filter_snr(signal=glitch) logger.info("Matched filter SNR of the glitch: {}".format(matched_filter_glitch_snr)) - def inject_glitch(self, glitch, glitch_onset_time, glitch_sampling_frequency, glitch_snr): """ Inject a glitch into the interferometer data. @@ -444,11 +442,12 @@ def inject_glitch(self, glitch, glitch_onset_time, glitch_sampling_frequency, gl glitch: numpy array Time-domain strain data of the glitch to inject. glitch_onset_time: float - GPS time at which to inject the glitch. The maxima of the glitch array will coincide with the glitch onset time. + GPS time at which to inject the glitch. The maxima of the glitch array will + coincide with the glitch onset time. glitch_sampling_frequency: float Sampling frequency of the glitch. glitch_snr: float - The glitch will be rescaled in such a way that it's optimal SNR matches with the glitch_snr. + The glitch will be rescaled in such a way that it's optimal SNR matches with the glitch_snr. """ glitch_time_array = np.arange(0, len(glitch), 1) / glitch_sampling_frequency @@ -459,14 +458,12 @@ def inject_glitch(self, glitch, glitch_onset_time, glitch_sampling_frequency, gl glitch = np.interp(_temporary_time_array, glitch_time_array, glitch) glitch_time_array = _temporary_time_array - scaled_glitch = self.adjust_optimal_snr_and_onset_time_of_glitch(glitch, glitch_snr, glitch_onset_time) + scaled_glitch = self.adjust_optimal_snr_and_onset_time_of_glitch(glitch, glitch_snr, glitch_onset_time) self.strain_data.frequency_domain_strain += scaled_glitch logger.info("Injected a glitch in: {}".format(self.name)) logger.info("Optimal SNR of the glitch: {}".format(glitch_snr)) self.check_matched_filter_glitch_snr(scaled_glitch) - - def inject_signal(self, parameters, injection_polarizations=None, waveform_generator=None, raise_error=True): """ General signal injection method. From 34220c99eac2fa3acdacc17a457b3872801c7526 Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Mon, 2 Mar 2026 00:20:17 +0100 Subject: [PATCH 03/10] Remove matched_filter_snr function --- bilby/gw/detector/interferometer.py | 42 ++++++++++++----------------- 1 file changed, 17 insertions(+), 25 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 4e5bb4d21..adc0c481d 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -392,58 +392,48 @@ def check_signal_duration(self, parameters, raise_error=True): else: logger.warning(msg) - def adjust_optimal_snr_and_onset_time_of_glitch(self, glitch, target_glitch_snr, glitch_onset_time): - """ Scale the time-domain strain data of the glitch to achieve a target optimal SNR. - Roll the glitch in such a way that the maxima coincides with the glitch_onset_time. + def adjust_optimal_snr_and_injection_time_of_glitch(self, glitch, glitch_snr, glitch_injection_time): + """ Scale the time-domain strain data of the glitch so that it's optimal SNR is equal to the the glitch_snr + Roll the glitch in such a way that the maxima coincides with the glitch_injection_time. Parameters ========== glitch: numpy array Time-domain strain data of the glitch to scale. - target_glitch_snr: float + glitch_snr: float Desired optimal SNR of the glitch. - glitch_onset_time: float + glitch_injection_time: float GPS time of the occurrence of the glitch Returns ======= scaled_glitch: numpy array frequency-domain strain of the glitch rescaled so that its optimal SNR matches - `target_glitch_snr` and rolled in such a way thatthe maxima coincides with the glitch_onset_time. + `target_glitch_snr` and rolled in such a way thatthe maxima coincides with the glitch_injection_time. """ n_full = len(self.strain_data.time_array) padded_glitch = np.zeros(n_full) padded_glitch[:len(glitch)] = glitch + # Roll glitch peak_idx = np.argmax(np.abs(glitch)) - relative_glitch_onset_time = glitch_onset_time - self.start_time - target_idx = int(relative_glitch_onset_time * self.sampling_frequency) + relative_glitch_injection_time = glitch_injection_time - self.start_time + target_idx = int(relative_glitch_injection_time * self.sampling_frequency) padded_glitch = np.roll(padded_glitch, target_idx - peak_idx) + # Perform fft glitch_fft, _ = utils.nfft(padded_glitch, sampling_frequency=self.sampling_frequency) temporary_snr_squared = np.real(self.optimal_snr_squared(signal=glitch_fft)) - scaled_glitch = glitch_fft * target_glitch_snr / np.sqrt(temporary_snr_squared) + scaled_glitch = glitch_fft * glitch_snr / np.sqrt(temporary_snr_squared) return scaled_glitch - def check_matched_filter_glitch_snr(self, glitch): - """ Compute matched filter SNR of a glitch in detector. - - Parameters - ========== - glitch: numpy array - frequency-domain strain data of the glitch. - """ - matched_filter_glitch_snr = self.matched_filter_snr(signal=glitch) - logger.info("Matched filter SNR of the glitch: {}".format(matched_filter_glitch_snr)) - - def inject_glitch(self, glitch, glitch_onset_time, glitch_sampling_frequency, glitch_snr): + def inject_glitch(self, glitch, glitch_injection_time, glitch_sampling_frequency, glitch_snr): """ Inject a glitch into the interferometer data. Parameters ========== glitch: numpy array Time-domain strain data of the glitch to inject. - glitch_onset_time: float - GPS time at which to inject the glitch. The maxima of the glitch array will - coincide with the glitch onset time. + glitch_injection_time: float + The maxima of the time domain strain of the glitch will conicide with the glitch_injection_time. glitch_sampling_frequency: float Sampling frequency of the glitch. glitch_snr: float @@ -458,10 +448,12 @@ def inject_glitch(self, glitch, glitch_onset_time, glitch_sampling_frequency, gl glitch = np.interp(_temporary_time_array, glitch_time_array, glitch) glitch_time_array = _temporary_time_array - scaled_glitch = self.adjust_optimal_snr_and_onset_time_of_glitch(glitch, glitch_snr, glitch_onset_time) + scaled_glitch = self.adjust_optimal_snr_and_injection_time_of_glitch(glitch, glitch_snr, glitch_injection_time) self.strain_data.frequency_domain_strain += scaled_glitch logger.info("Injected a glitch in: {}".format(self.name)) logger.info("Optimal SNR of the glitch: {}".format(glitch_snr)) + logger.info("Matched filter SNR of the glitch: {}".format(self.matched_filter_snr(signal=scaled_glitch))) + self.check_matched_filter_glitch_snr(scaled_glitch) def inject_signal(self, parameters, injection_polarizations=None, From 5b42e88dceb3788aa802e8909da9427c19bb459f Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Mon, 2 Mar 2026 00:22:22 +0100 Subject: [PATCH 04/10] Move to logger.info --- bilby/gw/detector/interferometer.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index adc0c481d..d67a6221b 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -454,8 +454,6 @@ def inject_glitch(self, glitch, glitch_injection_time, glitch_sampling_frequency logger.info("Optimal SNR of the glitch: {}".format(glitch_snr)) logger.info("Matched filter SNR of the glitch: {}".format(self.matched_filter_snr(signal=scaled_glitch))) - self.check_matched_filter_glitch_snr(scaled_glitch) - def inject_signal(self, parameters, injection_polarizations=None, waveform_generator=None, raise_error=True): """ General signal injection method. From 2545f45cc30560ee7346babcf6179d6d3dbfdf69 Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Mon, 2 Mar 2026 00:34:07 +0100 Subject: [PATCH 05/10] Add a check for the duration --- bilby/gw/detector/interferometer.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index d67a6221b..ac4dc63e3 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -435,7 +435,7 @@ def inject_glitch(self, glitch, glitch_injection_time, glitch_sampling_frequency glitch_injection_time: float The maxima of the time domain strain of the glitch will conicide with the glitch_injection_time. glitch_sampling_frequency: float - Sampling frequency of the glitch. + Sampling frequency of the input glitch array. glitch_snr: float The glitch will be rescaled in such a way that it's optimal SNR matches with the glitch_snr. """ @@ -443,6 +443,12 @@ def inject_glitch(self, glitch, glitch_injection_time, glitch_sampling_frequency glitch_time_array = np.arange(0, len(glitch), 1) / glitch_sampling_frequency glitch_duration = glitch_time_array[-1] + if glitch_duration > self.duration: + raise ValueError( + "Glitch duration ({} seconds) is longer than the duration of the time domain strain of the detector " + "({} seconds).".format(glitch_duration, self.duration) + ) + if glitch_sampling_frequency != self.sampling_frequency: _temporary_time_array = np.arange(0, glitch_duration, 1 / self.sampling_frequency) glitch = np.interp(_temporary_time_array, glitch_time_array, glitch) From 3bd3df0645df9f627a7dea1aff75a6c3c95bbf5d Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Mon, 2 Mar 2026 01:02:21 +0100 Subject: [PATCH 06/10] Example --- .../injection_examples/add_glitch_into_ifo.py | 80 +++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 examples/gw_examples/injection_examples/add_glitch_into_ifo.py diff --git a/examples/gw_examples/injection_examples/add_glitch_into_ifo.py b/examples/gw_examples/injection_examples/add_glitch_into_ifo.py new file mode 100644 index 000000000..4dbbc90e6 --- /dev/null +++ b/examples/gw_examples/injection_examples/add_glitch_into_ifo.py @@ -0,0 +1,80 @@ +import bilby +import gengli +import matplotlib.pyplot as plt +import numpy as np + +# To be installed from https://git.ligo.org/melissa.lopez/gengli +# or use glitch generator of your choice. + +duration = 4.0 +sampling_frequency = 1024.0 +glitch_snr = 12 +injection_parameters = dict( + mass_1=36.0, + mass_2=29.0, + a_1=0.4, + a_2=0.3, + tilt_1=0.5, + tilt_2=1.0, + phi_12=1.7, + phi_jl=0.3, + luminosity_distance=2000.0, + theta_jn=0.4, + psi=2.659, + phase=1.3, + geocent_time=1126259642.413, + ra=1.375, + dec=-1.2108, +) + +bilby.core.utils.random.seed(23) + +# Fixed arguments passed into the source model +waveform_arguments = dict( + waveform_approximant="IMRPhenomXPHM", + reference_frequency=50.0, + minimum_frequency=20.0, +) + +# Create the waveform_generator using a LAL BinaryBlackHole source function +# the generator will convert all the parameters +waveform_generator = bilby.gw.WaveformGenerator( + duration=duration, + sampling_frequency=sampling_frequency, + frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, + parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments=waveform_arguments, +) + +# Set up interferometers. In this case we'll use two interferometers +# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design +# sensitivity +ifos = bilby.gw.detector.InterferometerList(["H1", "L1"]) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=injection_parameters["geocent_time"] - 2, +) + +ifos.inject_signal( + waveform_generator=waveform_generator, parameters=injection_parameters +) + +generator = gengli.glitch_generator("L1") +glitch = generator.get_glitch(seed=200) * 1e-23 + +# Inject glitch +ifos[0].inject_glitch( + glitch, + glitch_injection_time=ifos[0].start_time + 2, + glitch_sampling_frequency=sampling_frequency, + glitch_snr=glitch_snr, +) + +fig, axes = plt.subplots(2, 1) +ax = axes[0] +ax.plot(ifos[0].time_array, ifos[0].time_domain_strain) + +ax = axes[1] +ax.plot(np.arange(len(glitch)) / sampling_frequency, glitch) +fig.savefig("./time_domain_strain.pdf") From b82a451d5d0525bd30f81311eb1aa83c7da1b065 Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Mon, 2 Mar 2026 01:03:18 +0100 Subject: [PATCH 07/10] Update description Time shift work in progress Time shifts fixed --- bilby/gw/detector/interferometer.py | 145 ++++++++++++------ .../injection_examples/add_glitch_into_ifo.py | 4 +- 2 files changed, 98 insertions(+), 51 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index ac4dc63e3..ec2130602 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -43,7 +43,8 @@ class Interferometer(object): minimum_frequency = PropertyAccessor('strain_data', 'minimum_frequency') maximum_frequency = PropertyAccessor('strain_data', 'maximum_frequency') frequency_mask = PropertyAccessor('strain_data', 'frequency_mask') - frequency_domain_strain = PropertyAccessor('strain_data', 'frequency_domain_strain') + frequency_domain_strain = PropertyAccessor( + 'strain_data', 'frequency_domain_strain') time_domain_strain = PropertyAccessor('strain_data', 'time_domain_strain') def __init__(self, name, power_spectral_density, minimum_frequency, maximum_frequency, length, latitude, longitude, @@ -108,10 +109,14 @@ def __repr__(self): 'maximum_frequency={}, length={}, latitude={}, longitude={}, elevation={}, ' \ 'xarm_azimuth={}, yarm_azimuth={}, xarm_tilt={}, yarm_tilt={})' \ .format(self.name, self.power_spectral_density, float(self.strain_data.minimum_frequency), - float(self.strain_data.maximum_frequency), float(self.geometry.length), - float(self.geometry.latitude), float(self.geometry.longitude), - float(self.geometry.elevation), float(self.geometry.xarm_azimuth), - float(self.geometry.yarm_azimuth), float(self.geometry.xarm_tilt), + float(self.strain_data.maximum_frequency), float( + self.geometry.length), + float(self.geometry.latitude), float( + self.geometry.longitude), + float(self.geometry.elevation), float( + self.geometry.xarm_azimuth), + float(self.geometry.yarm_azimuth), float( + self.geometry.xarm_tilt), float(self.geometry.yarm_tilt)) def set_strain_data_from_gwpy_timeseries(self, time_series): @@ -280,7 +285,8 @@ def antenna_response(self, ra, dec, time, psi, mode): """ if mode in ["plus", "cross", "x", "y", "breathing", "longitudinal"]: - polarization_tensor = get_polarization_tensor(ra, dec, time, psi, mode) + polarization_tensor = get_polarization_tensor( + ra, dec, time, psi, mode) return three_by_three_matrix_contraction(self.geometry.detector_tensor, polarization_tensor) elif mode == self.name: return 1 @@ -342,7 +348,8 @@ def get_detector_response(self, waveform_polarizations, parameters, frequencies= dt_geocent = parameters['geocent_time'] - self.strain_data.start_time dt = dt_geocent + time_shift - signal_ifo[mask] = signal_ifo[mask] * np.exp(-1j * 2 * np.pi * dt * frequencies) + signal_ifo[mask] = signal_ifo[mask] * \ + np.exp(-1j * 2 * np.pi * dt * frequencies) signal_ifo[mask] *= self.calibration_model.get_calibration_factor( frequencies, prefix='recalib_{}_'.format(self.name), **parameters @@ -371,7 +378,8 @@ def check_signal_duration(self, parameters, raise_error=True): if ("mass_1" not in parameters) and ("mass_2" not in parameters): if raise_error: - raise AttributeError("Unable to check signal duration as mass not given") + raise AttributeError( + "Unable to check signal duration as mass not given") else: return @@ -392,13 +400,13 @@ def check_signal_duration(self, parameters, raise_error=True): else: logger.warning(msg) - def adjust_optimal_snr_and_injection_time_of_glitch(self, glitch, glitch_snr, glitch_injection_time): + def adjust_optimal_snr(self, frequency_domain_strain, target_snr): """ Scale the time-domain strain data of the glitch so that it's optimal SNR is equal to the the glitch_snr Roll the glitch in such a way that the maxima coincides with the glitch_injection_time. Parameters ========== - glitch: numpy array + glitch_strain: numpy array Time-domain strain data of the glitch to scale. glitch_snr: float Desired optimal SNR of the glitch. @@ -411,21 +419,12 @@ def adjust_optimal_snr_and_injection_time_of_glitch(self, glitch, glitch_snr, gl frequency-domain strain of the glitch rescaled so that its optimal SNR matches `target_glitch_snr` and rolled in such a way thatthe maxima coincides with the glitch_injection_time. """ - n_full = len(self.strain_data.time_array) - padded_glitch = np.zeros(n_full) - padded_glitch[:len(glitch)] = glitch - # Roll glitch - peak_idx = np.argmax(np.abs(glitch)) - relative_glitch_injection_time = glitch_injection_time - self.start_time - target_idx = int(relative_glitch_injection_time * self.sampling_frequency) - padded_glitch = np.roll(padded_glitch, target_idx - peak_idx) - # Perform fft - glitch_fft, _ = utils.nfft(padded_glitch, sampling_frequency=self.sampling_frequency) - temporary_snr_squared = np.real(self.optimal_snr_squared(signal=glitch_fft)) - scaled_glitch = glitch_fft * glitch_snr / np.sqrt(temporary_snr_squared) - return scaled_glitch - - def inject_glitch(self, glitch, glitch_injection_time, glitch_sampling_frequency, glitch_snr): + temporary_snr_squared = np.real( + self.optimal_snr_squared(signal=frequency_domain_strain)) + return frequency_domain_strain * target_snr / np.sqrt(temporary_snr_squared) + + def inject_glitch(self, glitch_snr, glitch_parameters=None, glitch_time_domain_strain=None, + glitch_sample_times=None, glitch_injection_time=None, glitch_waveform_generator=None): """ Inject a glitch into the interferometer data. Parameters @@ -440,25 +439,64 @@ def inject_glitch(self, glitch, glitch_injection_time, glitch_sampling_frequency The glitch will be rescaled in such a way that it's optimal SNR matches with the glitch_snr. """ - glitch_time_array = np.arange(0, len(glitch), 1) / glitch_sampling_frequency - glitch_duration = glitch_time_array[-1] + if glitch_parameters is None: + glitch_parameters = {} + for key, val in [('ra', 0.0), ('dec', 0.0), ('psi', 0.0)]: + glitch_parameters.setdefault(key, val) - if glitch_duration > self.duration: + if glitch_time_domain_strain is None and glitch_waveform_generator is None: raise ValueError( - "Glitch duration ({} seconds) is longer than the duration of the time domain strain of the detector " - "({} seconds).".format(glitch_duration, self.duration) - ) + "inject_glitch needs one of glitch_waveform_generator or " + "glitch_time_domain_strain.") + elif glitch_time_domain_strain is not None: + glitch_sampling_frequency = 1.0 / (glitch_sample_times[1] - glitch_sample_times[0]) + # Resample the glitch to match with the interferometer's sampling frequency + if glitch_sampling_frequency != self.sampling_frequency: + temp_sample_times = np.arange( + len(self.time_array)) / self.sampling_frequency + glitch_time_domain_strain = np.interp( + temp_sample_times, glitch_sample_times, glitch_time_domain_strain) + glitch_sample_times = temp_sample_times + + # Padding the glitch to have the correct size for the fft. + padded_glitch = np.zeros(len(self.time_array)) + padded_glitch[:len(glitch_time_domain_strain) + ] += glitch_time_domain_strain + + # Since the maxima of the glitch is not at 0. + glitch_parameters['geocent_time'] = glitch_injection_time - \ + (glitch_sample_times[np.argmax(padded_glitch)]) + + # Convert to frequency domain glitch + glitch_frequency_domain_strain, _ = utils.nfft( + padded_glitch, self.sampling_frequency) + + # Adjust SNR + glitch_frequency_domain_strain = self.adjust_optimal_snr( + glitch_frequency_domain_strain, glitch_snr) + injection_polarizations = { + self.name: glitch_frequency_domain_strain} + + # Usual inject method + self.inject_signal_from_waveform_polarizations(parameters=glitch_parameters, + injection_polarizations=injection_polarizations) - if glitch_sampling_frequency != self.sampling_frequency: - _temporary_time_array = np.arange(0, glitch_duration, 1 / self.sampling_frequency) - glitch = np.interp(_temporary_time_array, glitch_time_array, glitch) - glitch_time_array = _temporary_time_array + elif glitch_waveform_generator is not None: + glitch_frequency_domain_strain = glitch_waveform_generator.frequency_domain_strain( + glitch_parameters) + glitch_frequency_domain_strain = self.adjust_optimal_snr( + glitch_frequency_domain_strain, glitch_snr) + injection_polarizations = { + self.name: glitch_frequency_domain_strain} + self.inject_signal_from_waveform_polarizations(parameters=glitch_parameters, + injection_polarizations=injection_polarizations) - scaled_glitch = self.adjust_optimal_snr_and_injection_time_of_glitch(glitch, glitch_snr, glitch_injection_time) - self.strain_data.frequency_domain_strain += scaled_glitch logger.info("Injected a glitch in: {}".format(self.name)) logger.info("Optimal SNR of the glitch: {}".format(glitch_snr)) - logger.info("Matched filter SNR of the glitch: {}".format(self.matched_filter_snr(signal=scaled_glitch))) + logger.info("Matched filter SNR of the glitch: {}".format( + self.matched_filter_snr(signal=glitch_frequency_domain_strain))) + + return glitch_time_domain_strain def inject_signal(self, parameters, injection_polarizations=None, waveform_generator=None, raise_error=True): @@ -558,7 +596,8 @@ def inject_signal_from_waveform_polarizations(self, parameters, injection_polari 'Injecting signal outside segment, start_time={}, merger time={}.' .format(self.strain_data.start_time, parameters['geocent_time'])) - signal_ifo = self.get_detector_response(injection_polarizations, parameters) + signal_ifo = self.get_detector_response( + injection_polarizations, parameters) self.strain_data.frequency_domain_strain += signal_ifo self.meta_data['optimal_SNR'] = ( @@ -568,8 +607,10 @@ def inject_signal_from_waveform_polarizations(self, parameters, injection_polari self.meta_data['parameters'] = parameters logger.info("Injected signal in {}:".format(self.name)) - logger.info(" optimal SNR = {:.2f}".format(self.meta_data['optimal_SNR'])) - logger.info(" matched filter SNR = {:.2f}".format(self.meta_data['matched_filter_SNR'])) + logger.info(" optimal SNR = {:.2f}".format( + self.meta_data['optimal_SNR'])) + logger.info(" matched filter SNR = {:.2f}".format( + self.meta_data['matched_filter_SNR'])) for key in parameters: logger.info(' {} = {}'.format(key, parameters[key])) @@ -580,7 +621,8 @@ def _window_power_correction(self): using the :code:`BILBY_INCORRECT_PSD_NORMALIZATION` environment variable. """ if string_to_boolean( - os.environ.get("BILBY_INCORRECT_PSD_NORMALIZATION", "FALSE").upper() + os.environ.get("BILBY_INCORRECT_PSD_NORMALIZATION", + "FALSE").upper() ): return self.strain_data.window_factor else: @@ -725,11 +767,12 @@ def matched_filter_snr(self, signal): """ return gwutils.matched_filter_snr( signal=signal[self.strain_data.frequency_mask], - frequency_domain_strain=self.strain_data.frequency_domain_strain[self.strain_data.frequency_mask], + frequency_domain_strain=self.strain_data.frequency_domain_strain[ + self.strain_data.frequency_mask], power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask], duration=self.strain_data.duration) - def whiten_frequency_series(self, frequency_series : np.array) -> np.array: + def whiten_frequency_series(self, frequency_series: np.array) -> np.array: """Whitens a frequency series with the noise properties of the detector .. math:: @@ -752,7 +795,7 @@ def whiten_frequency_series(self, frequency_series : np.array) -> np.array: def get_whitened_time_series_from_whitened_frequency_series( self, - whitened_frequency_series : np.array + whitened_frequency_series: np.array ) -> np.array: """Gets the whitened time series from a whitened frequency series. @@ -837,10 +880,12 @@ def save_data(self, outdir, label=None): if label is None: filename_asd = '{}/{}_asd.dat'.format(outdir, self.name) - filename_data = '{}/{}_frequency_domain_data.dat'.format(outdir, self.name) + filename_data = '{}/{}_frequency_domain_data.dat'.format( + outdir, self.name) else: filename_asd = '{}/{}_{}_asd.dat'.format(outdir, self.name, label) - filename_data = '{}/{}_{}_frequency_domain_data.dat'.format(outdir, self.name, label) + filename_data = '{}/{}_{}_frequency_domain_data.dat'.format( + outdir, self.name, label) np.savetxt(filename_data, np.array( [self.strain_data.frequency_array, @@ -859,7 +904,8 @@ def plot_data(self, signal=None, outdir='.', label=None): return fig, ax = plt.subplots() - df = self.strain_data.frequency_array[1] - self.strain_data.frequency_array[0] + df = self.strain_data.frequency_array[1] - \ + self.strain_data.frequency_array[0] asd = gwutils.asd_from_freq_series( freq_data=self.strain_data.frequency_domain_strain, df=df) @@ -992,7 +1038,8 @@ def _filename_from_outdir_label_extension(outdir, label, extension="h5"): )) def to_pickle(self, outdir="outdir", label=None): utils.check_directory_exists_and_if_not_mkdir('outdir') - filename = self._filename_from_outdir_label_extension(outdir, label, extension="pkl") + filename = self._filename_from_outdir_label_extension( + outdir, label, extension="pkl") safe_file_dump(self, filename, "dill") @classmethod diff --git a/examples/gw_examples/injection_examples/add_glitch_into_ifo.py b/examples/gw_examples/injection_examples/add_glitch_into_ifo.py index 4dbbc90e6..44c3d1285 100644 --- a/examples/gw_examples/injection_examples/add_glitch_into_ifo.py +++ b/examples/gw_examples/injection_examples/add_glitch_into_ifo.py @@ -3,8 +3,8 @@ import matplotlib.pyplot as plt import numpy as np -# To be installed from https://git.ligo.org/melissa.lopez/gengli -# or use glitch generator of your choice. +# gengli to be installed from https://git.ligo.org/melissa.lopez/gengli +# or use glitch generator of your choice or use an analytic glitch time series duration = 4.0 sampling_frequency = 1024.0 From f8a3dbc9464849bad490ce3dfa2817b29ba64aaf Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Thu, 5 Mar 2026 10:02:09 +0100 Subject: [PATCH 08/10] internally call inject_signal function for glitch injection task --- bilby/gw/detector/interferometer.py | 69 ++++++++++++++++------------- 1 file changed, 38 insertions(+), 31 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index ec2130602..e77e7927c 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -401,42 +401,42 @@ def check_signal_duration(self, parameters, raise_error=True): logger.warning(msg) def adjust_optimal_snr(self, frequency_domain_strain, target_snr): - """ Scale the time-domain strain data of the glitch so that it's optimal SNR is equal to the the glitch_snr - Roll the glitch in such a way that the maxima coincides with the glitch_injection_time. + """ Rescale a frequency-domain strain array to a target optimal SNR. Parameters ========== - glitch_strain: numpy array - Time-domain strain data of the glitch to scale. - glitch_snr: float - Desired optimal SNR of the glitch. - glitch_injection_time: float - GPS time of the occurrence of the glitch + frequency_domain_strain: numpy.array + Frequency-domain strain to rescale. + target_snr: float + Desired optimal SNR. Returns ======= - scaled_glitch: numpy array - frequency-domain strain of the glitch rescaled so that its optimal SNR matches - `target_glitch_snr` and rolled in such a way thatthe maxima coincides with the glitch_injection_time. + numpy.ndarray + Frequency-domain strain rescaled so that its optimal SNR equals ``target_snr``. """ temporary_snr_squared = np.real( self.optimal_snr_squared(signal=frequency_domain_strain)) return frequency_domain_strain * target_snr / np.sqrt(temporary_snr_squared) - def inject_glitch(self, glitch_snr, glitch_parameters=None, glitch_time_domain_strain=None, - glitch_sample_times=None, glitch_injection_time=None, glitch_waveform_generator=None): + def inject_glitch(self, glitch_parameters=None, glitch_time_domain_strain=None, + glitch_sample_times=None, glitch_waveform_generator=None): """ Inject a glitch into the interferometer data. Parameters ========== - glitch: numpy array + glitch_parameters: dict, optional + Dictionary of glitch parameters. + Must contain ``onset_time`` (the GPS time at which the glitch peak should occur). + Must contain ``snr``. The glitch will be rescaled in such a + way that it's optimal SNR matches to this value. + glitch_time_domain_strain: numpy.array, optional Time-domain strain data of the glitch to inject. - glitch_injection_time: float - The maxima of the time domain strain of the glitch will conicide with the glitch_injection_time. - glitch_sampling_frequency: float - Sampling frequency of the input glitch array. - glitch_snr: float - The glitch will be rescaled in such a way that it's optimal SNR matches with the glitch_snr. + glitch_sample_times: numpy.array + Array of sample times corresponding to ``glitch_time_domain_strain``. + glitch_waveform_generator: bilby.gw.WaveformGenerator, + A waveform generator used to produce the glitch frequency-domain + strain. """ if glitch_parameters is None: @@ -449,6 +449,7 @@ def inject_glitch(self, glitch_snr, glitch_parameters=None, glitch_time_domain_s "inject_glitch needs one of glitch_waveform_generator or " "glitch_time_domain_strain.") elif glitch_time_domain_strain is not None: + # Populate ra, dec, psi. Necessary to make an injection. glitch_sampling_frequency = 1.0 / (glitch_sample_times[1] - glitch_sample_times[0]) # Resample the glitch to match with the interferometer's sampling frequency if glitch_sampling_frequency != self.sampling_frequency: @@ -463,20 +464,20 @@ def inject_glitch(self, glitch_snr, glitch_parameters=None, glitch_time_domain_s padded_glitch[:len(glitch_time_domain_strain) ] += glitch_time_domain_strain - # Since the maxima of the glitch is not at 0. - glitch_parameters['geocent_time'] = glitch_injection_time - \ - (glitch_sample_times[np.argmax(padded_glitch)]) - # Convert to frequency domain glitch glitch_frequency_domain_strain, _ = utils.nfft( padded_glitch, self.sampling_frequency) # Adjust SNR glitch_frequency_domain_strain = self.adjust_optimal_snr( - glitch_frequency_domain_strain, glitch_snr) + glitch_frequency_domain_strain, glitch_parameters['snr']) injection_polarizations = { self.name: glitch_frequency_domain_strain} + # Roll the glitch since the glitch maxima of the glitch may not align in the same way as the signal maxima. + glitch_parameters['geocent_time'] = glitch_parameters["onset_time"] \ + - glitch_sample_times[np.argmax(glitch_time_domain_strain)] + # Usual inject method self.inject_signal_from_waveform_polarizations(parameters=glitch_parameters, injection_polarizations=injection_polarizations) @@ -485,16 +486,22 @@ def inject_glitch(self, glitch_snr, glitch_parameters=None, glitch_time_domain_s glitch_frequency_domain_strain = glitch_waveform_generator.frequency_domain_strain( glitch_parameters) glitch_frequency_domain_strain = self.adjust_optimal_snr( - glitch_frequency_domain_strain, glitch_snr) - injection_polarizations = { - self.name: glitch_frequency_domain_strain} + np.asarray(glitch_frequency_domain_strain[0]), glitch_parameters['snr']) + injection_polarizations = {self.name: glitch_frequency_domain_strain} + + # Populate ra, dec, psi. Necessary to make an injection. + for key, val in [('ra', 0.0), ('dec', 0.0), ('psi', 0.0)]: + glitch_parameters.setdefault(key, val) + # Roll the glitch since the glitch maxima of the glitch may not align in the same way as the signal maxima. + glitch_time_domain_strain = glitch_waveform_generator.time_domain_strain(glitch_parameters) + glitch_parameters['geocent_time'] = glitch_parameters["onset_time"] - \ + glitch_waveform_generator.time_array[np.argmax(glitch_time_domain_strain)] + self.inject_signal_from_waveform_polarizations(parameters=glitch_parameters, injection_polarizations=injection_polarizations) logger.info("Injected a glitch in: {}".format(self.name)) - logger.info("Optimal SNR of the glitch: {}".format(glitch_snr)) - logger.info("Matched filter SNR of the glitch: {}".format( - self.matched_filter_snr(signal=glitch_frequency_domain_strain))) + logger.info("Optimal SNR of the glitch: {}".format(glitch_parameters['snr'])) return glitch_time_domain_strain From ad1de7704e5fc538477296795f4da69f0ecccb7d Mon Sep 17 00:00:00 2001 From: Narola Harsh <51940885+narolaharsh@users.noreply.github.com> Date: Thu, 5 Mar 2026 10:49:20 +0100 Subject: [PATCH 09/10] Update example --- .../gw_examples/injection_examples/glitch.py | 139 ++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100644 examples/gw_examples/injection_examples/glitch.py diff --git a/examples/gw_examples/injection_examples/glitch.py b/examples/gw_examples/injection_examples/glitch.py new file mode 100644 index 000000000..fdc89f58a --- /dev/null +++ b/examples/gw_examples/injection_examples/glitch.py @@ -0,0 +1,139 @@ +import bilby +import matplotlib.pyplot as plt +import numpy as np + +duration = 8.0 +sampling_frequency = 2048.0 +glitch_snr = 12 +injection_parameters = dict( + mass_1=36.0, + mass_2=29.0, + a_1=0.4, + a_2=0.3, + tilt_1=0.5, + tilt_2=1.0, + phi_12=1.7, + phi_jl=0.3, + luminosity_distance=2100.0, + theta_jn=0.4, + psi=2.659, + phase=1.3, + geocent_time=205.0, + ra=1.375, + dec=-1.2108, +) + +bilby.core.utils.random.seed(23) + + +def sine_gaussian_glitch( + time_array, amplitude=1e-22, f0=100.0, tau=0.1, t0=0.0, phase=0.0 +): + """Generate a sine-Gaussian waveform. + + Parameters + ---------- + time_array : array_like + Time samples. + amplitude : float + Peak amplitude. + f0 : float + Central frequency in Hz. + tau : float + Decay time (Gaussian width) in seconds. + t0 : float + Central time of the burst in seconds. + phase : float + Phase offset in radians. + """ + envelope = np.exp(-((time_array - t0) ** 2) / (2 * tau**2)) + return amplitude * envelope * np.sin(2 * np.pi * f0 * (time_array - t0) + phase) + + +# Fixed arguments passed into the source model +waveform_arguments = dict( + waveform_approximant="IMRPhenomXPHM", + reference_frequency=50.0, + minimum_frequency=20.0, +) + +# Create the waveform_generator using a LAL BinaryBlackHole source function +# the generator will convert all the parameters +waveform_generator = bilby.gw.WaveformGenerator( + duration=duration, + sampling_frequency=sampling_frequency, + frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, + parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments=waveform_arguments, +) + +# Set up interferometers. In this case we'll use two interferometers +# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design +# sensitivity +ifos = bilby.gw.detector.InterferometerList(["H1", "L1"]) +ifos.set_strain_data_from_zero_noise( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=injection_parameters["geocent_time"] - 4, +) + +ifos.inject_signal( + waveform_generator=waveform_generator, parameters=injection_parameters +) + +inject_glitch_from_time_domain_strain = True + +if inject_glitch_from_time_domain_strain: + glitch_sample_times = np.arange(0, 2, 1 / 256) + from scipy.signal.windows import tukey + + glitch = sine_gaussian_glitch(glitch_sample_times, tau=0.1, t0=0.5) * tukey( + len(glitch_sample_times), alpha=0.2 + ) + + # Inject glitch + glitch_parameters = { + "onset_time": injection_parameters["geocent_time"] + 1, + "snr": glitch_snr, + } + ifos[0].inject_glitch( + glitch_snr, + glitch_parameters=glitch_parameters, + glitch_time_domain_strain=glitch, + glitch_sample_times=glitch_sample_times, + ) + + +inject_glitch_from_a_glitch_waveform_generator = False +if inject_glitch_from_a_glitch_waveform_generator: + + glitch_waveform_generator = bilby.gw.WaveformGenerator( + duration=duration, + sampling_frequency=sampling_frequency, + time_domain_source_model=sine_gaussian_glitch, + ) + glitch_parameters = { + "onset_time": injection_parameters["geocent_time"] + 1, + "amplitude": 1e-22, + "f0": 100.0, + "tau": 0.1, + "t0": 1.3, + "phase": 0.0, + "snr": glitch_snr, + } + + ifos[0].inject_glitch( + glitch_parameters=glitch_parameters, + glitch_waveform_generator=glitch_waveform_generator, + ) + + +fig, axes = plt.subplots(2, 1) +ax = axes[0] + +ax.axvline(x=injection_parameters["geocent_time"], color="black", label="Merger time") +ax.axvline(x=glitch_parameters["onset_time"], color="blue", label="glitch onset time") +ax.legend() + +ax.plot(ifos[0].time_array, ifos[0].time_domain_strain) +fig.savefig("./time_domain_strain.pdf") From 11f85b2e65c05885055d82e238b70e4cf5943c85 Mon Sep 17 00:00:00 2001 From: Harsh Narola Date: Thu, 5 Mar 2026 11:30:21 +0100 Subject: [PATCH 10/10] Delete examples/gw_examples/injection_examples/add_glitch_into_ifo.py --- .../injection_examples/add_glitch_into_ifo.py | 80 ------------------- 1 file changed, 80 deletions(-) delete mode 100644 examples/gw_examples/injection_examples/add_glitch_into_ifo.py diff --git a/examples/gw_examples/injection_examples/add_glitch_into_ifo.py b/examples/gw_examples/injection_examples/add_glitch_into_ifo.py deleted file mode 100644 index 44c3d1285..000000000 --- a/examples/gw_examples/injection_examples/add_glitch_into_ifo.py +++ /dev/null @@ -1,80 +0,0 @@ -import bilby -import gengli -import matplotlib.pyplot as plt -import numpy as np - -# gengli to be installed from https://git.ligo.org/melissa.lopez/gengli -# or use glitch generator of your choice or use an analytic glitch time series - -duration = 4.0 -sampling_frequency = 1024.0 -glitch_snr = 12 -injection_parameters = dict( - mass_1=36.0, - mass_2=29.0, - a_1=0.4, - a_2=0.3, - tilt_1=0.5, - tilt_2=1.0, - phi_12=1.7, - phi_jl=0.3, - luminosity_distance=2000.0, - theta_jn=0.4, - psi=2.659, - phase=1.3, - geocent_time=1126259642.413, - ra=1.375, - dec=-1.2108, -) - -bilby.core.utils.random.seed(23) - -# Fixed arguments passed into the source model -waveform_arguments = dict( - waveform_approximant="IMRPhenomXPHM", - reference_frequency=50.0, - minimum_frequency=20.0, -) - -# Create the waveform_generator using a LAL BinaryBlackHole source function -# the generator will convert all the parameters -waveform_generator = bilby.gw.WaveformGenerator( - duration=duration, - sampling_frequency=sampling_frequency, - frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, - parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters, - waveform_arguments=waveform_arguments, -) - -# Set up interferometers. In this case we'll use two interferometers -# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design -# sensitivity -ifos = bilby.gw.detector.InterferometerList(["H1", "L1"]) -ifos.set_strain_data_from_power_spectral_densities( - sampling_frequency=sampling_frequency, - duration=duration, - start_time=injection_parameters["geocent_time"] - 2, -) - -ifos.inject_signal( - waveform_generator=waveform_generator, parameters=injection_parameters -) - -generator = gengli.glitch_generator("L1") -glitch = generator.get_glitch(seed=200) * 1e-23 - -# Inject glitch -ifos[0].inject_glitch( - glitch, - glitch_injection_time=ifos[0].start_time + 2, - glitch_sampling_frequency=sampling_frequency, - glitch_snr=glitch_snr, -) - -fig, axes = plt.subplots(2, 1) -ax = axes[0] -ax.plot(ifos[0].time_array, ifos[0].time_domain_strain) - -ax = axes[1] -ax.plot(np.arange(len(glitch)) / sampling_frequency, glitch) -fig.savefig("./time_domain_strain.pdf")