From 48dae9a2c5fd333d8ad809a38828d4abd3c223f0 Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Thu, 20 Sep 2018 22:29:44 -0700
Subject: [PATCH] gstlal_compute_strain: fixed hoft_clean state vector bits.
---
gstlal-calibration/bin/gstlal_compute_strain | 4 +-
.../gst/lal/gstlal_transferfunction.c | 40 ++++----
.../gst/lal/gstlal_transferfunction.h | 1 +
.../python/calibration_parts.py | 7 +-
gstlal-calibration/tests/rms_test.py | 95 +++++++++++++++++++
5 files changed, 124 insertions(+), 23 deletions(-)
create mode 100755 gstlal-calibration/tests/rms_test.py
diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain
index 37facaa350..bc04e1fdef 100755
--- a/gstlal-calibration/bin/gstlal_compute_strain
+++ b/gstlal-calibration/bin/gstlal_compute_strain
@@ -2194,8 +2194,8 @@ if remove_cal_lines or remove_power_lines or witness_channel_list is not None:
#
if compute_calib_statevector and (remove_cal_lines or remove_power_lines or witness_channel_list is not None):
- low_rms_rate = pow(2, int(numpy.log(cleaning_check_range_low_max) / numpy.log(2) + 1.1))
- mid_rms_rate = pow(2, int(numpy.log(cleaning_check_range_mid_max) / numpy.log(2) + 1.1))
+ low_rms_rate = pow(2, int(numpy.log(2 * cleaning_check_range_low_max) / numpy.log(2) + 1.01))
+ 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)
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.c b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
index 6b249ba7a5..9a80f146c1 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.c
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
@@ -183,11 +183,11 @@ DEFINE_MAXIMUM(32);
DEFINE_MAXIMUM(64);
-static void write_transfer_functions(complex double *tfs, char *element_name, double df, gint64 rows, int columns, gboolean write_to_screen, char *filename, gboolean free_name) {
+static void write_transfer_functions(complex double *tfs, char *element_name, double df, gint64 rows, int columns, double t_start, double t_finish, gboolean write_to_screen, char *filename, gboolean free_name) {
gint64 i;
int j, j_stop;
if(write_to_screen) {
- g_print("\n\n==================== Transfer functions computed by %s ====================\nfrequency\t\t ", element_name);
+ g_print("\n\n==================== Transfer functions computed by %s from %f until %f ====================\nfrequency\t\t ", element_name, t_start, t_finish);
for(j = 1; j < columns; j++)
g_print("ch%d -> ch0\t\t\t\t ", j);
g_print("ch%d -> ch0\n\n", columns);
@@ -212,7 +212,7 @@ static void write_transfer_functions(complex double *tfs, char *element_name, do
if(filename) {
FILE *fp;
fp = fopen(filename, "a");
- g_fprintf(fp, "==================== Transfer functions computed by %s ====================\nfrequency\t\t ", element_name);
+ g_fprintf(fp, "==================== Transfer functions computed by %s from %f until %f ====================\nfrequency\t\t ", element_name, t_start, t_finish);
for(j = 1; j < columns; j++)
g_fprintf(fp, "ch%d -> ch0\t\t\t\t ", j);
g_fprintf(fp, "ch%d -> ch0\n\n", columns);
@@ -239,11 +239,11 @@ static void write_transfer_functions(complex double *tfs, char *element_name, do
}
-static void write_fir_filters(double *filters, char *element_name, gint64 rows, int columns, gboolean write_to_screen, char *filename, gboolean free_name) {
+static void write_fir_filters(double *filters, char *element_name, gint64 rows, int columns, double t_start, double t_finish, gboolean write_to_screen, char *filename, gboolean free_name) {
gint64 i;
int j, j_stop;
if(write_to_screen) {
- g_print("================== FIR filters computed by %s ==================\n", element_name);
+ g_print("================== FIR filters computed by %s from %f until %f ==================\n", element_name, t_start, t_finish);
for(j = 1; j < columns; j++)
g_print("ch%d -> ch0\t", j);
g_print("ch%d -> ch0\n\n", columns);
@@ -260,7 +260,7 @@ static void write_fir_filters(double *filters, char *element_name, gint64 rows,
if(filename) {
FILE *fp;
fp = fopen(filename, "a");
- g_fprintf(fp, "================== FIR filters computed by %s ==================\n", element_name);
+ g_fprintf(fp, "================== FIR filters computed by %s from %f until %f ==================\n", element_name, t_start, t_finish);
for(j = 1; j < columns; j++)
g_fprintf(fp, "ch%d -> ch0\t", j);
g_fprintf(fp, "ch%d -> ch0\n\n", columns);
@@ -837,14 +837,14 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
} \
} \
} \
- success &= update_transfer_functions_ ## DTYPE(element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix, num_tfs, fd_fft_length, fd_tf_length, element->workspace.w ## S_OR_D ## pf.sinc_table, element->workspace.w ## S_OR_D ## pf.sinc_length, element->workspace.w ## S_OR_D ## pf.sinc_taps_per_df, element->use_median ? 1 : element->num_ffts - element->workspace.w ## S_OR_D ## pf.num_ffts_dropped, element->workspace.w ## S_OR_D ## pf.transfer_functions_at_f, element->workspace.w ## S_OR_D ## pf.transfer_functions_solved_at_f, element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix_at_f, element->workspace.w ## S_OR_D ## pf.permutation, element->transfer_functions); \
+ success &= update_transfer_functions_ ## DTYPE(element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix, num_tfs, fd_fft_length, fd_tf_length, element->workspace.w ## S_OR_D ## pf.sinc_table, element->workspace.w ## S_OR_D ## pf.sinc_length, element->workspace.w ## S_OR_D ## pf.sinc_taps_per_df, element->use_median ? 1 : element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg - element->workspace.w ## S_OR_D ## pf.num_ffts_dropped, element->workspace.w ## S_OR_D ## pf.transfer_functions_at_f, element->workspace.w ## S_OR_D ## pf.transfer_functions_solved_at_f, element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix_at_f, element->workspace.w ## S_OR_D ## pf.permutation, element->transfer_functions); \
if(success) { \
GST_INFO_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 */ \
if(element->write_to_screen || element->filename) \
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->write_to_screen, element->filename, TRUE); \
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg * stride + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE); \
/* If this is this first transfer function after a gap, we may wish to store it */ \
if(element->use_first_after_gap && !element->num_tfs_since_gap) { \
for(i = 0; i < num_tfs * fd_tf_length; i++) \
@@ -871,7 +871,7 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]); \
/* 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->write_to_screen, element->filename, TRUE); \
+ 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) (element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg * stride + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE); \
/* If this is this first FIR filter after a gap, we may wish to store it */ \
if(element->use_first_after_gap && element->num_tfs_since_gap == 1) { \
for(i = 0; i < num_tfs * element->fir_length; i++) \
@@ -1018,7 +1018,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->write_to_screen, element->filename, TRUE);
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (element->workspace.wspf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
} else
GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. No transfer functions will be produced.");
/* Update FIR filters if we want */
@@ -1030,7 +1030,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
/* 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->write_to_screen, element->filename, TRUE);
+ 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) (element->workspace.wspf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
} else
GST_WARNING_OBJECT(element, "FIR filter(s) computation failed. No FIR filters will be produced.");
}
@@ -1093,7 +1093,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->write_to_screen, element->filename, TRUE);
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (element->workspace.wdpf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
} else
GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. No transfer functions will be produced.");
/* Update FIR filters if we want */
@@ -1105,7 +1105,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
/* 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->write_to_screen, element->filename, TRUE);
+ 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) (element->workspace.wdpf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
} else
GST_WARNING_OBJECT(element, "FIR filter(s) computation failed. No FIR filters will be produced.");
}
@@ -1763,7 +1763,7 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), (double) element->rate / element->fir_length, element->fir_length / 2 + 1, element->channels - 1, element->write_to_screen, element->filename, TRUE);
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), (double) element->rate / element->fir_length, element->fir_length / 2 + 1, element->channels - 1, 0, 0, element->write_to_screen, element->filename, TRUE);
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];
@@ -1772,7 +1772,7 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
/* 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, element->channels - 1, element->write_to_screen, element->filename, TRUE);
+ write_fir_filters(element->fir_filters, gst_element_get_name(element), element->fir_length, element->channels - 1, 0, 0, element->write_to_screen, element->filename, TRUE);
}
}
@@ -1812,16 +1812,18 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
gboolean success;
if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
/* If we are just beginning to compute new transfer functions with this data, initialize memory that we will fill to zero */
- if(!element->workspace.wspf.num_ffts_in_avg)
+ if(!element->workspace.wspf.num_ffts_in_avg) {
memset(element->workspace.wspf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * (element->fft_length / 2 + 1) * sizeof(*element->workspace.wspf.autocorrelation_matrix));
-
+ element->t_start_tf = (double) (gst_util_uint64_scale_int_round(GST_BUFFER_PTS(buffer) + GST_BUFFER_DURATION(buffer), element->rate, GST_SECOND) - element->sample_count + element->update_samples) / element->rate;
+ }
/* Send the data to a function to compute fft's and transfer functions */
success = find_transfer_functions_float(element, (float *) mapinfo.data, mapinfo.size, GST_BUFFER_PTS(buffer));
} else {
/* If we are just beginning to compute new transfer functions with this data, initialize memory that we will fill to zero */
- if(!element->workspace.wdpf.num_ffts_in_avg)
+ if(!element->workspace.wdpf.num_ffts_in_avg) {
memset(element->workspace.wdpf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * (element->fft_length / 2 + 1) * sizeof(*element->workspace.wdpf.autocorrelation_matrix));
-
+ element->t_start_tf = (double) (gst_util_uint64_scale_int_round(GST_BUFFER_PTS(buffer) + GST_BUFFER_DURATION(buffer), element->rate, GST_SECOND) - element->sample_count + element->update_samples) / element->rate;
+ }
/* Send the data to a function to compute fft's and transfer functions */
success = find_transfer_functions_double(element, (double *) mapinfo.data, mapinfo.size, GST_BUFFER_PTS(buffer));
}
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.h b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
index acf5aaac4f..1ecec76b18 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.h
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
@@ -78,6 +78,7 @@ struct _GSTLALTransferFunction {
/* Internal state */
gboolean computed_full_tfs;
+ double t_start_tf;
/* transfer function work space */
union {
diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py
index cb60d588ed..e3d1713122 100644
--- a/gstlal-calibration/python/calibration_parts.py
+++ b/gstlal-calibration/python/calibration_parts.py
@@ -379,9 +379,9 @@ def compute_rms(pipeline, head, rate, average_time, f_min = None, f_max = None,
if (f_min is not None) and (f_max is not None):
head = bandpass(pipeline, head, rate, f_low = f_min, f_high = f_max, filter_latency = filter_latency, td = td)
elif f_min is not None:
- head = highpass(pipeline, head, fcut = f_min, filter_latency = filter_latency, td = td)
+ head = highpass(pipeline, head, rate, fcut = f_min, filter_latency = filter_latency, td = td)
elif f_max is not None:
- head = lowpass(pipeline, head, fcut = f_max, filter_latency = filter_latency, td = td)
+ head = lowpass(pipeline, head, rate, fcut = f_max, filter_latency = filter_latency, td = td)
# Square it
head = mkpow(pipeline, head, exponent = 2.0)
@@ -392,6 +392,9 @@ def compute_rms(pipeline, head, rate, average_time, f_min = None, f_max = None,
# Compute running average
head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", default_kappa_re = 0.0, array_size = 1, avg_array_size = average_time * rate_out, filter_latency = filter_latency)
+ # Take the square root
+ head = mkpow(pipeline, head, exponent = 0.5)
+
return head
#
diff --git a/gstlal-calibration/tests/rms_test.py b/gstlal-calibration/tests/rms_test.py
new file mode 100755
index 0000000000..b37259af0e
--- /dev/null
+++ b/gstlal-calibration/tests/rms_test.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+# Copyright (C) 2018 Aaron Viets
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+#
+# =============================================================================
+#
+# Preamble
+#
+# =============================================================================
+#
+
+
+import numpy
+import sys
+from gstlal import pipeparts
+from gstlal import calibration_parts
+from gstlal import test_common
+from gi.repository import Gst
+
+
+#
+# =============================================================================
+#
+# Pipelines
+#
+# =============================================================================
+#
+
+def rms_01(pipeline, name):
+
+ #
+ # This test passes a random series of integers through rms
+ #
+
+ rate_in = 2048 # Hz
+ buffer_length = 1.0 # seconds
+ test_duration = 50.0 # seconds
+
+ #
+ # build pipeline
+ #
+
+ head = test_common.test_src(pipeline, buffer_length = buffer_length, rate = rate_in, width = 64, test_duration = test_duration, wave = 5, freq = 0)
+ white_noise = pipeparts.mktee(pipeline, head)
+ white_noise_over10 = pipeparts.mkaudioamplify(pipeline, white_noise, 0.1)
+ white_noise_bp100to150 = calibration_parts.bandpass(pipeline, white_noise, 2048, length = 1.0, f_low = 100, f_high = 150)
+ white_noise_bp250to300 = calibration_parts.bandpass(pipeline, white_noise, 2048, length = 1.0, f_low = 250, f_high = 300)
+ white_noise_bp450to500 = calibration_parts.bandpass(pipeline, white_noise, 2048, length = 1.0, f_low = 450, f_high = 500)
+ white_noise_bp400to500 = calibration_parts.bandpass(pipeline, white_noise, 2048, length = 1.0, f_low = 400, f_high = 500)
+ rms = calibration_parts.compute_rms(pipeline, white_noise, 1024, 1.0, f_min = 100, f_max = 500)
+ rms_over10 = calibration_parts.compute_rms(pipeline, white_noise_over10, 1024, 1.0, f_min = 100, f_max = 500)
+ rms_bp100to150 = calibration_parts.compute_rms(pipeline, white_noise_bp100to150, 1024, 1.0, f_min = 100, f_max = 500)
+ rms_bp250to300 = calibration_parts.compute_rms(pipeline, white_noise_bp250to300, 1024, 1.0, f_min = 100, f_max = 500)
+ rms_bp450to500 = calibration_parts.compute_rms(pipeline, white_noise_bp450to500, 1024, 1.0, f_min = 100, f_max = 500)
+ rms_bp400to500 = calibration_parts.compute_rms(pipeline, white_noise_bp400to500, 1024, 1.0, f_min = 100, f_max = 500)
+ pipeparts.mknxydumpsink(pipeline, rms, "rms.txt")
+ pipeparts.mknxydumpsink(pipeline, rms_over10, "rms_over10.txt")
+ pipeparts.mknxydumpsink(pipeline, rms_bp100to150, "rms_bp100to150.txt")
+ pipeparts.mknxydumpsink(pipeline, rms_bp250to300, "rms_bp250to300.txt")
+ pipeparts.mknxydumpsink(pipeline, rms_bp450to500, "rms_bp450to500.txt")
+ pipeparts.mknxydumpsink(pipeline, rms_bp400to500, "rms_bp400to500.txt")
+
+ #
+ # done
+ #
+
+ return pipeline
+
+
+#
+# =============================================================================
+#
+# Main
+#
+# =============================================================================
+#
+
+
+test_common.build_and_run(rms_01, "rms_01")
+
+
--
GitLab