From 7d26f03362a408bb230752d8a23e7eb16732fc0c Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Sat, 1 Dec 2018 16:37:05 -0800
Subject: [PATCH] gstlal-calibration: Numerous changes to make output cleaned
h(t) independent of start time, so that jobs run in parallel will produce
smooth output without kinks. lal_tdwhiten: Added ability to synchronize data
timestamps with filter updates.
---
gstlal-calibration/bin/gstlal_compute_strain | 123 +++---
.../tests/H1DCS_AllCorrections_Cleaning.ini | 28 +-
.../H1DCS_AllCorrections_Cleaning_TEST.ini | 28 +-
.../H1DCS_TEST_AllCorrections_Cleaning.ini | 363 ------------------
...GDS_LowLatency_AllCorrections_Cleaning.ini | 2 -
.../gst/lal/gstlal_adaptivefirfilt.c | 76 ++++
.../gst/lal/gstlal_adaptivefirfilt.h | 2 +
.../gst/lal/gstlal_transferfunction.c | 84 +++-
.../gst/lal/gstlal_transferfunction.h | 2 +
.../python/calibration_parts.py | 39 +-
.../tests/check_calibration/Makefile | 15 +-
gstlal-calibration/tests/lal_gate_test.py | 9 +-
gstlal-ugly/gst/lal/gstlal_tdwhiten.c | 178 ++++++++-
gstlal-ugly/gst/lal/gstlal_tdwhiten.h | 4 +
14 files changed, 455 insertions(+), 498 deletions(-)
delete mode 100644 gstlal-calibration/config_files/O2/H1/tests/H1DCS_TEST_AllCorrections_Cleaning.ini
diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain
index ba533a947b..b60c5c3bad 100755
--- a/gstlal-calibration/bin/gstlal_compute_strain
+++ b/gstlal-calibration/bin/gstlal_compute_strain
@@ -243,6 +243,7 @@ compute_calib_factors_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(
compute_calib_factors_complex_caps = "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % compute_factors_sr
# Short cut names for a few other items that appear numerous times
+input_frame_duration = int(InputConfigs["inputframeduration"])
filter_latency_factor = float(PipelineConfigs["filterlatency"])
demodulation_filter_time = int(TDCFConfigs["demodulationfiltertime"])
coherence_unc_threshold = float(TDCFConfigs["coherenceuncthreshold"])
@@ -821,7 +822,7 @@ if not verbose:
#
if InputConfigs["datasource"] == "lvshm": # Data is to be read from shared memory; "low-latency" mode
- src = pipeparts.mklvshmsrc(pipeline, shm_name = InputConfigs["shmpartition"], assumed_duration = int(InputConfigs["inputframeduration"]))
+ src = pipeparts.mklvshmsrc(pipeline, shm_name = InputConfigs["shmpartition"], assumed_duration = input_frame_duration)
elif InputConfigs["datasource"] == "frames": # Data is to be read from frame files; "offline" mode
src = pipeparts.mklalcachesrc(pipeline, location = options.frame_cache, cache_dsc_regex = instrument)
@@ -1372,19 +1373,19 @@ if compute_fs or compute_srcq:
if apply_complex_kappatst:
# We will apply an adaptive FIR filter to the TST component of the actuation that includes time-dependence in the gain and computational time delay
- adaptive_tst_filter = pipeparts.mkgeneric(pipeline, smooth_ktsttee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = esd_act_line_freq, static_filter = tstfilt, variable_filter_length = len(tstfilt), adaptive_filter_length = len(tstfilt), tukey_param = 0.5, filter_sample_rate = tstchainsr)
+ adaptive_tst_filter = pipeparts.mkgeneric(pipeline, smooth_ktsttee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = esd_act_line_freq, static_filter = tstfilt, variable_filter_length = len(tstfilt), adaptive_filter_length = len(tstfilt), tukey_param = 0.5, filter_sample_rate = tstchainsr, filter_timeshift = 1000000000 * actuation_filter_update_time, name = "adaptive_tst_filter")
if apply_complex_kappapum:
# We will apply an adaptive FIR filter to the PUM component of the actuation that includes time-dependence in the gain and computational time delay
- adaptive_pum_filter = pipeparts.mkgeneric(pipeline, smooth_kpumtee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = esd_act_line_freq, static_filter = pumfilt, variable_filter_length = len(pumfilt), adaptive_filter_length = len(pumfilt), tukey_param = 0.5, filter_sample_rate = pumchainsr)
+ adaptive_pum_filter = pipeparts.mkgeneric(pipeline, smooth_kpumtee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = esd_act_line_freq, static_filter = pumfilt, variable_filter_length = len(pumfilt), adaptive_filter_length = len(pumfilt), tukey_param = 0.5, filter_sample_rate = pumchainsr, filter_timeshift = 1000000000 * actuation_filter_update_time, name = "adaptive_pum_filter")
if apply_complex_kappauim:
# We will apply an adaptive FIR filter to the UIM component of the actuation that includes time-dependence in the gain and computational time delay
- adaptive_uim_filter = pipeparts.mkgeneric(pipeline, smooth_kuimtee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = esd_act_line_freq, static_filter = uimfilt, variable_filter_length = len(uimfilt), adaptive_filter_length = len(uimfilt), tukey_param = 0.5, filter_sample_rate = uimchainsr)
+ adaptive_uim_filter = pipeparts.mkgeneric(pipeline, smooth_kuimtee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = esd_act_line_freq, static_filter = uimfilt, variable_filter_length = len(uimfilt), adaptive_filter_length = len(uimfilt), tukey_param = 0.5, filter_sample_rate = uimchainsr, filter_timeshift = 1000000000 * actuation_filter_update_time, name = "adaptive_uim_filter")
if apply_complex_kappapu:
# We will apply an adaptive FIR filter to the PUM/UIM component of the actuation that includes time-dependence in the gain and computational time delay
- adaptive_pumuim_filter = pipeparts.mkgeneric(pipeline, smooth_kputee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = darm_ctrl_line_freq, static_filter = pumuimfilt, variable_filter_length = len(pumuimfilt), adaptive_filter_length = len(pumuimfilt), tukey_param = 0.5, filter_sample_rate = pumuimchainsr)
+ adaptive_pumuim_filter = pipeparts.mkgeneric(pipeline, smooth_kputee, "lal_adaptivefirfilt", update_samples = int(actuation_filter_update_time * compute_factors_sr), average_samples = int(actuation_filter_averaging_time * compute_factors_sr), phase_measurement_frequency = darm_ctrl_line_freq, static_filter = pumuimfilt, variable_filter_length = len(pumuimfilt), adaptive_filter_length = len(pumuimfilt), tukey_param = 0.5, filter_sample_rate = pumuimchainsr, filter_timeshift = 1000000000 * actuation_filter_update_time, name = "adaptive_pumuim_filter")
if apply_fcc or apply_fs or apply_srcq:
# We will apply an adaptive FIR filter to DARM_ERR that allows corrections for poles, zeros, and gain
@@ -1478,7 +1479,7 @@ if apply_fcc or apply_fs or apply_srcq:
# Now interleave the correction channels in tdep_zpk and feed them into lal_adaptivefirfilt to update the inverse sensing filter
tdep_zpk = calibration_parts.mkinterleave(pipeline, tdep_zpk, complex_data = True)
- adaptive_invsens_filter = pipeparts.mkgeneric(pipeline, tdep_zpk, "lal_adaptivefirfilt", update_samples = int(sensing_filter_update_time * compute_factors_sr), average_samples = int(sensing_filter_averaging_time * compute_factors_sr), num_zeros = variable_invsens_zeros, num_poles = 0, static_poles = static_invsens_poles, static_filter = reschainfilt, variable_filter_length = len(reschainfilt), adaptive_filter_length = len(reschainfilt), tukey_param = 0.5, filter_sample_rate = hoft_sr)
+ adaptive_invsens_filter = pipeparts.mkgeneric(pipeline, tdep_zpk, "lal_adaptivefirfilt", update_samples = int(sensing_filter_update_time * compute_factors_sr), average_samples = int(sensing_filter_averaging_time * compute_factors_sr), num_zeros = variable_invsens_zeros, num_poles = 0, static_poles = static_invsens_poles, static_filter = reschainfilt, variable_filter_length = len(reschainfilt), adaptive_filter_length = len(reschainfilt), tukey_param = 0.5, filter_sample_rate = hoft_sr, filter_timeshift = 1000000000 * sensing_filter_update_time, name = "adaptive_invsens_filter")
#
# CONTROL BRANCH
@@ -1573,13 +1574,17 @@ if any(act_highpass):
tst_filter_latency += float(act_highpass_delay)/tstchainsr
if apply_complex_kappatst:
- if filter_latency_factor > 0:
- # Make the TST path wait until a correction filter for kappa_TST is ready
- tst = calibration_parts.mkqueue(pipeline, tst, min_length = filter_latency_factor * (demodulation_filter_time + (median_smoothing_samples + factors_average_samples) / compute_factors_sr))
- # Filter the TST chain with an adaptive TST actuation filter that includes a linear-phase correction from kappa_tst
- tst = pipeparts.mkgeneric(pipeline, tst, "lal_tdwhiten", kernel = tstfilt[::-1], latency = tstdelay, taper_length = actuation_filter_taper_length)
- # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
- adaptive_tst_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, tst, "adaptive-filter", "kernel")
+ if filter_latency_factor == 0:
+ # Filter the TST chain with an adaptive TST actuation filter that includes a gain and linear-phase correction from kappa_tst
+ tst = pipeparts.mkgeneric(pipeline, tst, "lal_tdwhiten", kernel = tstfilt[::-1], latency = tstdelay, taper_length = actuation_filter_taper_length, name = "TST_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_tst_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, tst, "adaptive_filter", "kernel")
+ else:
+ # Filter the TST chain with an adaptive TST actuation filter that includes a gain and linear-phase correction from kappa_tst
+ tst = pipeparts.mkgeneric(pipeline, calibration_parts.mkqueue(pipeline, tst), "lal_tdwhiten", kernel = tstfilt[::-1], latency = tstdelay, taper_length = actuation_filter_taper_length, kernel_endtime = 0, name = "TST_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_tst_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, tst, "adaptive_filter", "kernel")
+ adaptive_tst_filter.connect("notify::filter-endtime", calibration_parts.update_property_simple, tst, "filter_endtime", "kernel_endtime")
else:
# Filter the TST chain with the static TST actuation filter
@@ -1617,25 +1622,33 @@ if apply_kappapum or apply_kappauim or apply_complex_kappapum or apply_complex_k
uim_filter_latency += float(act_highpass_delay)/uimchainsr
if apply_complex_kappapum:
- if filter_latency_factor > 0:
- # Make the PUM path wait until a correction filter for kappa_PUM is ready
- pum = calibration_parts.mkqueue(pipeline, pum, min_length = filter_latency_factor * (demodulation_filter_time + (median_smoothing_samples + factors_average_samples) / compute_factors_sr))
- # Filter the PUM chain with an adaptive PUM actuation filter that includes a linear-phase correction from kappa_pum
- pum = pipeparts.mkgeneric(pipeline, pum, "lal_tdwhiten", kernel = pumfilt[::-1], latency = pumdelay, taper_length = actuation_filter_taper_length)
- # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
- adaptive_pum_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, pum, "adaptive-filter", "kernel")
+ if filter_latency_factor == 0:
+ # Filter the PUM chain with an adaptive PUM actuation filter that includes a gain and linear-phase correction from kappa_pum
+ pum = pipeparts.mkgeneric(pipeline, pum, "lal_tdwhiten", kernel = pumfilt[::-1], latency = pumdelay, taper_length = actuation_filter_taper_length, name = "PUM_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_pum_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, pum, "adaptive_filter", "kernel")
+ else:
+ # Filter the PUM chain with an adaptive PUM actuation filter that includes a gain and linear-phase correction from kappa_pum
+ pum = pipeparts.mkgeneric(pipeline, calibration_parts.mkqueue(pipeline, pum), "lal_tdwhiten", kernel = pumfilt[::-1], latency = pumdelay, taper_length = actuation_filter_taper_length, kernel_endtime = 0, name = "PUM_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_pum_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, pum, "adaptive_filter", "kernel")
+ adaptive_pum_filter.connect("notify::filter-endtime", calibration_parts.update_property_simple, pum, "filter_endtime", "kernel_endtime")
else:
# Filter the PUM chain with the static PUM actuation filter
pum = pipeparts.mkfirbank(pipeline, pum, latency = pumdelay, fir_matrix = [pumfilt[::-1]], time_domain = td)
if apply_complex_kappauim:
- if filter_latency_factor > 0:
- # Make the UIM path wait until a correction filter for kappa_UIM is ready
- uim = calibration_parts.mkqueue(pipeline, uim, min_length = filter_latency_factor * (demodulation_filter_time + (median_smoothing_samples + factors_average_samples) / compute_factors_sr))
- # Filter the UIM chain with an adaptive UIM actuation filter that includes a linear-phase correction from kappa_uim
- uim = pipeparts.mkgeneric(pipeline, uim, "lal_tdwhiten", kernel = uimfilt[::-1], latency = uimdelay, taper_length = actuation_filter_taper_length)
- # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
- adaptive_uim_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, uim, "adaptive-filter", "kernel")
+ if filter_latency_factor == 0:
+ # Filter the UIM chain with an adaptive UIM actuation filter that includes a gain and linear-phase correction from kappa_uim
+ uim = pipeparts.mkgeneric(pipeline, uim, "lal_tdwhiten", kernel = uimfilt[::-1], latency = uimdelay, taper_length = actuation_filter_taper_length, name = "UIM_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_uim_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, uim, "adaptive_filter", "kernel")
+ else:
+ # Filter the UIM chain with an adaptive UIM actuation filter that includes a gain and linear-phase correction from kappa_uim
+ uim = pipeparts.mkgeneric(pipeline, calibration_parts.mkqueue(pipeline, uim), "lal_tdwhiten", kernel = uimfilt[::-1], latency = uimdelay, taper_length = actuation_filter_taper_length, kernel_endtime = 0, name = "UIM_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_uim_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, uim, "adaptive_filter", "kernel")
+ adaptive_uim_filter.connect("notify::filter-endtime", calibration_parts.update_property_simple, uim, "filter_endtime", "kernel_endtime")
else:
# Filter the UIM chain with the static UIM actuation filter
uim = pipeparts.mkfirbank(pipeline, uim, latency = uimdelay, fir_matrix = [uimfilt[::-1]], time_domain = td)
@@ -1679,13 +1692,17 @@ else:
pumuim_filter_latency += float(act_highpass_delay)/pumuimchainsr
if apply_complex_kappapu:
- if filter_latency_factor > 0:
- # Make the PUM/UIM path wait until a correction filter for kappa_PU is ready
- pumuim = calibration_parts.mkqueue(pipeline, pumuim, min_length = filter_latency_factor * (demodulation_filter_time + (median_smoothing_samples + factors_average_samples) / compute_factors_sr))
- # Filter the PUM/UIM chain with an adaptive PUM/UIM actuation filter that includes a linear-phase correction from kappa_pu
- pumuim = pipeparts.mkgeneric(pipeline, pumuim, "lal_tdwhiten", kernel = pumuimfilt[::-1], latency = pumuimdelay, taper_length = actuation_filter_taper_length)
- # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
- adaptive_pumuim_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, pumuim, "adaptive-filter", "kernel")
+ if filter_latency_factor == 0:
+ # Filter the PUM/UIM chain with an adaptive PUM/UIM actuation filter that includes a gain and linear-phase correction from kappa_pu
+ pumuim = pipeparts.mkgeneric(pipeline, pumuim, "lal_tdwhiten", kernel = pumuimfilt[::-1], latency = pumuimdelay, taper_length = actuation_filter_taper_length, name = "PUMUIM_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_pumuim_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, pumuim, "adaptive_filter", "kernel")
+ else:
+ # Filter the PUM/UIM chain with an adaptive PUM/UIM actuation filter that includes a gain and linear-phase correction from kappa_pu
+ pumuim = pipeparts.mkgeneric(pipeline, calibration_parts.mkqueue(pipeline, pumuim), "lal_tdwhiten", kernel = pumuimfilt[::-1], latency = pumuimdelay, taper_length = actuation_filter_taper_length, kernel_endtime = 0, name = "PUMUIM_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_pumuim_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, pumuim, "adaptive_filter", "kernel")
+ adaptive_pumuim_filter.connect("notify::filter-endtime", calibration_parts.update_property_simple, pumuim, "filter_endtime", "kernel_endtime")
else:
# Filter the PUM/UIM chain with the static PUM/UIM actuation filter
@@ -1754,13 +1771,18 @@ if any(invsens_highpass):
res_filter_latency += float(invsens_highpass_delay) / invsens_highpass_sr
if apply_fcc or apply_fs or apply_srcq:
- if filter_latency_factor > 0:
- # Make the residual path wait until a correction filter is ready
- res = calibration_parts.mkqueue(pipeline, res, min_length = filter_latency_factor * (demodulation_filter_time + (median_smoothing_samples + factors_average_samples) / compute_factors_sr))
- # Apply an adaptive filter to include the time-dependence of any sensing function parameters
- res = pipeparts.mkgeneric(pipeline, res, "lal_tdwhiten", kernel = reschainfilt[::-1], latency = reschaindelay, taper_length = sensing_filter_taper_length)
- # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
- adaptive_invsens_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, res, "adaptive-filter", "kernel")
+ if filter_latency_factor == 0:
+ # Apply an adaptive filter to include the time-dependence of any sensing function parameters
+ res = pipeparts.mkgeneric(pipeline, res, "lal_tdwhiten", kernel = reschainfilt[::-1], latency = reschaindelay, taper_length = sensing_filter_taper_length, name = "RES_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_invsens_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, res, "adaptive_filter", "kernel")
+ else:
+ # Apply an adaptive filter to include the time-dependence of any sensing function parameters
+ res = pipeparts.mkgeneric(pipeline, calibration_parts.mkqueue(pipeline, res), "lal_tdwhiten", kernel = reschainfilt[::-1], latency = reschaindelay, taper_length = sensing_filter_taper_length, kernel_endtime = 0, name = "RES_filter")
+ # Hook up the adaptive filter from lal_adaptivefirfilt to lal_tdwhiten so that the filter gets updated
+ adaptive_invsens_filter.connect("notify::adaptive-filter", calibration_parts.update_filter, res, "adaptive_filter", "kernel")
+ adaptive_invsens_filter.connect("notify::filter-endtime", calibration_parts.update_property_simple, res, "filter_endtime", "kernel_endtime")
+
# # Correct for time-dependence of f_cc
#if apply_fcc:
#
@@ -2217,12 +2239,8 @@ if remove_power_lines:
if witness_channel_list is not None:
# Remove initial data from computation of transfer functions and wait until the filters and kappas settle
witness_delay_time = witness_tf_time_shift if witness_tf_parallel_mode else (filter_settle_time + (1.0 - filter_latency_factor) * (demodulation_filter_time + median_smoothing_samples / compute_factors_sr + factors_average_samples / compute_factors_sr))
- # In high latency, make the witnesses wait to be filtered until new filters are ready
- witness_wait_time = (filter_settle_time + demodulation_filter_time + median_smoothing_samples / compute_factors_sr + factors_average_samples / compute_factors_sr + witness_channel_fft_time / 2.0 * (num_witness_ffts + 1.0)) if filter_latency_factor else 0.0
# How much does the "delay_time" need to increase per iteration of cleaning?
- witness_delay_increment = witness_filter_taper_time + (witness_channel_fft_time / 2.0 * (num_witness_ffts + 1.0) if not filter_latency_factor else 0.0)
- # How much does the "wait_time" need to increase per iteration of cleaning?
- witness_wait_increment = (witness_filter_taper_time + witness_channel_fft_time / 2.0 * (num_witness_ffts + 1.0)) if filter_latency_factor else 0.0
+ witness_delay_increment = witness_filter_taper_time + witness_channel_fft_time / 2.0 * (num_witness_ffts + 1.0) if not filter_latency_factor else 0.0
# If we haven't removed any lines, clean the regular h(t) data
if not (remove_cal_lines or remove_power_lines):
clean_strain = straintee
@@ -2245,9 +2263,8 @@ if witness_channel_list is not None:
witnesses.append(calibration_parts.caps_and_progress(pipeline, head_dict[key], "audio/x-raw, format=F64LE, channels=1, channel-mask=(bitmask)0x0", key))
if len(witnesses) != len(witness_channel_list[i]):
print "WARNING: Not all requested witness channels are being used"
- clean_strain = calibration_parts.clean_data(pipeline, clean_strain, hoft_sr, witnesses, witness_rates[i], witness_fft_samples, witness_fft_overlap, num_witness_ffts, min_witness_ffts, witness_tf_update_samples, witness_fir_samples, witness_frequency_resolution, witness_filter_taper_length, use_median = witness_tf_use_median, parallel_mode = witness_tf_parallel_mode, notch_frequencies = witness_notch_frequencies[i], noisesub_gate_bit = noisesubgatetee, delay_time = witness_delay_time, wait_time = witness_wait_time, critical_lock_loss_time = critical_lock_loss_time, filename = None if witness_tf_filename is None else "%s_%d.txt" % (witness_tf_filename, i))
+ clean_strain = calibration_parts.clean_data(pipeline, clean_strain, hoft_sr, witnesses, witness_rates[i], witness_fft_samples, witness_fft_overlap, num_witness_ffts, min_witness_ffts, witness_tf_update_samples, witness_fir_samples, witness_frequency_resolution, witness_filter_taper_length, use_median = witness_tf_use_median, parallel_mode = witness_tf_parallel_mode, notch_frequencies = witness_notch_frequencies[i], noisesub_gate_bit = noisesubgatetee, delay_time = witness_delay_time, critical_lock_loss_time = critical_lock_loss_time, filename = None if witness_tf_filename is None else "%s_%d.txt" % (witness_tf_filename, i))
witness_delay_time += witness_delay_increment
- witness_wait_time += witness_wait_increment
if remove_cal_lines or remove_power_lines or witness_channel_list is not None:
clean_strain = pipeparts.mkprogressreport(pipeline, clean_strain, "progress_hoft_cleaned_%s" % instrument)
@@ -2267,18 +2284,18 @@ if compute_calib_statevector and (remove_cal_lines or remove_power_lines or witn
mid_rms_rate = pow(2, int(numpy.log(2 * cleaning_check_range_mid_max) / numpy.log(2) + 1.01))
# Compute the RMS of the uncleaned strain in a low-frequency range to test subtraction of actuation lines
- strain_rms_lowfreq = calibration_parts.compute_rms(pipeline, straintee, low_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_low_min, f_max = cleaning_check_range_low_max, filter_latency = filter_latency_factor, rate_out = calibstate_sr, td = td)
+ strain_rms_lowfreq = calibration_parts.compute_rms(pipeline, straintee, low_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_low_min, f_max = cleaning_check_range_low_max, filter_latency = min(filter_latency_factor, 0.5), rate_out = calibstate_sr, td = td)
# Compute the RMS of the cleaned strain in a low-frequency range
- clean_strain_rms_lowfreq = calibration_parts.compute_rms(pipeline, clean_straintee, low_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_low_min, f_max = cleaning_check_range_low_max, filter_latency = filter_latency_factor, rate_out = calibstate_sr, td = td)
+ clean_strain_rms_lowfreq = calibration_parts.compute_rms(pipeline, clean_straintee, low_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_low_min, f_max = cleaning_check_range_low_max, filter_latency = min(filter_latency_factor, 0.5), rate_out = calibstate_sr, td = td)
# Require that ratio RMS(strain) / RMS(clean_strain) > 1.0
clean_hoft_ok_lowfreq = calibration_parts.complex_division(pipeline, strain_rms_lowfreq, clean_strain_rms_lowfreq)
clean_hoft_ok_lowfreq = pipeparts.mkbitvectorgen(pipeline, clean_hoft_ok_lowfreq, bit_vector=pow(2,26), threshold=1.0)
clean_hoft_ok_lowfreq = pipeparts.mkcapsfilter(pipeline, clean_hoft_ok_lowfreq, calibstate_caps)
# Compute the RMS of the uncleaned strain in a mid-frequency range to test subtraction of noise and/or the ~300 Hz pcal line
- strain_rms_midfreq = calibration_parts.compute_rms(pipeline, straintee, mid_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_mid_min, f_max = cleaning_check_range_mid_max, filter_latency = filter_latency_factor, rate_out = calibstate_sr, td = td)
+ strain_rms_midfreq = calibration_parts.compute_rms(pipeline, straintee, mid_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_mid_min, f_max = cleaning_check_range_mid_max, filter_latency = min(filter_latency_factor, 0.5), rate_out = calibstate_sr, td = td)
# Compute the RMS of the cleaned strain in a mid-frequency range
- clean_strain_rms_midfreq = calibration_parts.compute_rms(pipeline, clean_straintee, mid_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_mid_min, f_max = cleaning_check_range_mid_max, filter_latency = filter_latency_factor, rate_out = calibstate_sr, td = td)
+ clean_strain_rms_midfreq = calibration_parts.compute_rms(pipeline, clean_straintee, mid_rms_rate, cleaning_check_rms_time, f_min = cleaning_check_range_mid_min, f_max = cleaning_check_range_mid_max, filter_latency = min(filter_latency_factor, 0.5), rate_out = calibstate_sr, td = td)
# Require that ratio RMS(strain) / RMS(clean_strain) > 1.0
clean_hoft_ok_midfreq = calibration_parts.complex_division(pipeline, strain_rms_midfreq, clean_strain_rms_midfreq)
clean_hoft_ok_midfreq = pipeparts.mkbitvectorgen(pipeline, clean_hoft_ok_midfreq, bit_vector=pow(2,27), threshold=1.0)
@@ -2506,7 +2523,7 @@ elif InputConfigs["datasource"] == "lvshm":
start = int(lal.UTCToGPS(tm))
# start time of first frame file is the desired start time + either filter latency or kappa settling (if computing kappas), whichever is bigger
if compute_kappatst or compute_kappapu or compute_kappac or compute_fcc or compute_fs or compute_srcq:
- output_start = start + max(int(filter_settle_time), demodulation_filter_time + int(TDCFConfigs["mediansmoothingtime"]) + int(TDCFConfigs["tdcfaveragingtime"]))
+ output_start = start + max(int(filter_settle_time), int((1.0 - filter_latency_factor) * (demodulation_filter_time + int(TDCFConfigs["mediansmoothingtime"]) + int(TDCFConfigs["tdcfaveragingtime"]))))
else:
output_start = start + int(filter_settle_time)
diff --git a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini b/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
index 21e36a9e93..f88b99b492 100644
--- a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
+++ b/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
@@ -1,6 +1,4 @@
[InputConfigurations]
-# Track what "version" of config file this is, so that the pipeline knows which options are present in this file
-ConfigVersion: 2
# Filters file containing calibration FIR filters, relative to the directory gstlal-calibration/tests/check_calibration/, from which the pipeline is expected to be run
FiltersFileName: H1DCS_newsrcline_1173225472.npz
# Data source should be set to frames or lvshm
@@ -8,6 +6,8 @@ DataSource: frames
FileChecksum: No
# Right now, SkipBadFiles needs to be off when reading from frames
SkipBadFiles: No
+# Assumed duration of input frames in seconds
+InputFrameDuration: 64
############################################
# If reading from frames use these options #
############################################
@@ -17,8 +17,6 @@ SkipBadFiles: No
# If reading from shared memory use these options #
###################################################
SHMPartition: LHO_Online
-# Assumed duration of input frames in seconds
-InputFrameDuration: 1
[OutputConfigurations]
CompressionScheme: 6
@@ -161,13 +159,13 @@ FactorsFromFiltersFile: Yes
# Updating Sensing and Actuation filters with all frequency-dependent corrections parameters #
##############################################################################################
# Length of time (in seconds) between when inverse-sensing FIR filter is updated
-SensingFilterUpdateTime: 60
+SensingFilterUpdateTime: 64
# Length of time (in seconds) over which the smoothed time-dependent parameters of the sensing function are averaged before updating the filter
SensingFilterAveragingTime: 1
# Number of samples to be used when tapering old inverse sensing filter and ramping in new filter
SensingFilterTaperLength: 32768
# Length of time (in seconds) between when the actuation FIR filters are updated
-ActuationFilterUpdateTime: 60
+ActuationFilterUpdateTime: 64
# Length of time (in seconds) over which the smoothed time-dependent parameters of the actuation function are averaged before updating the filter
ActuationFilterAveragingTime: 1
# Number of samples to be used when tapering old actuation filters and ramping in new filters
@@ -210,7 +208,7 @@ CohUncDARMLine1Channel: CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY
PowerLinesChannel: PEM-EY_MAINSMON_EBAY_1_DQ
# Comma-separated list of witness channels to use to subtract noise from h(t)
# Set to None if no witness channels are to be used
-WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ
+WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ;PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ;ASC-DHARD_P_OUT_DQ,ASC-DHARD_Y_OUT_DQ,ASC-CHARD_P_OUT_DQ,ASC-CHARD_Y_OUT_DQ;LSC-SRCL_IN1_DQ,LSC-MICH_IN1_DQ,LSC-PRCL_IN1_DQ
###############################
# EPICS Records Channel Names #
###############################
@@ -286,7 +284,7 @@ EPICSRefSR: 16
# Sample rate for power lines channel
PowerLinesChannelSR: 1024
# Sample rates at which transfer functions will be computed and witness channels will be filtered, given as a semicolon-separated list, e.g., 2048;2048;512;2048. This must be given if WitnessChannelList is not None, and it must be the same length.
-WitnessChannelSR: 2048
+WitnessChannelSR: 2048;2048;512;2048
# Sample rates at which to compute and record TDCFs
ComputeFactorsSR: 16
RecordFactorsSR: 16
@@ -330,29 +328,29 @@ PowerLinesTFAveragingTime: 128
# The length in seconds of the fast Fourier transforms used to compute transfer functions between witness channels and h(t). The fft's are windowed with Hann windows and overlapped.
WitnessChannelFFTTime: 4.0
# The number of ffts to take before averaging the witness -> h(t) transfer functions calculation. The average is taken after the ratio h(f) / witness(f).
-NumWitnessFFTs: 510
+NumWitnessFFTs: 509
# Sets the minimum number of FFTs necessary to produce the first transfer functions and clean data after data flow starts.
-MinWitnessFFTs: 510
+MinWitnessFFTs: 509
# The length in seconds of the filters applied to the witness channels before subtracting from h(t)
WitnessFIRLength: 0.5
# The frequency resolution of the filters applied to the witness channels before subtracting from h(t). It can be advantageous to lower the frequency resolution in order to average over excess noise.
WitnessFrequencyResolution: 1.0
# List of minima and maxima of frequency ranges where the Fourier transform of h(t) will be replaced by a straight line in the calculation of transfer functions between witness channels and h(t) for noise subtraction. Semicolons separate lists for different sets of witness channels. If no notches are desired, use zeros, e.g., \'0;0;0\'. Here is an example using the expected format: \'495.0,515.0,985.0,1015.0;59,60,119,121;0\' This can be useful, e.g., if there are loud lines in the signal that are not present in the witness channels.
-WitnessNotchFrequencies: 495.0,515.0,985.0,1015.0
+WitnessNotchFrequencies: 12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0
# The amount of time after transfer functions between witness channels and h(t) are finished to begin the calculation of the next set of transfer functions
-WitnessTFUpdateTime: 2
+WitnessTFUpdateTime: 4
# If lock-loss lasts at least this many seconds, transfer functions will revert to those computed at the beginning of a lock stretch. Set to zero to disbale.
CriticalLockLossTime: 0
# The amount of time to use to taper in newly computed FIR filters for witness channels being used for noise subtraction.
-WitnessFilterTaperTime: 10
+WitnessFilterTaperTime: 2
# If writing transfer functions to file, this sets the name. If transfer functions should not be written to file, this should be set to None
WitnessTFFilename: transfer_functions_DCS
# Should the transfer function calculation use a median? If not, an average (mean) is used.
-WitnessTFUseMedian: No
+WitnessTFUseMedian: Yes
# Should transfer functions be computed on a fixed schedule, so that the output does not depend on start time? This is useful for running jobs in parallel. Otherwise, they are computed asap.
WitnessTFParallelMode: Yes
# When using parallel mode, how many seconds later should we shift the time when transfer functions start being computed from a multiple of the cycle period?
-WitnessTFTimeShift: 350
+WitnessTFTimeShift: 360
###############################
# Options for HOFT_CLEAN bits #
###############################
diff --git a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning_TEST.ini b/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning_TEST.ini
index c369ec1364..12f37ab6b1 100644
--- a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning_TEST.ini
+++ b/gstlal-calibration/config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning_TEST.ini
@@ -1,6 +1,4 @@
[InputConfigurations]
-# Track what "version" of config file this is, so that the pipeline knows which options are present in this file
-ConfigVersion: 2
# Filters file containing calibration FIR filters, relative to the directory gstlal-calibration/tests/check_calibration/, from which the pipeline is expected to be run
FiltersFileName: H1DCS_newsrcline_1173225472.npz
# Data source should be set to frames or lvshm
@@ -8,6 +6,8 @@ DataSource: frames
FileChecksum: No
# Right now, SkipBadFiles needs to be off when reading from frames
SkipBadFiles: No
+# Assumed duration of input frames in seconds
+InputFrameDuration: 64
############################################
# If reading from frames use these options #
############################################
@@ -17,8 +17,6 @@ SkipBadFiles: No
# If reading from shared memory use these options #
###################################################
SHMPartition: LHO_Online
-# Assumed duration of input frames in seconds
-InputFrameDuration: 1
[OutputConfigurations]
CompressionScheme: 6
@@ -161,13 +159,13 @@ FactorsFromFiltersFile: Yes
# Updating Sensing and Actuation filters with all frequency-dependent corrections parameters #
##############################################################################################
# Length of time (in seconds) between when inverse-sensing FIR filter is updated
-SensingFilterUpdateTime: 60
+SensingFilterUpdateTime: 64
# Length of time (in seconds) over which the smoothed time-dependent parameters of the sensing function are averaged before updating the filter
SensingFilterAveragingTime: 1
# Number of samples to be used when tapering old inverse sensing filter and ramping in new filter
SensingFilterTaperLength: 32768
# Length of time (in seconds) between when the actuation FIR filters are updated
-ActuationFilterUpdateTime: 60
+ActuationFilterUpdateTime: 64
# Length of time (in seconds) over which the smoothed time-dependent parameters of the actuation function are averaged before updating the filter
ActuationFilterAveragingTime: 1
# Number of samples to be used when tapering old actuation filters and ramping in new filters
@@ -210,7 +208,7 @@ CohUncDARMLine1Channel: CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY
PowerLinesChannel: PEM-EY_MAINSMON_EBAY_1_DQ
# Comma-separated list of witness channels to use to subtract noise from h(t)
# Set to None if no witness channels are to be used
-WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ
+WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ;PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ;ASC-DHARD_P_OUT_DQ,ASC-DHARD_Y_OUT_DQ,ASC-CHARD_P_OUT_DQ,ASC-CHARD_Y_OUT_DQ;LSC-SRCL_IN1_DQ,LSC-MICH_IN1_DQ,LSC-PRCL_IN1_DQ
###############################
# EPICS Records Channel Names #
###############################
@@ -286,7 +284,7 @@ EPICSRefSR: 16
# Sample rate for power lines channel
PowerLinesChannelSR: 1024
# Sample rates at which transfer functions will be computed and witness channels will be filtered, given as a semicolon-separated list, e.g., 2048;2048;512;2048. This must be given if WitnessChannelList is not None, and it must be the same length.
-WitnessChannelSR: 2048
+WitnessChannelSR: 2048;2048;512;2048
# Sample rates at which to compute and record TDCFs
ComputeFactorsSR: 16
RecordFactorsSR: 16
@@ -330,29 +328,29 @@ PowerLinesTFAveragingTime: 128
# The length in seconds of the fast Fourier transforms used to compute transfer functions between witness channels and h(t). The fft's are windowed with Hann windows and overlapped.
WitnessChannelFFTTime: 4.0
# The number of ffts to take before averaging the witness -> h(t) transfer functions calculation. The average is taken after the ratio h(f) / witness(f).
-NumWitnessFFTs: 510
+NumWitnessFFTs: 509
# Sets the minimum number of FFTs necessary to produce the first transfer functions and clean data after data flow starts.
-MinWitnessFFTs: 510
+MinWitnessFFTs: 509
# The length in seconds of the filters applied to the witness channels before subtracting from h(t)
WitnessFIRLength: 0.5
# The frequency resolution of the filters applied to the witness channels before subtracting from h(t). It can be advantageous to lower the frequency resolution in order to average over excess noise.
WitnessFrequencyResolution: 1.0
# List of minima and maxima of frequency ranges where the Fourier transform of h(t) will be replaced by a straight line in the calculation of transfer functions between witness channels and h(t) for noise subtraction. Semicolons separate lists for different sets of witness channels. If no notches are desired, use zeros, e.g., \'0;0;0\'. Here is an example using the expected format: \'495.0,515.0,985.0,1015.0;59,60,119,121;0\' This can be useful, e.g., if there are loud lines in the signal that are not present in the witness channels.
-WitnessNotchFrequencies: 495.0,515.0,985.0,1015.0
+WitnessNotchFrequencies: 12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0
# The amount of time after transfer functions between witness channels and h(t) are finished to begin the calculation of the next set of transfer functions
-WitnessTFUpdateTime: 2
+WitnessTFUpdateTime: 4
# If lock-loss lasts at least this many seconds, transfer functions will revert to those computed at the beginning of a lock stretch. Set to zero to disbale.
CriticalLockLossTime: 0
# The amount of time to use to taper in newly computed FIR filters for witness channels being used for noise subtraction.
-WitnessFilterTaperTime: 10
+WitnessFilterTaperTime: 2
# If writing transfer functions to file, this sets the name. If transfer functions should not be written to file, this should be set to None
WitnessTFFilename: transfer_functions_DCS_TEST
# Should the transfer function calculation use a median? If not, an average (mean) is used.
-WitnessTFUseMedian: No
+WitnessTFUseMedian: Yes
# Should transfer functions be computed on a fixed schedule, so that the output does not depend on start time? This is useful for running jobs in parallel. Otherwise, they are computed asap.
WitnessTFParallelMode: Yes
# When using parallel mode, how many seconds later should we shift the time when transfer functions start being computed from a multiple of the cycle period?
-WitnessTFTimeShift: 350
+WitnessTFTimeShift: 360
###############################
# Options for HOFT_CLEAN bits #
###############################
diff --git a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_TEST_AllCorrections_Cleaning.ini b/gstlal-calibration/config_files/O2/H1/tests/H1DCS_TEST_AllCorrections_Cleaning.ini
deleted file mode 100644
index 286ecb850a..0000000000
--- a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_TEST_AllCorrections_Cleaning.ini
+++ /dev/null
@@ -1,363 +0,0 @@
-[InputConfigurations]
-# Track what "version" of config file this is, so that the pipeline knows which options are present in this file
-ConfigVersion: 0
-# Filters file containing calibration FIR filters, relative to the directory gstlal-calibration/tests/check_calibration/, from which the pipeline is expected to be run
-FiltersFileName: H1DCS_newsrcline_1173225472.npz
-# Data source should be set to frames or lvshm
-DataSource: frames
-FileChecksum: No
-# Right now, SkipBadFiles needs to be off when reading from frames
-SkipBadFiles: No
-############################################
-# If reading from frames use these options #
-############################################
-# None
-
-###################################################
-# If reading from shared memory use these options #
-###################################################
-SHMPartition: LHO_Online
-# Assumed duration of input frames in seconds
-InputFrameDuration: 1
-
-[OutputConfigurations]
-CompressionScheme: 6
-CompressionLevel: 3
-ChanPrefix: DCS-
-# Set to "None" if you do not want a channel suffix
-ChanSuffix: None
-# Data sink should be set to frames or lvshm
-DataSink: frames
-#################################################
-# If writing to shared memory use these options #
-#################################################
-OutputSHMPartition: hoft_test
-BufferMode: 2
-# Use this to approximate the frame size (in bytes) when writing to shared memory
-FrameSize: 405338
-NumBuffers: 10
-###############################################
-# If writing to frame files use these options #
-###############################################
-FrameType: H1DCS_TEST
-
-[CalibrationConfigurations]
-IFO: H1
-# Set calibration mode to Full or Partial
-CalibrationMode: Full
-ComputeCalibStateVector: Yes
-
-[DebuggingConfigurations]
-# If you want to write a pipeline graph, provide the graph name. Otherwise, set name equal to None
-PipelineGraphFilename: gstlal_compute_strain
-Verbose: No
-# Turn this on to write data presentation timestamps and real-time unix timestamps to file at the beginning and end of the pipeline, to measure latency
-TestLatency: No
-
-[TDCFConfigurations]
-#########################################################
-# Options related to time dependent correction factors #
-#########################################################
-ComputeKappaTST: Yes
-ApplyKappaTST: No
-# Set this to have the \kappa_tst factors filter the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
-ApplyComplexKappaTST: Yes
-
-ComputeKappaPU: Yes
-ApplyKappaPU: No
-# Set this to have the \kappa_pu factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors
-ApplyComplexKappaPU: Yes
-
-ComputeKappaPUM: No
-ApplyKappaPUM: No
-# Set this to have the \kappa_p factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
-ApplyComplexKappaPUM: No
-
-ComputeKappaUIM: No
-ApplyKappaUIM: No
-# Set this to have the \kappa_u factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
-ApplyComplexKappaUIM: No
-
-# Set this to use a calibration line injected using the UIM stage of actuation to compute \kappa_U. Otherwise, the DARM_ctrl line is used.
-UseUIMLine: Yes
-
-ComputeKappaC: Yes
-ApplyKappaC: Yes
-
-ComputeFcc: Yes
-ApplyFcc: Yes
-
-ComputeSRCQ: Yes
-ApplySRCQ: Yes
-
-ComputeFs: Yes
-ApplyFs: Yes
-
-###########################################
-# Options related to the coherence gating #
-###########################################
-UseCoherence: Yes
-CoherenceUncThreshold: 0.004
-# Amount of time used in front-end to compute coherence
-CoherenceTime: 130
-###################################################################
-# Options related to the computation configurations for the TDCFs #
-###################################################################
-ComputeFactorsSR: 16
-RecordFactorsSR: 16
-# Length in seconds of low-pass FIR filter used in demodulation of the calibration lines
-DemodulationFilterTime: 20
-# Time (in seconds) to smooth out \kappas with a median-like method
-MedianSmoothingTime: 128
-TDCFAveragingTime: 10
-#If set to yes, bad computed kappas will be replaced by the previous computed median in the running median array. Otherwise, they are replaced with the default value
-TDCFDefaultToMedian: Yes
-##################################################
-# Options related to updating cavity pole filter #
-##################################################
-# Duration of the Fcc filter in the time domain in seconds
-FccFilterDuration: 0.01
-# Number of seconds to average Fcc values before creating a new Fcc filter
-FccAveragingTime: 60
-# Number of samples to be used when tapering old filter and ramping in new filter
-FccFilterTaperLength: 32768
-############################
-# Nominal values for TDCFs #
-############################
-ExpectedKappaTSTReal: 1.0
-ExpectedKappaTSTImag: 0.0
-ExpectedKappaPUMReal: 1.0
-ExpectedKappaPUMImag: 0.0
-ExpectedKappaUIMReal: 1.0
-ExpectedKappaUIMImag: 0.0
-ExpectedKappaPUReal: 1.0
-ExpectedKappaPUImag: 0.0
-ExpectedKappaC: 1.0
-ExpectedFcc: 360.0
-ExpectedFs: 6.91
-ExpectedSRCQ: 21.739
-################################
-# Acceptable variance in TDCFs #
-################################
-KappaTSTRealVar: 0.2
-KappaTSTImagVar: 0.2
-KappaPURealVar: 0.2
-KappaPUImagVar: 0.2
-KappaPUMRealVar: 0.2
-KappaPUMImagVar: 0.2
-KappaUIMRealVar: 0.2
-KappaUIMImagVar: 0.2
-KappaCVar: 0.2
-FccVar: 50.0
-FsVar: 5.0
-SRCQInvMin: 0.0
-SRCQInvMax: 0.5
-#######################
-# EPICS records input #
-#######################
-# Set to Yes if EPICS records for TDCF computations should be read from filters file. If set to No, they will be read from frames
-FactorsFromFiltersFile: Yes
-##############################################################################################
-# Updating Sensing and Actuation filters with all frequency-dependent corrections parameters #
-##############################################################################################
-# Length of time (in seconds) between when inverse-sensing FIR filter is updated
-SensingFilterUpdateTime: 60
-# Length of time (in seconds) over which the smoothed time-dependent parameters of the sensing function are averaged before updating the filter
-SensingFilterAveragingTime: 1
-# Number of samples to be used when tapering old inverse sensing filter and ramping in new filter
-SensingFilterTaperLength: 32768
-# Length of time (in seconds) between when the actuation FIR filters are updated
-ActuationFilterUpdateTime: 60
-# Length of time (in seconds) over which the smoothed time-dependent parameters of the actuation function are averaged before updating the filter
-ActuationFilterAveragingTime: 1
-# Number of samples to be used when tapering old actuation filters and ramping in new filters
-ActuationFilterTaperLength: 32768
-
-[ChannelNames]
-#############################
-# Calibration Channel Names #
-#############################
-DARMCtrlChannel: CAL-DARM_CTRL_WHITEN_OUT_DBL_DQ
-DARMErrChannel: CAL-DARM_ERR_WHITEN_OUT_DBL_DQ
-DeltaLTSTChannel: CAL-DELTAL_CTRL_TST_DBL_DQ
-DeltaLPUMChannel: CAL-DELTAL_CTRL_PUM_DBL_DQ
-DeltaLUIMChannel: CAL-DELTAL_CTRL_UIM_DBL_DQ
-DeltaLResChannel: CAL-DELTAL_RESIDUAL_DBL_DQ
-####################################
-# Data Quality Vector Channel Name #
-####################################
-InputDQChannel: ODC-MASTER_CHANNEL_OUT_DQ
-##################################
-# Calibration Line Channel Names #
-##################################
-DARMExcChannel: CAL-CS_LINE_SUM_DQ
-TSTExcChannel: SUS-ETMY_L3_CAL_LINE_OUT_DQ
-PUMExcChannel: SUS-ETMY_L2_CAL_LINE_OUT_DQ
-UIMExcChannel: SUS-ETMY_L1_CAL_LINE_OUT_DQ
-PCALChannel: CAL-PCALY_TX_PD_OUT_DQ
-#######################################
-# Coherence Uncertainty Channel Names #
-#######################################
-CohUncSusLine1Channel: CAL-CS_TDEP_SUS_LINE1_UNCERTAINTY
-CohUncSusLine2Channel: CAL-CS_TDEP_SUS_LINE2_UNCERTAINTY
-CohUncSusLine3Channel: CAL-CS_TDEP_SUS_LINE1_UNCERTAINTY
-CohUncPcalyLine1Channel: CAL-CS_TDEP_PCALY_LINE1_UNCERTAINTY
-CohUncPcalyLine2Channel: CAL-CS_TDEP_PCALY_LINE2_UNCERTAINTY
-CohUncDARMLine1Channel: CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY
-###################################
-# Noise Subtraction Channel Names #
-###################################
-PowerLinesChannel: PEM-EY_MAINSMON_EBAY_1_DQ
-# Comma-separated list of witness channels to use to subtract noise from h(t)
-# Set to None if no witness channels are to be used
-WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ;PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ
-###############################
-# EPICS Records Channel Names #
-###############################
-EP1RealChannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_REAL
-EP1ImagChannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_IMAG
-EP2RealChannel: CAL-CS_TDEP_REF_CLGRATIO_CTRL_REAL
-EP2ImagChannel: CAL-CS_TDEP_REF_CLGRATIO_CTRL_IMAG
-EP3RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_REAL
-EP3ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_IMAG
-EP4RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_TST_REAL
-EP4ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_TST_IMAG
-EP5RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_REAL
-EP5ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_IMAG
-EP6RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_REAL
-EP6ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_IMAG
-EP7RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_D_REAL
-EP7ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_D_IMAG
-EP8RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_REAL
-EP8Imagchannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_IMAG
-EP9RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_REAL
-EP9ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_IMAG
-EP10RealChannel: CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_REAL
-EP10ImagChannel: CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_IMAG
-EP11RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_REAL
-EP11ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_IMAG
-EP12RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_D_REAL
-EP12ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_D_IMAG
-EP13RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_REAL
-EP13ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_IMAG
-EP14RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_REAL
-EP14ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_IMAG
-EP15RealChannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_PUM_REAL
-EP15Imagchannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_PUM_IMAG
-EP16RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_UIM_INV_REAL
-EP16ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_UIM_INV_IMAG
-EP17RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_PUM_REAL
-EP17ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_PUM_IMAG
-EP18RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_PUM_REAL
-EP18ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_PUM_IMAG
-EP19RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_UIM_REAL
-EP19ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_UIM_IMAG
-EP20RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_PUM_REAL
-EP20ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_PUM_IMAG
-EP21RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_UIM_REAL
-EP21ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_UIM_IMAG
-EP22RealChannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_UIM_REAL
-EP22ImagChannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_UIM_IMAG
-EP23RealChannel: CAL-CS_TDEP_PUM_LINE1_REF_A_PUM_NOLOCK_REAL
-EP23ImagChannel: CAL-CS_TDEP_PUM_LINE1_REF_A_PUM_NOLOCK_IMAG
-EP24RealChannel: CAL-CS_TDEP_UIM_LINE1_REF_A_UIM_NOLOCK_REAL
-EP24ImagChannel: CAL-CS_TDEP_UIM_LINE1_REF_A_UIM_NOLOCK_IMAG
-
-[SampleRates]
-# Sample rate at which to compute h(t)
-HoftSR: 16384
-# Sample rate at which to compute CALIB_STATE_VECTOR
-CalibStateSR: 16
-# Sample rate of control channel
-# Should be 16384 if using DARM_CTRL and 4096 if using DELTAL_CTRL
-CtrlSR: 16384
-# Sample rate of ODC channel
-ODCSR: 16384
-# Sample rate of TST excitation channel
-TSTExcSR: 512
-# Sample rate of PUM excitation channel
-PUMExcSR: 512
-# Sample rate of UIM excitation channel
-UIMExcSR: 512
-# Sample rate of coherence channels
-CohSR: 16
-# Sample rate for the EPICS reference channels
-EPICSRefSR: 16
-# Sample rate for power lines channel
-PowerLinesChannelSR: 1024
-# Sample rates at which transfer functions will be computed and witness channels will be filtered, given as a semicolon-separated list, e.g., 2048;2048;512;2048. This must be given if WitnessChannelList is not None, and it must be the same length.
-WitnessChannelSR: 2048;2048
-# Sample rates at which to compute and record TDCFs
-ComputeFactorsSR: 16
-RecordFactorsSR: 16
-
-[Bitmasks]
-ObsReadyBitmask: 4
-ObsIntentBitmask: 2
-CBCHWInjBitmask: 16777216
-BurstHWInjBitmask: 33554432
-DetCharHWInjBitmask: 67108864
-StochHWInjBitmask: 8388608
-
-[PipelineConfigurations]
-BufferLength: 1.0
-FrequencyDomainFiltering: No
-Dewhitening: No
-# Latency of all filtering/averaging/median processes (other than calibration model filters) as a fraction of filter length. Value should be set between 0.0 and 1.0.
-FilterLatency: 0.0
-
-[DataCleaningConfigurations]
-####################################################
-# Options for turning on and off line subtraction #
-####################################################
-# Remove the DC component from the residual and control channels before filtering
-RemoveDC: No
-# Subtract the calibration lines from the h(t) spectrum
-RemoveCalLines: Yes
-# Subtract the power lines from the h(t) spectrum
-RemovePowerLines: Yes
-####################################################
-# Options for running power mains line subtraction #
-####################################################
-# Amount by which frequency of power lines varies with time
-PowerLinesFreqVar: 0.02
-# Time over which to average the transfer function between the power mains witness channel and h(t) at 60 Hz and harmonics
-PowerLinesTFAveragingTime: 128
-#######################################
-# Options for broadband noise removal #
-#######################################
-# The length in seconds of the fast Fourier transforms used to compute transfer functions between witness channels and h(t). The fft's are windowed with Hann windows and overlapped.
-WitnessChannelFFTTime: 4.0
-# The number of ffts to take before averaging the witness -> h(t) transfer functions calculation. The average is taken after the ratio h(f) / witness(f).
-NumWitnessFFTs: 1800
-# Sets the minimum number of FFTs necessary to produce the first transfer functions and clean data after data flow starts.
-MinWitnessFFTs: 400
-# The length in seconds of the filters applied to the witness channels before subtracting from h(t)
-WitnessFIRLength: 0.5
-# The frequency resolution of the filters applied to the witness channels before subtracting from h(t). It can be advantageous to lower the frequency resolution in order to average over excess noise.
-WitnessFrequencyResolution: 1.0
-# List of minima and maxima of frequency ranges where the Fourier transform of h(t) will be replaced by a straight line in the calculation of transfer functions between witness channels and h(t) for noise subtraction. Semicolons separate lists for different sets of witness channels. If no notches are desired, use zeros, e.g., \'0;0;0\'. Here is an example using the expected format: \'495.0,515.0,985.0,1015.0;59,60,119,121;0\' This can be useful, e.g., if there are loud lines in the signal that are not present in the witness channels.
-WitnessNotchFrequencies: 495.0,515.0,985.0,1015.0;495.0,515.0,985.0,1015.0
-# The amount of time after transfer functions between witness channels and h(t) are finished to begin the calculation of the next set of transfer functions
-WitnessTFUpdateTime: 3600
-# If lock-loss lasts at least this many seconds, transfer functions will revert to those computed at the beginning of a lock stretch. Set to zero to disbale.
-CriticalLockLossTime: 1800
-# The amount of time to use to taper in newly computed FIR filters for witness channels being used for noise subtraction.
-WitnessFilterTaperTime: 10
-# If writing transfer functions to file, this sets the name. If transfer functions should not be written to file, this should be set to None
-WitnessTFFilename: transfer_functions
-# Should the transfer function calculation use a median? If not, an average (mean) is used.
-WitnessTFUseMedian: Yes
-###############################
-# Options for HOFT_CLEAN bits #
-###############################
-# The amount of data from h(t) and cleaned h(t) that is used to compute and compare the rms. This comparison between cleaned and uncleaned h(t) determines whether the HOFT_CLEAN bits of the calibration state vector are on or off.
-CleaningCheckRMSTime: 20.0
-# Minimum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_LOWFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
-CleaningCheckRangeLowMin: 15
-# Maximum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_LOWFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
-CleaningCheckRangeLowMax: 40
-# Minimum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_MIDFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
-CleaningCheckRangeMidMin: 100
-# Maximum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_MIDFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
-CleaningCheckRangeMidMax: 500
diff --git a/gstlal-calibration/config_files/O2/H1/tests/H1GDS_LowLatency_AllCorrections_Cleaning.ini b/gstlal-calibration/config_files/O2/H1/tests/H1GDS_LowLatency_AllCorrections_Cleaning.ini
index 5d08ce9f4e..65af1e997b 100644
--- a/gstlal-calibration/config_files/O2/H1/tests/H1GDS_LowLatency_AllCorrections_Cleaning.ini
+++ b/gstlal-calibration/config_files/O2/H1/tests/H1GDS_LowLatency_AllCorrections_Cleaning.ini
@@ -1,6 +1,4 @@
[InputConfigurations]
-# Track what "version" of config file this is, so that the pipeline knows which options are present in this file
-ConfigVersion: 0
# Filters file containing calibration FIR filters
FiltersFileName: H1GDS_minlatency_1175954418.npz
# Data source should be set to frames or lvshm
diff --git a/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.c b/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.c
index 0e1d9b2967..330b9d841b 100644
--- a/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.c
+++ b/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.c
@@ -112,6 +112,8 @@ enum property {
ARG_ADAPTIVE_FILTER_LENGTH,
ARG_TUKEY_PARAM,
ARG_FILTER_SAMPLE_RATE,
+ ARG_FILTER_TIMESHIFT,
+ ARG_FILTER_ENDTIME,
ARG_WRITE_TO_SCREEN,
ARG_FILENAME,
ARG_FAKE
@@ -419,6 +421,14 @@ static void average_input_data_ ## DTYPE(GSTLALAdaptiveFIRFilt *element, complex
\
/* Let other elements know about the update */ \
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_ADAPTIVE_FILTER]); \
+ /* Provide a timestamp indicating when the filter becomes invalid if requested */ \
+ if(element->filter_timeshift < G_MAXINT64) { \
+ if(element->filter_timeshift < 0 && (guint64) (-element->filter_timeshift) > pts) \
+ element->filter_endtime = 0; \
+ else \
+ element->filter_endtime = pts + element->filter_timeshift; \
+ g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FILTER_ENDTIME]); \
+ } \
\
/* Write FIR filter to the screen or a file if we want */ \
if(element->write_to_screen || element->filename) { \
@@ -507,6 +517,29 @@ static gboolean start(GstBaseSink *sink) {
}
+/*
+ * event()
+ */
+
+
+static gboolean event(GstBaseSink *sink, GstEvent *event) {
+
+ GSTLALAdaptiveFIRFilt *element = GSTLAL_ADAPTIVEFIRFILT(sink);
+ gboolean success;
+ GST_DEBUG_OBJECT(element, "Got %s event on sink pad", GST_EVENT_TYPE_NAME(event));
+
+ if(GST_EVENT_TYPE(event) == GST_EVENT_EOS && element->filter_timeshift < G_MAXINT64) {
+ /* These filters should remain usable as long as possible */
+ element->filter_endtime = G_MAXUINT64 - 1;
+ g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FILTER_ENDTIME]);
+ }
+
+ success = GST_BASE_SINK_CLASS(gstlal_adaptivefirfilt_parent_class)->event(sink, event);
+
+ return success;
+}
+
+
/*
* set_caps()
*/
@@ -779,6 +812,10 @@ static void set_property(GObject *object, enum property id, const GValue *value,
element->filter_sample_rate = g_value_get_int(value);
break;
+ case ARG_FILTER_TIMESHIFT:
+ element->filter_timeshift = g_value_get_int64(value);
+ break;
+
case ARG_WRITE_TO_SCREEN:
element->write_to_screen = g_value_get_boolean(value);
break;
@@ -871,6 +908,14 @@ static void get_property(GObject *object, enum property id, GValue *value, GPara
g_value_set_int(value, element->filter_sample_rate);
break;
+ case ARG_FILTER_TIMESHIFT:
+ g_value_set_int64(value, element->filter_timeshift);
+ break;
+
+ case ARG_FILTER_ENDTIME:
+ g_value_set_uint64(value, element->filter_endtime);
+ break;
+
case ARG_WRITE_TO_SCREEN:
g_value_set_boolean(value, element->write_to_screen);
break;
@@ -927,6 +972,7 @@ static void gstlal_adaptivefirfilt_class_init(GSTLALAdaptiveFIRFiltClass *klass)
gstbasesink_class->set_caps = GST_DEBUG_FUNCPTR(set_caps);
gstbasesink_class->start = GST_DEBUG_FUNCPTR(start);
+ gstbasesink_class->event = GST_DEBUG_FUNCPTR(event);
gstbasesink_class->render = GST_DEBUG_FUNCPTR(render);
gstbasesink_class->stop = GST_DEBUG_FUNCPTR(stop);
@@ -1113,6 +1159,26 @@ static void gstlal_adaptivefirfilt_class_init(GSTLALAdaptiveFIRFiltClass *klass)
1, G_MAXINT, 16384,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
);
+ properties[ARG_FILTER_TIMESHIFT] = g_param_spec_int64(
+ "filter-timeshift",
+ "Filter time-shift",
+ "The number of nanoseconds after the completion of a FIR filter calculation\n\t\t\t"
+ "that the FIR filter remains valid for use on the filtered data. This is\n\t\t\t"
+ "added to the presentation timestamp when the filter is completed to compute\n\t\t\t"
+ "the filter-endtime property. Default is to disable.",
+ G_MININT64, G_MAXINT64, G_MAXINT64,
+ G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+ );
+ properties[ARG_FILTER_ENDTIME] = g_param_spec_uint64(
+ "filter-endtime",
+ "Filter end time",
+ "The time when a computed FIR filter ceases to be valid for use on\n\t\t\t"
+ "filtered data. This can be compared to the presentation timestamps of the\n\t\t\t"
+ "filtered data to determine whether the filter is still valid. Default is\n\t\t\t"
+ "to disable.",
+ 0, G_MAXUINT64, G_MAXUINT64,
+ G_PARAM_READABLE | G_PARAM_STATIC_STRINGS
+ );
properties[ARG_WRITE_TO_SCREEN] = g_param_spec_boolean(
"write-to-screen",
"Write to Screen",
@@ -1200,6 +1266,16 @@ static void gstlal_adaptivefirfilt_class_init(GSTLALAdaptiveFIRFiltClass *klass)
ARG_FILTER_SAMPLE_RATE,
properties[ARG_FILTER_SAMPLE_RATE]
);
+ g_object_class_install_property(
+ gobject_class,
+ ARG_FILTER_TIMESHIFT,
+ properties[ARG_FILTER_TIMESHIFT]
+ );
+ g_object_class_install_property(
+ gobject_class,
+ ARG_FILTER_ENDTIME,
+ properties[ARG_FILTER_ENDTIME]
+ );
g_object_class_install_property(
gobject_class,
ARG_WRITE_TO_SCREEN,
diff --git a/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.h b/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.h
index b069998812..2c3a5dd8d6 100644
--- a/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.h
+++ b/gstlal-calibration/gst/lal/gstlal_adaptivefirfilt.h
@@ -98,6 +98,8 @@ struct _GSTLALAdaptiveFIRFilt {
double *tukey;
gint64 tukey_length;
gint filter_sample_rate;
+ gint64 filter_timeshift;
+ guint64 filter_endtime;
gboolean write_to_screen;
char *filename;
};
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.c b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
index f5fd051264..cf50d1dbb7 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.c
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
@@ -126,8 +126,10 @@ enum property {
ARG_HIGH_PASS,
ARG_LOW_PASS,
ARG_NOTCH_FREQUENCIES,
+ ARG_FIR_TIMESHIFT,
ARG_TRANSFER_FUNCTIONS,
ARG_FIR_FILTERS,
+ ARG_FIR_ENDTIME,
ARG_FAKE
};
@@ -860,12 +862,12 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
} \
element->num_tfs_since_gap++; \
} else if (element->parallel_mode) \
- GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. Trying again."); \
- else \
GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. Waiting for the next cycle."); \
+ else \
+ GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. Trying again."); \
\
if(num_ffts_in_avg_if_nogap == element->num_ffts) { \
- element->sample_count = (gint64) (gst_util_uint64_scale_int_round(pts, element->rate, GST_SECOND) + src_size - element->update_delay_samples) % (element->update_samples + element->num_ffts * stride + element->fft_overlap); \
+ element->sample_count = (gint64) (gst_util_uint64_scale_int_round(pts, element->rate, GST_SECOND) + src_size + element->update_samples - element->update_delay_samples) % (element->update_samples + element->num_ffts * stride + element->fft_overlap); \
if(element->sample_count > element->update_samples) \
element->sample_count -= element->update_samples + element->num_ffts * stride + element->fft_overlap; \
element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg = 0; \
@@ -876,9 +878,17 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
if(success && element->make_fir_filters) { \
success &= update_fir_filters_ ## DTYPE(element->transfer_functions, num_tfs, element->fir_length, element->rate, element->workspace.w ## S_OR_D ## pf.fir_filter, element->workspace.w ## S_OR_D ## pf.fir_plan, element->workspace.w ## S_OR_D ## pf.fir_window, element->workspace.w ## S_OR_D ## pf.tukey, element->fir_filters); \
if(success) { \
- GST_INFO_OBJECT(element, "Just computed new FIR filters"); \
+ GST_LOG_OBJECT(element, "Just computed new FIR filters"); \
/* Let other elements know about the update */ \
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]); \
+ /* Provide a timestamp indicating when the filter becomes invalid if requested */ \
+ if(element->fir_timeshift < G_MAXINT64) { \
+ if(element->fir_timeshift < 0 && (guint64) (-element->fir_timeshift) > pts) \
+ element->fir_endtime = 0; \
+ else \
+ element->fir_endtime = pts + element->fir_timeshift; \
+ g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_ENDTIME]); \
+ } \
/* Write FIR filters to the screen or a file if we want */ \
if(element->write_to_screen || element->filename) \
write_fir_filters(element->fir_filters, gst_element_get_name(element), element->fir_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (num_ffts_in_avg_if_nogap * stride + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE); \
@@ -887,7 +897,9 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
for(i = 0; i < num_tfs * element->fir_length; i++) \
element->post_gap_fir_filters[i] = element->fir_filters[i]; \
} \
- } else \
+ } else if (element->parallel_mode) \
+ GST_WARNING_OBJECT(element, "FIR filter(s) computation failed. Waiting for the next cycle."); \
+ else \
GST_WARNING_OBJECT(element, "FIR filter(s) computation failed. Trying again."); \
} \
} \
@@ -1025,7 +1037,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
}
success &= update_transfer_functions_float(element->workspace.wspf.autocorrelation_matrix, num_tfs, fd_fft_length, fd_tf_length, element->workspace.wspf.sinc_table, element->workspace.wspf.sinc_length, element->workspace.wspf.sinc_taps_per_df, element->use_median ? 1 : element->num_ffts - element->workspace.wspf.num_ffts_dropped, element->workspace.wspf.transfer_functions_at_f, element->workspace.wspf.transfer_functions_solved_at_f, element->workspace.wspf.autocorrelation_matrix_at_f, element->workspace.wspf.permutation, element->transfer_functions);
if(success) {
- GST_INFO_OBJECT(element, "Just computed new transfer functions");
+ GST_LOG_OBJECT(element, "Just computed new transfer functions");
/* Let other elements know about the update */
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
@@ -1037,7 +1049,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
if(success && element->make_fir_filters) {
success &= update_fir_filters_float(element->transfer_functions, num_tfs, element->fir_length, element->rate, element->workspace.wspf.fir_filter, element->workspace.wspf.fir_plan, element->workspace.wspf.fir_window, element->workspace.wspf.tukey, element->fir_filters);
if(success) {
- GST_INFO_OBJECT(element, "Just computed new FIR filters");
+ GST_LOG_OBJECT(element, "Just computed new FIR filters");
/* Let other elements know about the update */
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
/* Write FIR filters to the screen or a file if we want */
@@ -1100,7 +1112,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
}
success &= update_transfer_functions_double(element->workspace.wdpf.autocorrelation_matrix, num_tfs, fd_fft_length, fd_tf_length, element->workspace.wdpf.sinc_table, element->workspace.wdpf.sinc_length, element->workspace.wdpf.sinc_taps_per_df, element->use_median ? 1 : element->num_ffts - element->workspace.wdpf.num_ffts_dropped, element->workspace.wdpf.transfer_functions_at_f, element->workspace.wdpf.transfer_functions_solved_at_f, element->workspace.wdpf.autocorrelation_matrix_at_f, element->workspace.wdpf.permutation, element->transfer_functions);
if(success) {
- GST_INFO_OBJECT(element, "Just computed new transfer functions");
+ GST_LOG_OBJECT(element, "Just computed new transfer functions");
/* Let other elements know about the update */
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
@@ -1112,7 +1124,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
if(success && element->make_fir_filters) {
success &= update_fir_filters_double(element->transfer_functions, num_tfs, element->fir_length, element->rate, element->workspace.wdpf.fir_filter, element->workspace.wdpf.fir_plan, element->workspace.wdpf.fir_window, element->workspace.wdpf.tukey, element->fir_filters);
if(success) {
- GST_INFO_OBJECT(element, "Just computed new FIR filters");
+ GST_LOG_OBJECT(element, "Just computed new FIR filters");
/* Let other elements know about the update */
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
/* Write FIR filters to the screen or a file if we want */
@@ -1125,6 +1137,12 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
}
}
+ if(GST_EVENT_TYPE(event) == GST_EVENT_EOS && element->fir_timeshift < G_MAXINT64) {
+ /* These filters should remain usable as long as possible */
+ element->fir_endtime = G_MAXUINT64 - 1;
+ g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_ENDTIME]);
+ }
+
success = GST_BASE_SINK_CLASS(gstlal_transferfunction_parent_class)->event(sink, event);
return success;
@@ -1735,7 +1753,7 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
element->offset0 = GST_BUFFER_OFFSET(buffer);
if(element->parallel_mode) {
/* In this case, we want to compute the transfer functions on a schedule, not asap */
- element->sample_count = (gint64) (gst_util_uint64_scale_int_round(element->t0, element->rate, GST_SECOND) - element->update_delay_samples) % (element->update_samples + element->num_ffts * (element->fft_length - element->fft_overlap) + element->fft_overlap);
+ element->sample_count = (gint64) (gst_util_uint64_scale_int_round(element->t0, element->rate, GST_SECOND) + element->update_samples - element->update_delay_samples) % (element->update_samples + element->num_ffts * (element->fft_length - element->fft_overlap) + element->fft_overlap);
if(element->sample_count > element->update_samples)
element->sample_count -= element->update_samples + element->num_ffts * (element->fft_length - element->fft_overlap) + element->fft_overlap;
} else if(element->sample_count > element->update_samples) {
@@ -1774,7 +1792,7 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
gint64 i;
for(i = 0; i < (element->channels - 1) * (element->fir_length / 2 + 1); i++)
element->transfer_functions[i] = element->post_gap_transfer_functions[i];
- GST_INFO_OBJECT(element, "Just reverted to post-gap transfer functions");
+ GST_LOG_OBJECT(element, "Just reverted to post-gap transfer functions");
/* Let other elements know about the update */
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
@@ -1783,7 +1801,7 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
if(element->make_fir_filters) {
for(i = 0; i < (element->channels - 1) * element->fir_length; i++)
element->fir_filters[i] = element->post_gap_fir_filters[i];
- GST_INFO_OBJECT(element, "Just reverted to post-gap FIR filters");
+ GST_LOG_OBJECT(element, "Just reverted to post-gap FIR filters");
/* Let other elements know about the update */
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
/* Write FIR filters to the screen or a file if we want */
@@ -2086,6 +2104,10 @@ static void set_property(GObject *object, enum property id, const GValue *value,
element->num_notches /= 2;
break;
+ case ARG_FIR_TIMESHIFT:
+ element->fir_timeshift = g_value_get_int64(value);
+ break;
+
case ARG_TRANSFER_FUNCTIONS:
if(element->transfer_functions) {
g_free(element->transfer_functions);
@@ -2192,6 +2214,10 @@ static void get_property(GObject *object, enum property id, GValue *value, GPara
g_value_take_boxed(value, gstlal_g_value_array_from_doubles(element->notch_frequencies, 2 * element->num_notches));
break;
+ case ARG_FIR_TIMESHIFT:
+ g_value_set_int64(value, element->fir_timeshift);
+ break;
+
case ARG_TRANSFER_FUNCTIONS:
if(element->transfer_functions) {
GValueArray *va;
@@ -2222,6 +2248,10 @@ static void get_property(GObject *object, enum property id, GValue *value, GPara
}
break;
+ case ARG_FIR_ENDTIME:
+ g_value_set_uint64(value, element->fir_endtime);
+ break;
+
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, id, pspec);
break;
@@ -2476,6 +2506,16 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
),
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
);
+ properties[ARG_FIR_TIMESHIFT] = g_param_spec_int64(
+ "fir-timeshift",
+ "FIR time-shift",
+ "The number of nanoseconds after the completion of a FIR filter calculation\n\t\t\t"
+ "that the FIR filter remains valid for use on the filtered data. This is\n\t\t\t"
+ "added to the presentation timestamp when the filter is completed to compute\n\t\t\t"
+ "the fir-endtime property. Default is to disable.",
+ G_MININT64, G_MAXINT64, G_MAXINT64,
+ G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+ );
properties[ARG_TRANSFER_FUNCTIONS] = g_param_spec_value_array(
"transfer-functions",
"Transfer Functions",
@@ -2514,6 +2554,16 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
),
G_PARAM_READABLE | G_PARAM_STATIC_STRINGS
);
+ properties[ARG_FIR_ENDTIME] = g_param_spec_uint64(
+ "fir-endtime",
+ "FIR end time",
+ "The time when a computed FIR filter ceases to be valid for use on\n\t\t\t"
+ "filtered data. This can be compared to the presentation timestamps of the\n\t\t\t"
+ "filtered data to determine whether the filter is still valid. Default is\n\t\t\t"
+ "to disable.",
+ 0, G_MAXUINT64, G_MAXUINT64,
+ G_PARAM_READABLE | G_PARAM_STATIC_STRINGS
+ );
g_object_class_install_property(
@@ -2606,6 +2656,11 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
ARG_NOTCH_FREQUENCIES,
properties[ARG_NOTCH_FREQUENCIES]
);
+ g_object_class_install_property(
+ gobject_class,
+ ARG_FIR_TIMESHIFT,
+ properties[ARG_FIR_TIMESHIFT]
+ );
g_object_class_install_property(
gobject_class,
ARG_TRANSFER_FUNCTIONS,
@@ -2616,6 +2671,11 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
ARG_FIR_FILTERS,
properties[ARG_FIR_FILTERS]
);
+ g_object_class_install_property(
+ gobject_class,
+ ARG_FIR_ENDTIME,
+ properties[ARG_FIR_ENDTIME]
+ );
}
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.h b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
index 1f80a41a6e..c3705ec4cb 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.h
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
@@ -165,10 +165,12 @@ struct _GSTLALTransferFunction {
double *notch_frequencies;
gint64 *notch_indices;
int num_notches;
+ gint64 fir_timeshift;
complex double *post_gap_transfer_functions;
double *post_gap_fir_filters;
complex double *transfer_functions;
double *fir_filters;
+ guint64 fir_endtime;
};
diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py
index 0def9ba8dc..7bb7c7a1bd 100644
--- a/gstlal-calibration/python/calibration_parts.py
+++ b/gstlal-calibration/python/calibration_parts.py
@@ -75,25 +75,25 @@ def mkcomplexfirbank2(pipeline, src, latency = None, fir_matrix = None, time_dom
def mkpow(pipeline, src, **properties):
return pipeparts.mkgeneric(pipeline, src, "cpow", **properties)
-def mkmultiplier(pipeline, srcs, sync = True, queue_length = 0):
+def mkmultiplier(pipeline, srcs, sync = True, queue_length = [0]):
elem = pipeparts.mkgeneric(pipeline, None, "lal_adder", sync=sync, mix_mode="product")
if srcs is not None:
- for src in srcs:
- mkqueue(pipeline, src, length = queue_length).link(elem)
+ for i in range(0, len(srcs)):
+ mkqueue(pipeline, srcs[i], length = queue_length[min(i, len(queue_length) - 1)]).link(elem)
return elem
-def mkadder(pipeline, srcs, sync = True, queue_length = 0):
+def mkadder(pipeline, srcs, sync = True, queue_length = [0]):
elem = pipeparts.mkgeneric(pipeline, None, "lal_adder", sync=sync)
if srcs is not None:
- for src in srcs:
- mkqueue(pipeline, src, length = queue_length).link(elem)
+ for i in range(0, len(srcs)):
+ mkqueue(pipeline, srcs[i], length = queue_length[min(i, len(queue_length) - 1)]).link(elem)
return elem
def mkgate(pipeline, src, control, threshold, queue_length = 0, **properties):
elem = pipeparts.mkgate(pipeline, mkqueue(pipeline, src, length = queue_length), control = mkqueue(pipeline, control, length = queue_length), threshold = threshold, **properties)
return elem
-def mkinterleave(pipeline, srcs, complex_data = False):
+def mkinterleave(pipeline, srcs, complex_data = False, queue_length = [0]):
complex_factor = 1 + int(complex_data)
num_srcs = complex_factor * len(srcs)
i = 0
@@ -103,7 +103,7 @@ def mkinterleave(pipeline, srcs, complex_data = False):
matrix[0][i] = 1
mixed_srcs.append(pipeparts.mkmatrixmixer(pipeline, src, matrix=matrix))
i += complex_factor
- elem = mkadder(pipeline, tuple(mixed_srcs))
+ elem = mkadder(pipeline, tuple(mixed_srcs), queue_length = queue_length)
#chan1 = pipeparts.mkmatrixmixer(pipeline, src1, matrix=[[1,0]])
#chan2 = pipeparts.mkmatrixmixer(pipeline, src2, matrix=[[0,1]])
@@ -172,6 +172,9 @@ def demodulate(pipeline, head, freq, td, rate, filter_time, filter_latency, pref
head = pipeparts.mkgeneric(pipeline, head, "lal_demodulate", line_frequency = freq, prefactor_real = prefactor_real, prefactor_imag = prefactor_imag)
head = mkresample(pipeline, head, 5, filter_latency == 0.0, rate)
+ if filter_latency != 0:
+ # Remove the first several seconds of output, which depend on start time
+ head = pipeparts.mkgeneric(pipeline, head, "lal_insertgap", chop_length = 7000000000)
head = lowpass(pipeline, head, rate, length = filter_time, fcut = 0, filter_latency = filter_latency, td = td)
return head
@@ -204,7 +207,6 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics,
upsample_quality = 4
resample_shift = 16.0 + 16.5
zero_latency = filter_latency == 0.0
- filter_latency = min(0.5, filter_latency)
witness = pipeparts.mktee(pipeline, witness)
signal = pipeparts.mktee(pipeline, signal)
@@ -912,6 +914,10 @@ def compute_Xi_split_act(pipeline, pcalfpcal4, darmfpcal4, fpcal4, EP11, EP12, E
return Xi
+def update_property_simple(prop_maker, arg, prop_taker, maker_prop_name, taker_prop_name):
+ prop = prop_maker.get_property(maker_prop_name)
+ prop_taker.set_property(taker_prop_name, prop)
+
def update_filter(filter_maker, arg, filter_taker, maker_prop_name, taker_prop_name):
firfilter = filter_maker.get_property(maker_prop_name)[::-1]
filter_taker.set_property(taker_prop_name, firfilter)
@@ -920,7 +926,7 @@ def update_filters(filter_maker, arg, filter_taker, maker_prop_name, taker_prop_
firfilter = filter_maker.get_property(maker_prop_name)[filter_number][::-1]
filter_taker.set_property(taker_prop_name, firfilter)
-def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, min_ffts, update_samples, fir_length, frequency_resolution, filter_taper_length, use_median = False, parallel_mode = False, notch_frequencies = [], noisesub_gate_bit = None, delay_time = 0.0, wait_time = 0.0, critical_lock_loss_time = 0, filename = None):
+def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, min_ffts, update_samples, fir_length, frequency_resolution, filter_taper_length, use_median = False, parallel_mode = False, notch_frequencies = [], noisesub_gate_bit = None, delay_time = 0.0, critical_lock_loss_time = 0, filename = None):
#
# Use witness channels that monitor the environment to remove environmental noise
@@ -928,8 +934,6 @@ def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_lengt
# between witness channels.
#
- default_fir_filter = numpy.zeros(fir_length)
-
signal_tee = pipeparts.mktee(pipeline, signal)
witnesses = list(witnesses)
witness_tees = []
@@ -942,11 +946,16 @@ def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_lengt
transfer_functions = mkinterleave(pipeline, numpy.insert(witness_tees, 0, resampled_signal, axis = 0))
if noisesub_gate_bit is not None:
transfer_functions = mkgate(pipeline, transfer_functions, noisesub_gate_bit, 1)
- transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, min_ffts = min_ffts, update_samples = update_samples, make_fir_filters = -1, fir_length = fir_length, frequency_resolution = frequency_resolution, high_pass = 15, update_after_gap = True, use_median = use_median, parallel_mode = parallel_mode, notch_frequencies = notch_frequencies, use_first_after_gap = critical_lock_loss_time * witness_rate, update_delay_samples = int(delay_time * witness_rate), filename = filename)
+ transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, min_ffts = min_ffts, update_samples = update_samples, make_fir_filters = -1, fir_length = fir_length, frequency_resolution = frequency_resolution, high_pass = 15, update_after_gap = True, use_median = use_median, parallel_mode = parallel_mode, notch_frequencies = notch_frequencies, use_first_after_gap = critical_lock_loss_time * witness_rate, update_delay_samples = int(delay_time * witness_rate), fir_timeshift = 0, filename = filename)
signal_minus_noise = [signal_tee]
for i in range(0, len(witnesses)):
- minus_noise = pipeparts.mkgeneric(pipeline, mkqueue(pipeline, witness_tees[i], min_length = wait_time), "lal_tdwhiten", kernel = default_fir_filter, latency = fir_length / 2, taper_length = filter_taper_length)
- transfer_functions.connect("notify::fir-filters", update_filters, minus_noise, "fir_filters", "kernel", i)
+ if parallel_mode:
+ minus_noise = pipeparts.mkgeneric(pipeline, mkqueue(pipeline, witness_tees[i]), "lal_tdwhiten", kernel = numpy.zeros(fir_length), latency = fir_length / 2, taper_length = filter_taper_length, kernel_endtime = 0)
+ transfer_functions.connect("notify::fir-filters", update_filters, minus_noise, "fir_filters", "kernel", i)
+ transfer_functions.connect("notify::fir-endtime", update_property_simple, minus_noise, "fir_endtime", "kernel_endtime")
+ else:
+ minus_noise = pipeparts.mkgeneric(pipeline, witness_tees[i], "lal_tdwhiten", kernel = numpy.zeros(fir_length), latency = fir_length / 2, taper_length = filter_taper_length)
+ transfer_functions.connect("notify::fir-filters", update_filters, minus_noise, "fir_filters", "kernel", i)
signal_minus_noise.append(mkresample(pipeline, minus_noise, 5, False, signal_rate))
return mkadder(pipeline, tuple(signal_minus_noise))
diff --git a/gstlal-calibration/tests/check_calibration/Makefile b/gstlal-calibration/tests/check_calibration/Makefile
index 6495718750..fe6e06d375 100644
--- a/gstlal-calibration/tests/check_calibration/Makefile
+++ b/gstlal-calibration/tests/check_calibration/Makefile
@@ -6,16 +6,19 @@
# which interferometer (H or L)
IFO = H
# determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
-OBSRUN = PreER13
+OBSRUN = O2
-START = 1225967424
+START = $(shell echo 1185767424 - 715 | bc)
+#1225967424
#$(shell echo 1185763328 - 700 | bc)
-END = 1225968448
+END = $(shell echo 1185771520 + 715 | bc)
+#1225968448
#$(shell echo 1185767424 + 700 | bc)
+# 1185771520
SHMRUNTIME = 400
# How much time does the calibration need to settle at the start and end?
-PLOT_WARMUP_TIME = 256
-PLOT_COOLDOWN_TIME = 64
+PLOT_WARMUP_TIME = 715
+PLOT_COOLDOWN_TIME = 715
GDSCONFIGS = ../../config_files/PreER13/H1/tests/H1GDS_1227048227.ini
DCSCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
@@ -24,7 +27,7 @@ GDSTESTCONFIGS = ../../config_files/PreER13/H1/H1GDS_TEST_1225558818.ini
DCSTESTCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning_TEST.ini
GDSSHMCONFIGS = ../../config_files/PreER13/H1/tests/H1GDS_1222058826_shm2frames.ini
-all: CALCS_GDS_ASD GDS_over_CALCS
+all: noise_subtraction_ASD_DCS_TEST
###############################################
### These commands should change less often ###
diff --git a/gstlal-calibration/tests/lal_gate_test.py b/gstlal-calibration/tests/lal_gate_test.py
index 215834ecf1..f1691fe307 100755
--- a/gstlal-calibration/tests/lal_gate_test.py
+++ b/gstlal-calibration/tests/lal_gate_test.py
@@ -47,7 +47,7 @@ def lal_gate_01(pipeline, name):
rate = 512 # Hz
buffer_length = 1.0 # seconds
- test_duration = 100 # seconds
+ test_duration = 16 # seconds
frequency = 0.1 # Hz
attack_length = -3 # seconds
hold_length = 0 # seconds
@@ -60,13 +60,18 @@ def lal_gate_01(pipeline, name):
# Make a sine wave
src = test_common.test_src(pipeline, buffer_length = buffer_length, rate = rate, test_duration = test_duration, wave = 0, volume = 1.0, freq = frequency, width = 64)
+
# Add a DC offset
head = pipeparts.mkgeneric(pipeline, src, "lal_add_constant", value = DC_offset)
head = pipeparts.mktee(pipeline, head)
pipeparts.mknxydumpsink(pipeline, head, "%s_in.txt" % name)
+ control = calibration_parts.mkqueue(pipeline, head, min_length = 8)
+ control = pipeparts.mkgeneric(pipeline, control, "splitcounter", name = "control")
+ head = pipeparts.mkgeneric(pipeline, head, "splitcounter", name = "before")
# Gate it
- head = calibration_parts.mkgate(pipeline, head, head, threshold, attack_length = int(attack_length * rate), hold_length = int(hold_length * rate))
+ head = calibration_parts.mkgate(pipeline, head, control, threshold, attack_length = int(attack_length * rate), hold_length = int(hold_length * rate))
+ head = pipeparts.mkgeneric(pipeline, head, "splitcounter", name = "after")
pipeparts.mknxydumpsink(pipeline, head, "%s_out.txt" % name)
#
diff --git a/gstlal-ugly/gst/lal/gstlal_tdwhiten.c b/gstlal-ugly/gst/lal/gstlal_tdwhiten.c
index 99708e7ba1..a080f1018c 100644
--- a/gstlal-ugly/gst/lal/gstlal_tdwhiten.c
+++ b/gstlal-ugly/gst/lal/gstlal_tdwhiten.c
@@ -19,7 +19,7 @@
/*
* ============================================================================
*
- * Preamble
+ * Preamble
*
* ============================================================================
*/
@@ -32,6 +32,7 @@
#include <string.h>
#include <math.h>
+#include <unistd.h>
/*
@@ -69,7 +70,7 @@
/*
* ============================================================================
*
- * Kernel Info
+ * Kernel Info
*
* ============================================================================
*/
@@ -82,6 +83,12 @@ struct kernelinfo_t {
guint64 offset;
GstClockTime timestamp;
+ /*
+ * timestamp when the kernel is no longer valid
+ */
+
+ GstClockTime endtime;
+
/*
* the kernel
*/
@@ -104,6 +111,7 @@ static struct kernelinfo_t *kernelinfo_new(gsize length)
kernelinfo->offset = GST_BUFFER_OFFSET_NONE;
kernelinfo->timestamp = GST_CLOCK_TIME_NONE;
+ kernelinfo->endtime = GST_CLOCK_TIME_NONE;
kernelinfo->length = length;
kernelinfo->kernel = kernel;
@@ -142,7 +150,7 @@ static gsize kernelinfo_longest(GQueue *kernels)
/*
* ============================================================================
*
- * Utilities
+ * Utilities
*
* ============================================================================
*/
@@ -155,9 +163,9 @@ static unsigned get_input_length(GSTLALTDwhiten *element, unsigned output_length
/*
-* how many output samples can be generated from the given number of input
-* samples
-*/
+ * how many output samples can be generated from the given number of input
+ * samples
+ */
static unsigned get_output_length(const GSTLALTDwhiten *element, unsigned input_length)
@@ -221,7 +229,8 @@ static void set_metadata(GSTLALTDwhiten *element, GstBuffer *buf, guint64 outsam
element->next_out_offset += outsamples;
GST_BUFFER_OFFSET_END(buf) = element->next_out_offset;
GST_BUFFER_PTS(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET(buf) - element->offset0, GST_SECOND, GST_AUDIO_INFO_RATE(&element->audio_info));
- GST_BUFFER_DURATION(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, GST_AUDIO_INFO_RATE(&element->audio_info)) - GST_BUFFER_PTS(buf);
+ element->next_pts = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, GST_AUDIO_INFO_RATE(&element->audio_info));
+ GST_BUFFER_DURATION(buf) = element->next_pts - GST_BUFFER_PTS(buf);
if(element->need_discont) {
GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_DISCONT);
element->need_discont = FALSE;
@@ -441,7 +450,7 @@ static GstFlowReturn filter_and_push(GSTLALTDwhiten *element, guint64 output_len
/*
* ============================================================================
*
- * Signals
+ * Signals
*
* ============================================================================
*/
@@ -466,7 +475,7 @@ static void rate_changed(GstElement *element, gint rate, void *data)
/*
* ============================================================================
*
- * GStreamer Boiler Plate
+ * GStreamer Boiler Plate
*
* ============================================================================
*/
@@ -482,7 +491,7 @@ G_DEFINE_TYPE(
/*
* ============================================================================
*
- * Gst BaseTransform Method Overriddes
+ * Gst BaseTransform Method Overriddes
*
* ============================================================================
*/
@@ -590,9 +599,11 @@ static gboolean start(GstBaseTransform *trans)
GSTLALTDwhiten *element = GSTLAL_TDWHITEN(trans);
element->adapter = g_object_new(GST_TYPE_AUDIOADAPTER, "unit-size", 8, NULL);
element->t0 = GST_CLOCK_TIME_NONE;
+ element->next_pts = 0;
element->offset0 = GST_BUFFER_OFFSET_NONE;
element->next_out_offset = GST_BUFFER_OFFSET_NONE;
element->need_discont = TRUE;
+ g_mutex_init(&element->kernel_lock);
return TRUE;
}
@@ -655,6 +666,7 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
*/
element->t0 = GST_BUFFER_PTS(inbuf);
+ element->next_pts = GST_BUFFER_PTS(inbuf);
element->offset0 = GST_BUFFER_OFFSET(inbuf);
element->next_out_offset = element->offset0 + kernelinfo_longest(element->kernels) - 1 - element->latency;
@@ -667,6 +679,51 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
g_assert_cmpuint(GST_BUFFER_PTS(inbuf), ==, gst_audioadapter_expected_timestamp(element->adapter));
element->next_in_offset = GST_BUFFER_OFFSET_END(inbuf);
+ /*
+ * Check if we are past the end time of the latest kernel in use.
+ * If so, put a waiting kernel into use if available.
+ */
+
+ /*
+ * FIXME: the purpose of the wait loop below is to make the element
+ * wait for a new filter kernel if the current one is "expired."
+ * This should instead be done using the controller:
+ * https://gstreamer.freedesktop.org/documentation/application-development/advanced/dparams.html
+ * https://gstreamer.freedesktop.org/data/doc/gstreamer/head/gstreamer-libs/html/GstTimedValueControlSource.html
+ */
+
+ kernelinfo = g_queue_peek_tail(element->kernels);
+ if(kernelinfo->endtime != GST_CLOCK_TIME_NONE && element->next_pts >= kernelinfo->endtime) {
+ GstClockTime endtime = kernelinfo->endtime;
+ /* Wait until an appropriate filter arrives before continuing */
+ while(element->next_pts >= endtime) {
+ /* First check if we changed our mind about the expiration date fo the current kernel */
+ kernelinfo = g_queue_peek_tail(element->kernels);
+ if(kernelinfo->endtime > element->next_pts)
+ goto setoffsets;
+ kernelinfo = g_queue_peek_head(element->waiting_kernels);
+ if(kernelinfo) {
+ if(kernelinfo->endtime != GST_CLOCK_TIME_NONE) {
+ /* Throw away ones we won't use */
+ if(kernelinfo->endtime <= element->next_pts)
+ kernelinfo_free(g_queue_pop_head(element->waiting_kernels));
+ else
+ endtime = kernelinfo->endtime;
+ }
+ }
+ /* Don't waste CPUs */
+ sleep(1);
+ }
+ /* Dump any kernels we haven't used and won't use */
+ g_mutex_lock(&element->kernel_lock);
+ kernelinfo = g_queue_peek_tail(element->kernels);
+ if(kernelinfo->offset == GST_BUFFER_OFFSET_NONE)
+ kernelinfo_free(g_queue_pop_head(element->kernels));
+ g_queue_push_tail(element->kernels, g_queue_pop_head(element->waiting_kernels));
+ g_mutex_unlock(&element->kernel_lock);
+ }
+setoffsets:
+
/*
* set offsets and timestamps for kernels that are new.
* set_property() will ensure that there is only one new kernel at
@@ -813,7 +870,8 @@ done:
enum property {
PROP_TAPER_LENGTH = 1,
PROP_KERNEL,
- PROP_LATENCY
+ PROP_LATENCY,
+ PROP_KERNEL_ENDTIME
};
@@ -830,16 +888,24 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
case PROP_KERNEL: {
GValueArray *va = g_value_get_boxed(value);
- struct kernelinfo_t *kernelinfo;
+ struct kernelinfo_t *kernelinfo, *old_kernelinfo;
/* dump any kernels that haven't yet been used */
- for(kernelinfo = g_queue_peek_tail(element->kernels); kernelinfo && kernelinfo->offset == GST_BUFFER_OFFSET_NONE; kernelinfo = g_queue_peek_tail(element->kernels))
+ for(kernelinfo = g_queue_peek_tail(element->kernels); kernelinfo && kernelinfo->endtime == GST_CLOCK_TIME_NONE && kernelinfo->offset == GST_BUFFER_OFFSET_NONE; kernelinfo = g_queue_peek_tail(element->kernels))
kernelinfo_free(g_queue_pop_tail(element->kernels));
/* construct new kernelinfo object */
kernelinfo = kernelinfo_new(va->n_values);
gstlal_doubles_from_g_value_array(va, kernelinfo->kernel, NULL);
/* push the new kernel into the queue */
- g_queue_push_tail(element->kernels, kernelinfo);
-
+ old_kernelinfo = g_queue_peek_tail(element->kernels);
+ if(!old_kernelinfo) {
+ if(element->kernel_endtime < G_MAXUINT64)
+ kernelinfo->endtime = element->kernel_endtime;
+ g_queue_push_tail(element->kernels, kernelinfo);
+ } else if(old_kernelinfo->endtime == GST_CLOCK_TIME_NONE)
+ g_queue_push_tail(element->kernels, kernelinfo);
+ else
+ /* this kernel should wait until the end time of the previous one */
+ g_queue_push_tail(element->waiting_kernels, kernelinfo);
break;
}
@@ -849,6 +915,63 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
break;
}
+ case PROP_KERNEL_ENDTIME: {
+ /*
+ * FIXME: This should instead be handled with the controller. Then this property would be unnecessary.
+ * https://gstreamer.freedesktop.org/documentation/application-development/advanced/dparams.html
+ * https://gstreamer.freedesktop.org/data/doc/gstreamer/head/gstreamer-libs/html/GstTimedValueControlSource.html
+ */
+ GstClockTime kernel_endtime = (GstClockTime) g_value_get_uint64(value);
+ if(!g_queue_is_empty(element->kernels) && kernel_endtime < G_MAXUINT64) {
+ g_mutex_lock(&element->kernel_lock);
+ struct kernelinfo_t *kernelinfo;
+ /*
+ * If this endtime is less than the current incoming data timestamps
+ * or less than the end time of the last kernel, we will throw the
+ * kernel away.
+ */
+ GstClockTime min_endtime = element->next_pts;
+ guint i;
+ for(i = 0; i < g_queue_get_length(element->waiting_kernels); i++) {
+ kernelinfo = g_queue_peek_nth(element->waiting_kernels, i);
+ if(kernelinfo->endtime != GST_CLOCK_TIME_NONE)
+ min_endtime = min_endtime > kernelinfo->endtime ? min_endtime : kernelinfo->endtime;
+ }
+ for(i = 0; i < g_queue_get_length(element->kernels); i++) {
+ kernelinfo = g_queue_peek_nth(element->kernels, i);
+ if(kernelinfo->endtime != GST_CLOCK_TIME_NONE)
+ min_endtime = min_endtime > kernelinfo->endtime ? min_endtime : kernelinfo->endtime;
+ }
+ /* First check if there are kernels waiting to be used */
+ if(!g_queue_is_empty(element->waiting_kernels)) {
+ if(kernel_endtime > min_endtime) {
+ kernelinfo = g_queue_peek_tail(element->waiting_kernels);
+ kernelinfo->endtime = kernel_endtime;
+ element->kernel_endtime = kernel_endtime;
+ /* In this case, also throw away any waiting kernels without end times */
+ for(i = 0; i < g_queue_get_length(element->waiting_kernels); i++) {
+ kernelinfo = g_queue_peek_nth(element->waiting_kernels, i);
+ if(kernelinfo->endtime == GST_CLOCK_TIME_NONE)
+ kernelinfo_free(g_queue_pop_nth(element->waiting_kernels, i));
+ }
+ } else
+ kernelinfo_free(g_queue_pop_tail(element->waiting_kernels));
+ } else if(!g_queue_is_empty(element->kernels)) {
+ /* If there are no waiting kernels, apply this to the most recent kernel being used */
+ if(kernel_endtime > min_endtime) {
+ kernelinfo = g_queue_peek_tail(element->kernels);
+ kernelinfo->endtime = kernel_endtime;
+ element->kernel_endtime = kernel_endtime;
+ } else
+ kernelinfo_free(g_queue_pop_tail(element->kernels));
+ }
+ g_mutex_unlock(&element->kernel_lock);
+ } else if (kernel_endtime < G_MAXUINT64)
+ /* The first kernel has not arrived yet, so store it for when it does */
+ element->kernel_endtime = kernel_endtime;
+ break;
+ }
+
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
break;
@@ -888,6 +1011,12 @@ static void get_property(GObject *object, enum property prop_id, GValue *value,
g_value_set_int64(value, element->latency);
break;
+ case PROP_KERNEL_ENDTIME: {
+ /* returns the most recent end time in the queue, or G_MAXUINT64 if no times are present */
+ g_value_set_uint64(value, element->kernel_endtime);
+ break;
+ }
+
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
break;
@@ -908,6 +1037,9 @@ static void finalize(GObject *object)
g_queue_free_full(element->kernels, (GDestroyNotify) kernelinfo_free);
element->kernels = NULL;
+ g_queue_free_full(element->waiting_kernels, (GDestroyNotify) kernelinfo_free);
+ element->waiting_kernels = NULL;
+ g_mutex_clear(&element->kernel_lock);
if(element->adapter) {
g_object_unref(element->adapter);
element->adapter = NULL;
@@ -1019,6 +1151,20 @@ static void gstlal_tdwhiten_class_init(GSTLALTDwhitenClass *klass)
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT | GST_PARAM_CONTROLLABLE
)
);
+ g_object_class_install_property(
+ gobject_class,
+ PROP_KERNEL_ENDTIME,
+ g_param_spec_uint64(
+ "kernel-endtime",
+ "Kernel end time",
+ "Timestamp at which most recently acquired filter kernel becomes invalid.\n\t\t\t"
+ "At this time, element will wait for another filter before processing\n\t\t\t"
+ "more data. This feature can be enabled at any time during data flow,\n\t\t\t"
+ "but should not be disabled once enabled. Default is to leave disabled.",
+ 0, G_MAXUINT64, G_MAXUINT64,
+ G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT | GST_PARAM_CONTROLLABLE
+ )
+ );
signals[SIGNAL_RATE_CHANGED] = g_signal_new(
"rate-changed",
@@ -1049,5 +1195,7 @@ static void gstlal_tdwhiten_init(GSTLALTDwhiten *filter)
filter->adapter = NULL;
filter->latency = 0;
filter->kernels = g_queue_new();
+ filter->waiting_kernels = g_queue_new();
+ filter->kernel_endtime = G_MAXUINT64;
gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(filter), TRUE);
}
diff --git a/gstlal-ugly/gst/lal/gstlal_tdwhiten.h b/gstlal-ugly/gst/lal/gstlal_tdwhiten.h
index 9018eef5cc..d9c684aa77 100644
--- a/gstlal-ugly/gst/lal/gstlal_tdwhiten.h
+++ b/gstlal-ugly/gst/lal/gstlal_tdwhiten.h
@@ -63,13 +63,17 @@ typedef struct {
guint32 taper_length;
GQueue *kernels;
+ GQueue *waiting_kernels;
gint64 latency;
+ guint64 kernel_endtime;
+ GMutex kernel_lock;
/*
* timestamp book-keeping
*/
GstClockTime t0;
+ GstClockTime next_pts;
guint64 offset0;
guint64 next_out_offset;
guint64 next_in_offset;
--
GitLab