From 6ee2014bd4329fe9b2ccb937a6c17f415ec24a01 Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Fri, 9 Aug 2019 15:07:14 -0700
Subject: [PATCH] lal_sensingtdcfs:  New element in gstlal-calibration to solve
 for the time-dependent correction factors of the sensing function

---
 gstlal-calibration/gst/lal/Makefile.am        |   3 +-
 .../gst/lal/gstlal_matrixsolver.h             |   2 +-
 .../gst/lal/gstlal_sensingtdcfs.c             | 876 ++++++++++++++++++
 .../gst/lal/gstlal_sensingtdcfs.h             |  95 ++
 .../gst/lal/gstlalcalibration.c               |   2 +
 .../python/calibration_parts.py               |   2 +-
 .../tests/check_calibration/Makefile          |  17 +-
 .../plot_transfer_function.py                 |  26 +-
 8 files changed, 999 insertions(+), 24 deletions(-)
 create mode 100644 gstlal-calibration/gst/lal/gstlal_sensingtdcfs.c
 create mode 100644 gstlal-calibration/gst/lal/gstlal_sensingtdcfs.h

diff --git a/gstlal-calibration/gst/lal/Makefile.am b/gstlal-calibration/gst/lal/Makefile.am
index 5db3b65871..a605b517c9 100644
--- a/gstlal-calibration/gst/lal/Makefile.am
+++ b/gstlal-calibration/gst/lal/Makefile.am
@@ -20,7 +20,8 @@ lib@GSTPLUGINPREFIX@gstlalcalibration_la_SOURCES = \
 	gstlal_dqtukey.c gstlal_dqtukey.h \
 	gstlal_property.c gstlal_property.h \
 	gstlal_typecast.c gstlal_typecast.h \
-	gstlal_matrixsolver.c gstlal_matrixsolver.h
+	gstlal_matrixsolver.c gstlal_matrixsolver.h \
+	gstlal_sensingtdcfs.c gstlal_sensingtdcfs.h
 lib@GSTPLUGINPREFIX@gstlalcalibration_la_CPPFLAGS = $(AM_CPPFLAGS) $(PYTHON_CPPFLAGS)
 lib@GSTPLUGINPREFIX@gstlalcalibration_la_CFLAGS = $(AM_CFLAGS) $(LAL_CFLAGS) $(GSTLAL_CFLAGS) $(gstreamer_CFLAGS) $(gstreamer_audio_CFLAGS)
 lib@GSTPLUGINPREFIX@gstlalcalibration_la_LDFLAGS = $(AM_LDFLAGS) $(LAL_LIBS) $(GSTLAL_LIBS) $(PYTHON_LIBS) $(gstreamer_LIBS) $(gstreamer_audio_LIBS) $(GSTLAL_PLUGIN_LDFLAGS)
diff --git a/gstlal-calibration/gst/lal/gstlal_matrixsolver.h b/gstlal-calibration/gst/lal/gstlal_matrixsolver.h
index 88b3516aa7..767986d113 100644
--- a/gstlal-calibration/gst/lal/gstlal_matrixsolver.h
+++ b/gstlal-calibration/gst/lal/gstlal_matrixsolver.h
@@ -61,7 +61,7 @@ struct _GSTLALMatrixSolver {
 		GSTLAL_MATRIXSOLVER_F32 = 0,
 		GSTLAL_MATRIXSOLVER_F64,
 		GSTLAL_MATRIXSOLVER_Z64,
-		GSTLAL_MATRIXSOLVER_Z128,
+		GSTLAL_MATRIXSOLVER_Z128
 	} data_type;
 
 	/* timestamp book-keeping */
diff --git a/gstlal-calibration/gst/lal/gstlal_sensingtdcfs.c b/gstlal-calibration/gst/lal/gstlal_sensingtdcfs.c
new file mode 100644
index 0000000000..a16d08122c
--- /dev/null
+++ b/gstlal-calibration/gst/lal/gstlal_sensingtdcfs.c
@@ -0,0 +1,876 @@
+/*
+ * Copyright (C) 2019  Aaron Viets <aaron.viets@ligo.org>
+ *
+ * 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
+ *
+ * ============================================================================
+ */
+
+
+/*
+ * stuff from C
+ */
+
+
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+
+/*
+ *  stuff from gobject/gstreamer
+ */
+
+
+#include <glib.h>
+#include <gst/gst.h>
+#include <gst/audio/audio.h>
+#include <gst/base/gstbasetransform.h>
+
+
+/*
+ * our own stuff
+ */
+
+
+#include <gstlal/gstlal_audio_info.h>
+#include <gstlal_sensingtdcfs.h>
+
+
+/*
+ * ============================================================================
+ *
+ *			   GStreamer Boiler Plate
+ *
+ * ============================================================================
+ */
+
+
+#define INCAPS \
+	"audio/x-raw, " \
+	"format = (string) {"GST_AUDIO_NE(Z64)", "GST_AUDIO_NE(Z128)"}, " \
+	"rate = " GST_AUDIO_RATE_RANGE ", " \
+	"channels = (int) [4, 5], " \
+	"layout = (string) interleaved, " \
+	"channel-mask = (bitmask) 0"
+
+#define OUTCAPS \
+	"audio/x-raw, " \
+	"format = (string) {"GST_AUDIO_NE(F32)", "GST_AUDIO_NE(F64)"}, " \
+	"rate = " GST_AUDIO_RATE_RANGE ", " \
+	"channels = (int) [4, 6], " \
+	"layout = (string) interleaved, " \
+	"channel-mask = (bitmask) 0"
+
+
+static GstStaticPadTemplate sink_factory = GST_STATIC_PAD_TEMPLATE(
+	GST_BASE_TRANSFORM_SINK_NAME,
+	GST_PAD_SINK,
+	GST_PAD_ALWAYS,
+	GST_STATIC_CAPS(INCAPS)
+);
+
+
+static GstStaticPadTemplate src_factory = GST_STATIC_PAD_TEMPLATE(
+	GST_BASE_TRANSFORM_SRC_NAME,
+	GST_PAD_SRC,
+	GST_PAD_ALWAYS,
+	GST_STATIC_CAPS(OUTCAPS)
+);
+
+
+G_DEFINE_TYPE(
+	GSTLALSensingTDCFs,
+	gstlal_sensingtdcfs,
+	GST_TYPE_BASE_TRANSFORM
+);
+
+
+/*
+ * ============================================================================
+ *
+ *				Utilities
+ *
+ * ============================================================================
+ */
+
+
+/*
+ * set the metadata on an output buffer
+ */
+
+
+static void set_metadata(GSTLALSensingTDCFs *element, GstBuffer *buf, guint64 outsamples, gboolean gap) {
+
+	GST_BUFFER_OFFSET(buf) = element->next_out_offset;
+	element->next_out_offset += outsamples;
+	GST_BUFFER_OFFSET_END(buf) = element->next_out_offset;
+	GST_BUFFER_TIMESTAMP(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET(buf) - element->offset0, GST_SECOND, element->rate);
+	GST_BUFFER_DURATION(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, element->rate) - GST_BUFFER_TIMESTAMP(buf);
+	GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP);
+	if(G_UNLIKELY(element->need_discont)) {
+		GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_DISCONT);
+		element->need_discont = FALSE;
+	}
+	if(gap)
+		GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_GAP);
+	else
+		GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP);
+}
+
+
+/*
+ * Macro functions to solve for kappa_C using three different sensing function models: 0, 1, and 2
+ */
+
+
+#define DEFINE_KAPPA_C_0(DTYPE, F_OR_NOT) \
+DTYPE kappa_C_0_ ## DTYPE(DTYPE Gres1_real, DTYPE Gres1_imag, DTYPE Gres2_real, DTYPE Gres2_imag, DTYPE Y1_real, DTYPE Y1_imag, DTYPE Y2_real, DTYPE Y2_imag, DTYPE f1, DTYPE f2) { \
+ \
+	DTYPE H1, H2, I1, I2, Xi, Zeta, a, b, c, d, e, Delta0, Delta1, p, q, R0; \
+	complex DTYPE Q0, S0; \
+	H1 = f1 * f1 * (Gres1_real - Y1_real); \
+	H2 = f2 * f2 * (Gres2_real - Y2_real); \
+	I1 = f1 * f1 * (Gres1_imag - Y1_imag); \
+	I2 = f2 * f2 * (Gres2_imag - Y2_imag); \
+	Xi = (f2 * f2 * H1 - f1 * f1 * H2) / (f1 * f1 - f2 * f2); \
+	Zeta = I1 / f1 - I2 / f2; \
+	a = f2 * f2 * Xi * (H2 + Xi) * Zeta * Zeta; \
+	b = f2 * (f2 * f2 - f1 * f1) * (H2 + Xi) * Zeta * I2 + pow ## F_OR_NOT(f2, 4) * (H2 + 2 * Xi) * Zeta * Zeta; \
+	c = pow ## F_OR_NOT(f2, 3) * (f2 * f2 - f1 * f1) * Zeta * I2 + pow ## F_OR_NOT(f2 * f2 - f1 * f1, 2) * pow ## F_OR_NOT(H2 + Xi, 2) + pow ## F_OR_NOT(f2, 6) * Zeta * Zeta; \
+	d = 2 * f2 * f2 * pow ## F_OR_NOT(f2 * f2 - f1 * f1, 2) * (H2 + Xi); \
+	e = pow ## F_OR_NOT(f2, 4) * pow ## F_OR_NOT(f2 * f2 - f1 * f1, 2); \
+	Delta0 = c * c - 3 * b * d + 12 * a * e; \
+	Delta1 = 2 * pow ## F_OR_NOT(c, 3) - 9 * b * c * d + 27 * b * b * e + 27 * a * d * d - 72 * a * c * e; \
+	p = (8 * a * c - 3 * b * b) / (8 * a * a); \
+	q = (pow ## F_OR_NOT(b, 3) - 4 * a * b * c + 8 * a * a * d) / (8 * pow ## F_OR_NOT(a, 3)); \
+	Q0 = cpow ## F_OR_NOT(0.5 * ((complex DTYPE) Delta1 + cpow ## F_OR_NOT((complex DTYPE) (Delta1 * Delta1 - 4 * pow ## F_OR_NOT(Delta0, 3)), 0.5)), 1.0 / 3.0); \
+	S0 = 0.5 * cpow ## F_OR_NOT((-2 * p) / 3.0 + (Q0 + Delta0 / Q0) / (3 * a), 0.5); \
+	R0 = pow ## F_OR_NOT(b, 3) + 8 * d * a * a - 4 * a * b * c; \
+	if(R0 <= 0.0) \
+		return creal ## F_OR_NOT(-b / (4 * a) - S0 + 0.5 * cpow ## F_OR_NOT(-4 * S0 * S0 - 2 * p + q / S0, 0.5)); \
+	else \
+		return creal ## F_OR_NOT(-b / (4 * a) + S0 + 0.5 * cpow ## F_OR_NOT(-4 * S0 * S0 - 2 * p - q / S0, 0.5)); \
+}
+
+
+DEFINE_KAPPA_C_0(float, f);
+DEFINE_KAPPA_C_0(double, );
+
+
+/*
+ * Macro functions to solve for f_cc using three different sensing function models: 0, 1, and 2
+ */
+
+
+#define DEFINE_F_CC_0(DTYPE) \
+DTYPE f_cc_0_ ## DTYPE(DTYPE Gres1_imag, DTYPE Gres2_imag, DTYPE Y1_imag, DTYPE Y2_imag, DTYPE f1, DTYPE f2, DTYPE kappa_C) { \
+ \
+	return (f2 * f2 - f1 * f1) / (kappa_C * (f1 * (Gres1_imag - Y1_imag) - f2 * (Gres2_imag - Y2_imag))); \
+}
+
+
+DEFINE_F_CC_0(float);
+DEFINE_F_CC_0(double);
+
+
+/*
+ * Macro functions to solve for f_s^2 using three different sensing function models: 0, 1, and 2
+ */
+
+
+#define DEFINE_F_S_SQUARED_0(DTYPE) \
+DTYPE f_s_squared_0_ ## DTYPE(DTYPE Gres1_real, DTYPE Gres2_real, DTYPE Y1_real, DTYPE Y2_real, DTYPE f1, DTYPE f2, DTYPE kappa_C) { \
+ \
+	return kappa_C * (Y2_real - Gres2_real - Y1_real + Gres1_real) / (1.0 / (f2 * f2) - 1.0 / (f1 * f1)); \
+}
+
+
+DEFINE_F_S_SQUARED_0(float);
+DEFINE_F_S_SQUARED_0(double);
+
+
+/*
+ * Macro functions to solve for f_s / Q using three different sensing function models: 0, 1, and 2
+ */
+
+
+#define DEFINE_FS_OVER_Q_0(DTYPE) \
+DTYPE f_s_over_Q_0_ ## DTYPE(DTYPE Gres1_real, DTYPE Y1_real, DTYPE f1, DTYPE kappa_C, DTYPE f_cc, DTYPE f_s_squared) { \
+ \
+	return f_cc * (kappa_C * Y1_real - 1.0 - f_s_squared / (f1 * f1) - kappa_C * Gres1_real); \
+}
+
+
+DEFINE_FS_OVER_Q_0(float);
+DEFINE_FS_OVER_Q_0(double);
+
+
+/*
+ * ============================================================================
+ *
+ *		     GstBaseTransform Method Overrides
+ *
+ * ============================================================================
+ */
+
+
+/*
+ * get_unit_size()
+ */
+
+
+static gboolean get_unit_size(GstBaseTransform *trans, GstCaps *caps, gsize *size)
+{
+	GstAudioInfo info;
+	gboolean success = TRUE;
+
+	success &= gstlal_audio_info_from_caps(&info, caps);
+
+	if(success) {
+		*size = GST_AUDIO_INFO_BPF(&info);
+	} else
+		GST_WARNING_OBJECT(trans, "unable to parse caps %" GST_PTR_FORMAT, caps);
+
+	return success;
+}
+
+
+/*
+ * transform_caps()
+ */
+
+
+static GstCaps *transform_caps(GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps, GstCaps *filter) {
+
+	guint n;
+	caps = gst_caps_normalize(gst_caps_copy(caps));
+
+	switch(direction) {
+	case GST_PAD_SRC:
+		/* 
+		 * We know the caps on the source pad, and we want to put constraints on
+		 * the sink pad caps.  The sink pad caps are complex, while the source pad
+		 * caps are real.  If there are 4 channels on the source pad, there should
+		 * be 4 on the sink pad as well.  If there are 6 on the source pad, there
+		 * should be 5 on the sink pad.  There are no other possibilities.
+		 */
+		for(n = 0; n < gst_caps_get_size(caps); n++) {
+			GstStructure *str = gst_caps_get_structure(caps, n);
+			const GValue *v = gst_structure_get_value(str, "channels");
+			if(GST_VALUE_HOLDS_INT_RANGE(v)) {
+				gint channels_out_min, channels_out_max;
+				channels_out_min = gst_value_get_int_range_min(v);
+				channels_out_max = gst_value_get_int_range_max(v);
+				g_assert_cmpint(channels_out_min, <=, 6);
+				g_assert_cmpint(channels_out_max, >=, 4);
+				if(channels_out_min >= 5)
+					/* Then we know there must be 5 on the sink pad */
+					gst_structure_set(str, "channels", G_TYPE_INT, 5, NULL);
+				else if(channels_out_max <=5)
+					/* Then we know there must be 4 on the sink pad */
+					gst_structure_set(str, "channels", G_TYPE_INT, 4, NULL);
+				else
+					gst_structure_set(str, "channels", GST_TYPE_INT_RANGE, 4, 5, NULL);
+
+			} else if(G_VALUE_HOLDS_INT(v)) {
+				gint channels_out = g_value_get_int(v);
+				if(channels_out == 4)
+					gst_structure_set(str, "channels", G_TYPE_INT, 4, NULL);
+				else if(channels_out == 6)
+					gst_structure_set(str, "channels", G_TYPE_INT, 5, NULL);
+				else
+					GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid number of channels in caps"));
+
+			} else
+				GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid type for channels in caps"));
+
+			const gchar *format = gst_structure_get_string(str, "format");
+			if(!format) {
+				GST_DEBUG_OBJECT(trans, "unrecognized caps %" GST_PTR_FORMAT, caps);
+				goto error;
+			} else if(!strcmp(format, GST_AUDIO_NE(F32)))
+				gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(Z64), NULL);
+			else if(!strcmp(format, GST_AUDIO_NE(F64)))
+				gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(Z128), NULL);
+			else {
+				GST_DEBUG_OBJECT(trans, "unrecognized format %s in %" GST_PTR_FORMAT, format, caps);
+				goto error;
+			}
+		}
+		break;
+
+	case GST_PAD_SINK:
+		/*
+		 * We know the caps on the sink pad, and we want to put constraints on
+		 * the source pad caps.  The sink pad caps are complex, while the source pad
+		 * caps are real.  If there are 4 channels on the sink pad, there should
+		 * be 4 on the source pad as well.  If there are 5 on the sink pad, there
+		 * should be 6 on the source pad.  There are no other possibilities.
+		 */
+		for(n = 0; n < gst_caps_get_size(caps); n++) {
+			GstStructure *str = gst_caps_get_structure(caps, n);
+			const GValue *v = gst_structure_get_value(str, "channels");
+			if(GST_VALUE_HOLDS_INT_RANGE(v)) {
+				gint channels_in_min, channels_in_max;
+				channels_in_min = gst_value_get_int_range_min(v);
+				channels_in_max = gst_value_get_int_range_max(v);
+				g_assert_cmpint(channels_in_min, <=, 5);
+				g_assert_cmpint(channels_in_max, >=, 4);
+				if(channels_in_min == 5)
+					/* Then we know there must be 6 channels on the source pad */
+					gst_structure_set(str, "channels", G_TYPE_INT, 6, NULL);
+				else if(channels_in_max == 4)
+					/* Then we know there must be 4 channels on the source pad */
+					gst_structure_set(str, "channels", G_TYPE_INT, 4, NULL);
+				else
+					gst_structure_set(str, "channels", GST_TYPE_INT_RANGE, 4, 6, NULL);
+
+			} else if(G_VALUE_HOLDS_INT(v)) {
+				gint channels_in;
+				channels_in = g_value_get_int(v);
+				if(channels_in == 4)
+					gst_structure_set(str, "channels", G_TYPE_INT, 4, NULL);
+				else if(channels_in == 5)
+					gst_structure_set(str, "channels", G_TYPE_INT, 6, NULL);
+				else
+					GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid number of channels in caps"));
+
+			} else
+				GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid type for channels in caps"));
+
+			const gchar *format = gst_structure_get_string(str, "format");
+			if(!format) {
+				GST_DEBUG_OBJECT(trans, "unrecognized caps %" GST_PTR_FORMAT, caps);
+				goto error;
+			} else if(!strcmp(format, GST_AUDIO_NE(Z64)))
+				gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(F32), NULL);
+			else if(!strcmp(format, GST_AUDIO_NE(Z128)))
+				gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(F64), NULL);
+			else {
+				GST_DEBUG_OBJECT(trans, "unrecognized format %s in %" GST_PTR_FORMAT, format, caps);
+				goto error;
+			}
+		}
+		break;
+
+	case GST_PAD_UNKNOWN:
+		GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid direction GST_PAD_UNKNOWN"));
+		goto error;
+	default:
+		g_assert_not_reached();
+	}
+
+	if(filter) {
+		GstCaps *intersection = gst_caps_intersect(caps, filter);
+		gst_caps_unref(caps);
+		caps = intersection;
+	}
+	return gst_caps_simplify(caps);
+
+error:
+	gst_caps_unref(caps);
+	return GST_CAPS_NONE;
+}
+
+
+/*
+ * set_caps()
+ */
+
+
+static gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, GstCaps *outcaps)
+{
+	GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(trans);
+	gint rate, channels_in, channels_out;
+	gsize unit_size_in, unit_size_out;
+
+	/*
+ 	 * parse the caps
+ 	 */
+
+	GstStructure *outstr = gst_caps_get_structure(outcaps, 0);
+	GstStructure *instr = gst_caps_get_structure(incaps, 0);
+	const gchar *name = gst_structure_get_string(outstr, "format");
+	if(!name) {
+		GST_DEBUG_OBJECT(element, "unable to parse format from %" GST_PTR_FORMAT, outcaps);
+		return FALSE;
+	}
+	if(!get_unit_size(trans, outcaps, &unit_size_out)) {
+		GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed");
+		return FALSE;
+	}
+	if(!get_unit_size(trans, incaps, &unit_size_in)) {
+		GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed");
+		return FALSE;
+	}
+	if(!gst_structure_get_int(outstr, "rate", &rate)) {
+		GST_DEBUG_OBJECT(element, "unable to parse rate from %" GST_PTR_FORMAT, outcaps);
+		return FALSE;
+	}
+	if(!gst_structure_get_int(outstr, "channels", &channels_out)) {
+		GST_DEBUG_OBJECT(element, "unable to parse channels from %" GST_PTR_FORMAT, outcaps);
+		return FALSE;
+	}
+	if(!gst_structure_get_int(instr, "channels", &channels_in)) {
+		GST_DEBUG_OBJECT(element, "unable to parse channels from %" GST_PTR_FORMAT, incaps);
+		return FALSE;
+	}
+
+	/* Requirements for channels */
+	if(channels_in == 4) {
+		g_assert_cmpint(channels_out, ==, 4);
+		if(element->sensing_model != 0 && element->sensing_model != 1) {
+			GST_WARNING_OBJECT(element, "When there are 4 input channels, sensing-model must be either 0 or 1.  Resetting sensing-model to 0.");
+			element->sensing_model = 0;
+		}
+	} else if(channels_in == 5) {
+		g_assert_cmpint(channels_out, ==, 6);
+		if(element->sensing_model != 2) {
+			GST_WARNING_OBJECT(element, "When there are 5 input channels, sensing-model must be 2.  Resetting sensing-model to 2.");
+			element->sensing_model = 2;
+		}
+	} else
+		g_assert_not_reached();
+
+	/*
+ 	 * record stream parameters
+ 	 */
+
+	if(!strcmp(name, GST_AUDIO_NE(F32))) {
+		element->data_type = GSTLAL_SENSINGTDCFS_FLOAT;
+		g_assert_cmpuint(unit_size_out, ==, 4 * (guint) channels_out);
+		g_assert_cmpuint(unit_size_in, ==, 4 * (guint) channels_in);
+	} else if(!strcmp(name, GST_AUDIO_NE(F64))) {
+		element->data_type = GSTLAL_SENSINGTDCFS_DOUBLE;
+		g_assert_cmpuint(unit_size_out, ==, 8 * (guint) channels_out);
+		g_assert_cmpuint(unit_size_in, ==, 8 * (guint) channels_in);
+	} else
+		g_assert_not_reached();
+
+	element->rate = rate;
+	element->channels_out = channels_out;
+	element->channels_in = channels_in;
+	element->unit_size_out = unit_size_out;
+	element->unit_size_in = unit_size_in;
+
+	return TRUE;
+}
+
+
+/*
+ * transform_size{}
+ */
+
+
+static gboolean transform_size(GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps, gsize size, GstCaps *othercaps, gsize *othersize) {
+
+	GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(trans);
+
+	/*
+	 * The data types of inputs and outputs are the same, but the number of channels differs.
+	 * For N output channels, there are N(N+1) input channels.
+	 */
+
+	switch(direction) {
+	case GST_PAD_SRC:
+		/*
+		 * We know the size of the output buffer and want to compute the size of the input buffer.
+		 * The size of the output buffer should be a multiple of the unit_size.
+		 */
+
+		if(G_UNLIKELY(size % element->unit_size_out)) {
+			GST_DEBUG_OBJECT(element, "buffer size %" G_GSIZE_FORMAT " is not a multiple of %" G_GSIZE_FORMAT, size, (gsize) element->unit_size_out);
+			return FALSE;
+		}
+
+		*othersize = size * element->unit_size_in / element->unit_size_out;
+
+		break;
+
+	case GST_PAD_SINK:
+		/*
+		 * We know the size of the input buffer and want to compute the size of the output buffer.
+		 * The size of the output buffer should be a multiple of unit_size * (N+1).
+		 */
+
+		if(G_UNLIKELY(size % element->unit_size_in)) {
+			GST_ERROR_OBJECT(element, "buffer size %" G_GSIZE_FORMAT " is not a multiple of %" G_GSIZE_FORMAT, size, (gsize) element->unit_size_in);
+			return FALSE;
+		}
+
+		*othersize = size * element->unit_size_out / element->unit_size_in;
+
+		break;
+
+	case GST_PAD_UNKNOWN:
+		GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid direction GST_PAD_UNKNOWN"));
+		return FALSE;
+	}
+
+	return TRUE;
+}
+
+
+/*
+ * start()
+ */
+
+
+static gboolean start(GstBaseTransform *trans) {
+
+	GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(trans);
+
+	element->t0 = GST_CLOCK_TIME_NONE;
+	element->offset0 = GST_BUFFER_OFFSET_NONE;
+	element->next_in_offset = GST_BUFFER_OFFSET_NONE;
+	element->next_out_offset = GST_BUFFER_OFFSET_NONE;
+	element->need_discont = TRUE;
+
+	/* Sanity checks */
+	if(element->freq1 == G_MAXDOUBLE)
+		GST_WARNING_OBJECT(element, "freq1 was not set. It must be set in order to produce sensible output.");
+	if(element->freq2 == G_MAXDOUBLE)
+		GST_WARNING_OBJECT(element, "freq2 was not set. It must be set in order to produce sensible output.");
+	if(element->sensing_model == 2 && element->freq2 == G_MAXDOUBLE)
+		GST_WARNING_OBJECT(element, "freq4 was not set. When using sensing-model is 2, freq4 must be set in order to produce sensible output.");
+
+	return TRUE;
+}
+
+
+/*
+ * transform()
+ */
+
+
+static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuffer *outbuf) {
+
+	GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(trans);
+	GstMapInfo inmap, outmap;
+	GstFlowReturn result = GST_FLOW_OK;
+
+	/*
+	 * check for discontinuity
+	 */
+
+	if(G_UNLIKELY(GST_BUFFER_IS_DISCONT(inbuf) || GST_BUFFER_OFFSET(inbuf) != element->next_in_offset || !GST_CLOCK_TIME_IS_VALID(element->t0))) {
+		GST_DEBUG_OBJECT(element, "pushing discontinuous buffer at input timestamp %lu", (long unsigned) GST_TIME_AS_SECONDS(GST_BUFFER_PTS(inbuf)));
+		element->t0 = GST_BUFFER_PTS(inbuf);
+		element->offset0 = element->next_out_offset = GST_BUFFER_OFFSET(inbuf);
+		element->need_discont = TRUE;
+	}
+	element->next_in_offset = GST_BUFFER_OFFSET_END(inbuf);
+
+	/*
+	 * process buffer
+	 */
+
+	if(!GST_BUFFER_FLAG_IS_SET(inbuf, GST_BUFFER_FLAG_GAP)) {
+
+		/*
+		 * input is not gap.
+		 */
+
+		gst_buffer_map(inbuf, &inmap, GST_MAP_READ);
+		gst_buffer_map(outbuf, &outmap, GST_MAP_WRITE);
+		if(element->data_type == GSTLAL_SENSINGTDCFS_FLOAT) {
+			complex float *indata = (complex float *) inmap.data;
+			float *outdata = (float *) outmap.data;
+			guint64 i, samples = outmap.size / element->unit_size_out;
+			float f1 = (float) element->freq1;
+			float f2 = (float) element->freq2;
+			float kappa_C, f_cc;
+			switch(element->sensing_model) {
+			case 0: ;
+				float f_s_squared, f_s_over_Q;
+				for(i = 0; i < samples; i++) {
+					kappa_C = kappa_C_0_float(crealf(indata[4 * i]), cimagf(indata[4 * i]), crealf(indata[4 * i + 1]), cimagf(indata[4 * i + 1]), crealf(indata[4 * i + 2]), cimagf(indata[4 * i + 2]), crealf(indata[4 * i + 3]), cimagf(indata[4 * i + 3]), f1, f2);
+					f_cc = f_cc_0_float(cimagf(indata[4 * i]), cimagf(indata[4 * i + 1]), cimagf(indata[4 * i + 2]), cimagf(indata[4 * i + 3]), f1, f2, kappa_C);
+					f_s_squared = f_s_squared_0_float(crealf(indata[4 * i]), crealf(indata[4 * i + 1]), crealf(indata[4 * i + 2]), crealf(indata[4 * i + 3]), f1, f2, kappa_C);
+					f_s_over_Q = f_s_over_Q_0_float(crealf(indata[4 * i]), crealf(indata[4 * i + 2]), f1, kappa_C, f_cc, f_s_squared);
+					outdata[4 * i] = kappa_C;
+					outdata[4 * i + 1] = f_cc;
+					outdata[4 * i + 2] = f_s_squared;
+					outdata[4 * i + 3] = f_s_over_Q;
+				}
+				break;
+			case 1:
+				/* Solution not developed */
+				break;
+			case 2:
+				/* Solution not developed */
+				break;
+			default:
+				g_assert_not_reached();
+			}
+		} else if (element->data_type == GSTLAL_SENSINGTDCFS_DOUBLE) {
+			complex double *indata = (complex double *) inmap.data;
+			double *outdata = (double *) outmap.data;
+			guint64 i, samples = outmap.size / element->unit_size_out;
+			double f1 = element->freq1;
+			double f2 = element->freq2;
+			double kappa_C, f_cc;
+			switch(element->sensing_model) {
+			case 0: ;
+				double f_s_squared, f_s_over_Q;
+				for(i = 0; i < samples; i++) {
+					kappa_C = kappa_C_0_double(crealf(indata[4 * i]), cimagf(indata[4 * i]), crealf(indata[4 * i + 1]), cimagf(indata[4 * i + 1]), crealf(indata[4 * i + 2]), cimagf(indata[4 * i + 2]), crealf(indata[4 * i + 3]), cimagf(indata[4 * i + 3]), f1, f2);
+					f_cc = f_cc_0_double(cimagf(indata[4 * i]), cimagf(indata[4 * i + 1]), cimagf(indata[4 * i + 2]), cimagf(indata[4 * i + 3]), f1, f2, kappa_C);
+					f_s_squared = f_s_squared_0_double(crealf(indata[4 * i]), crealf(indata[4 * i + 1]), crealf(indata[4 * i + 2]), crealf(indata[4 * i + 3]), f1, f2, kappa_C);
+					f_s_over_Q = f_s_over_Q_0_double(crealf(indata[4 * i]), crealf(indata[4 * i + 2]), f1, kappa_C, f_cc, f_s_squared);
+					outdata[4 * i] = kappa_C;
+					outdata[4 * i + 1] = f_cc;
+					outdata[4 * i + 2] = f_s_squared;
+					outdata[4 * i + 3] = f_s_over_Q;
+				}
+				break;
+			case 1:
+				/* Solution not developed */
+				break;
+			case 2:
+				/* Solution not developed */
+				break;
+			default:
+				g_assert_not_reached();
+			}
+		} else {
+			g_assert_not_reached();
+		}
+		set_metadata(element, outbuf, outmap.size / element->unit_size_out, FALSE);
+		gst_buffer_unmap(outbuf, &outmap);
+		gst_buffer_unmap(inbuf, &inmap);
+
+	} else {
+
+		/*
+		 * input is gap.
+		 */
+
+		gst_buffer_map(outbuf, &outmap, GST_MAP_WRITE);
+		memset(outmap.data, 0, outmap.size);
+		set_metadata(element, outbuf, outmap.size / element->unit_size_out, TRUE);
+		gst_buffer_unmap(outbuf, &outmap);
+	}
+
+	/*
+	 * done
+	 */
+
+	return result;
+}
+
+
+/*
+ * ============================================================================
+ *
+ *			  GObject Method Overrides
+ *
+ * ============================================================================
+ */
+
+
+/*
+ * properties
+ */
+
+
+enum property {
+	ARG_SENSING_MODEL = 1,
+	ARG_FREQ1,
+	ARG_FREQ2,
+	ARG_FREQ4
+};
+
+
+static void set_property(GObject *object, enum property prop_id, const GValue *value, GParamSpec *pspec) {
+
+	GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(object);
+
+	GST_OBJECT_LOCK(element);
+
+	switch (prop_id) {
+	case ARG_SENSING_MODEL:
+		element->sensing_model = g_value_get_int(value);
+		break;
+	case ARG_FREQ1:
+		element->freq1 = g_value_get_double(value);
+		break;
+	case ARG_FREQ2:
+		element->freq2 = g_value_get_double(value);
+		break;
+	case ARG_FREQ4:
+		element->freq4 = g_value_get_double(value);
+		break;
+	default:
+		G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
+		break;
+	}
+
+	GST_OBJECT_UNLOCK(element);	
+}
+
+
+static void get_property(GObject *object, enum property prop_id, GValue *value, GParamSpec *pspec) {
+
+	GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(object);
+
+	GST_OBJECT_LOCK(element);
+
+	switch (prop_id) {
+	case ARG_SENSING_MODEL:
+		g_value_set_int(value, element->sensing_model);
+		break;
+	case ARG_FREQ1:
+		g_value_set_double(value, element->freq1);
+		break;
+	case ARG_FREQ2:
+		g_value_set_double(value, element->freq2);
+		break;
+	case ARG_FREQ4:
+		g_value_set_double(value, element->freq4);
+		break;
+	default:
+		G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
+		break;
+	}
+
+	GST_OBJECT_UNLOCK(element);
+}
+
+
+/*
+ * class_init()
+ */
+
+
+static void gstlal_sensingtdcfs_class_init(GSTLALSensingTDCFsClass *klass)
+{
+	GstBaseTransformClass *transform_class = GST_BASE_TRANSFORM_CLASS(klass);
+	GstElementClass *element_class = GST_ELEMENT_CLASS(klass);
+	GObjectClass *gobject_class = G_OBJECT_CLASS(klass);
+
+	gst_element_class_set_details_simple(
+		element_class,
+		"Sensing TDCFs",
+		"Filter/Audio",
+		"Solves for the time-dependent correction factors of the sensing function using\n\t\t\t   "
+		"the solution described in LIGO DCC document P1900052, Section 5.2.6.  It takes\n\t\t\t   "
+		"the complex inputs Gres^1, Gres^2, Y^1, and Y^2 (and Y^3 if sensing-model is 2),\n\t\t\t   "
+		"in that order.  The outputs are kappa_C, f_cc, f_s, and Q, in that order.\n\t\t\t   "
+		"Currently, sensing-model=0 is the only sensing model that has been developed.",
+		"Aaron Viets <aaron.viets@ligo.org>"
+	);
+
+	gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&src_factory));
+	gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&sink_factory));
+
+	transform_class->get_unit_size = GST_DEBUG_FUNCPTR(get_unit_size);
+	transform_class->set_caps = GST_DEBUG_FUNCPTR(set_caps);
+	transform_class->transform_caps = GST_DEBUG_FUNCPTR(transform_caps);
+	transform_class->transform_size = GST_DEBUG_FUNCPTR(transform_size);
+	transform_class->start = GST_DEBUG_FUNCPTR(start);
+	transform_class->transform = GST_DEBUG_FUNCPTR(transform);
+
+	gobject_class->set_property = GST_DEBUG_FUNCPTR(set_property);
+	gobject_class->get_property = GST_DEBUG_FUNCPTR(get_property);
+
+	g_object_class_install_property(
+		gobject_class,
+		ARG_SENSING_MODEL,
+		g_param_spec_int(
+			"sensing-model",
+			"Sensing Model",
+			"Which model of the sensing function to use.  Each model is described below:\n\n\t\t\t"
+			"sensing-model=0:\n\t\t\t"
+			"    kappa_C\t\t    f^2\n\t\t\t"
+			" -------------- * ------------------------- * C_res(f)\n\t\t\t"
+			" 1 + i f / f_cc   f^2 + f_s^2 - i f f_s / Q\n\n\t\t\t"
+			"Complex inputs: Gres1, Gres2, Y1, Y2.  See P1900052, Sec. 5.2.6.\n\t\t\t"
+			"Real outputs: kappa_C, f_cc, f_s^2, f_s / Q\n\n\t\t\t"
+			"sensing-model=1:\n\t\t\t"
+			"    kappa_C\n\t\t\t"
+			" -------------- * (\?\?\?)... not yet modeled\n\t\t\t"
+			" 1 + i f / f_cc\n\n\t\t\t"
+			"Complex inputs: Gres1, Gres2, Y1, Y2.\n\t\t\t"
+			"Real outputs: kappa_C, f_cc, \?\?, \?\?\n\n\t\t\t"
+			"sensing-model=2:\n\t\t\t"
+			"    kappa_C\t\t    f^2\n\t\t\t"
+			" -------------- * ------------------------- * (\?\?\?)... not yet modeled.\n\t\t\t"
+			" 1 + i f / f_cc   f^2 + f_s^2 - i f f_s / Q\n\n\t\t\t"
+			"Complex inputs: Gres1, Gres2, Y1, Y2, Y3\?\n\t\t\t"
+			"Real outputs: kappa_C, f_cc, f_s^2, f_s / Q, \?\?, \?\?",
+			0, 2, 0,
+			G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+		)
+	);
+	g_object_class_install_property(
+		gobject_class,
+		ARG_FREQ1,
+		g_param_spec_double(
+			"freq1",
+			"Frequency 1",
+			"First Pcal line frequency, typically around 15-20 Hz.",
+			-G_MAXDOUBLE, G_MAXDOUBLE, G_MAXDOUBLE,
+			G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+		)
+	);
+	g_object_class_install_property(
+		gobject_class,
+		ARG_FREQ2,
+		g_param_spec_double(
+			"freq2",
+			"Frequency 2",
+			"Second Pcal line frequency, typically around 400 Hz.",
+			-G_MAXDOUBLE, G_MAXDOUBLE, G_MAXDOUBLE,
+			G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+		)
+	);
+	g_object_class_install_property(
+		gobject_class,
+		ARG_FREQ4,
+		g_param_spec_double(
+			"freq4",
+			"Frequency 4",
+			"Fourth Pcal line frequency, typically around 10 Hz.",
+			-G_MAXDOUBLE, G_MAXDOUBLE, G_MAXDOUBLE,
+			G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+		)
+	);
+}
+
+
+/*
+ * init()
+ */
+
+
+static void gstlal_sensingtdcfs_init(GSTLALSensingTDCFs *element)
+{
+	gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(element), TRUE);
+	element->rate = 0;
+	element->channels_in = 0;
+	element->channels_out = 0;
+	element->unit_size_in = 0;
+	element->unit_size_out = 0;
+}
+
diff --git a/gstlal-calibration/gst/lal/gstlal_sensingtdcfs.h b/gstlal-calibration/gst/lal/gstlal_sensingtdcfs.h
new file mode 100644
index 0000000000..2f3230d7cd
--- /dev/null
+++ b/gstlal-calibration/gst/lal/gstlal_sensingtdcfs.h
@@ -0,0 +1,95 @@
+/*
+ * Copyright (C) 2019  Aaron Viets <aaron.viets@ligo.org>
+ *
+ * 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.
+ */
+
+
+#ifndef __GSTLAL_SENSINGTDCFS_H__
+#define __GSTLAL_SENSINGTDCFS_H__
+
+#include <glib.h>
+#include <gst/gst.h>
+#include <gst/base/gstbasetransform.h>
+
+G_BEGIN_DECLS
+#define GSTLAL_SENSINGTDCFS_TYPE \
+	(gstlal_sensingtdcfs_get_type())
+#define GSTLAL_SENSINGTDCFS(obj) \
+	(G_TYPE_CHECK_INSTANCE_CAST((obj), GSTLAL_SENSINGTDCFS_TYPE, GSTLALSensingTDCFs))
+#define GSTLAL_SENSINGTDCFS_CLASS(klass) \
+	(G_TYPE_CHECK_CLASS_CAST((klass), GSTLAL_SENSINGTDCFS_TYPE, GSTLALSensingTDCFsClass))
+#define GST_IS_GSTLAL_SENSINGTDCFS(obj) \
+	(G_TYPE_CHECK_INSTANCE_TYPE((obj), GSTLAL_SENSINGTDCFS_TYPE))
+#define GST_IS_GSTLAL_SENSINGTDCFS_CLASS(klass) \
+	(G_TYPE_CHECK_CLASS_TYPE((klass), GSTLAL_SENSINGTDCFS_TYPE))
+
+
+typedef struct _GSTLALSensingTDCFs GSTLALSensingTDCFs;
+typedef struct _GSTLALSensingTDCFsClass GSTLALSensingTDCFsClass;
+
+
+/**
+ * GSTLALSensingTDCFs:
+ */
+
+
+struct _GSTLALSensingTDCFs {
+	GstBaseTransform element;
+
+	/* stream info */
+	gint rate;
+	gint channels_in;
+	gint channels_out;
+	gint unit_size_in;
+	gint unit_size_out;
+	enum gstlal_sensingtdcfs_data_type {
+		GSTLAL_SENSINGTDCFS_FLOAT = 0,
+		GSTLAL_SENSINGTDCFS_DOUBLE
+	} data_type;
+
+	/* timestamp book-keeping */
+	GstClockTime t0;
+	guint64 offset0;
+	guint64 next_in_offset;
+	guint64 next_out_offset;
+	gboolean need_discont;
+
+	/* properties */
+	gint sensing_model;
+	double freq1;
+	double freq2;
+	double freq4;
+};
+
+
+/**
+ * GSTLALSensingTDCFsClass:
+ * @parent_class:  the parent class
+ */
+
+
+struct _GSTLALSensingTDCFsClass {
+	GstBaseTransformClass parent_class;
+};
+
+
+GType gstlal_sensingtdcfs_get_type(void);
+
+
+G_END_DECLS
+
+
+#endif	/* __GSTLAL_SENSINGTDCFS_H__ */
diff --git a/gstlal-calibration/gst/lal/gstlalcalibration.c b/gstlal-calibration/gst/lal/gstlalcalibration.c
index 70881ac0a7..0f9da3f097 100644
--- a/gstlal-calibration/gst/lal/gstlalcalibration.c
+++ b/gstlal-calibration/gst/lal/gstlalcalibration.c
@@ -72,6 +72,7 @@
 #include <gstlal_property.h>
 #include <gstlal_typecast.h>
 #include <gstlal_matrixsolver.h>
+#include <gstlal_sensingtdcfs.h>
 
 
 /*
@@ -108,6 +109,7 @@ static gboolean plugin_init(GstPlugin *plugin)
 		{"lal_property", GSTLAL_PROPERTY_TYPE},
 		{"lal_typecast", GSTLAL_TYPECAST_TYPE},
 		{"lal_matrixsolver", GSTLAL_MATRIXSOLVER_TYPE},
+		{"lal_sensingtdcfs", GSTLAL_SENSINGTDCFS_TYPE},
 		{NULL, 0},
 	};
 
diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py
index a49914db07..fca19f5443 100644
--- a/gstlal-calibration/python/calibration_parts.py
+++ b/gstlal-calibration/python/calibration_parts.py
@@ -1062,7 +1062,7 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate):
 	MV_matrix = mkinterleave(pipeline, MV_matrix)
 	MV_matrix = pipeparts.mkcapsfilter(pipeline, MV_matrix, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=%d" % (rate, (2 * num_stages) * (2 * num_stages + 1)))
 	kappas = pipeparts.mkgeneric(pipeline, MV_matrix, "lal_matrixsolver")
-	kappas = list(mkdeinterleave(pipeline, kappas, 2 * num_stages))
+	kappas = list(mkdeinterleave(pipeline, pipeparts.mkcapsfilter(pipeline, kappas, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=%d" % (rate, 2 * num_stages)), 2 * num_stages))
 	for i in range(len(kappas)):
 		kappas[i] = pipeparts.mkcapsfilter(pipeline, kappas[i], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)
 		if i >= len(kappas) / 2:
diff --git a/gstlal-calibration/tests/check_calibration/Makefile b/gstlal-calibration/tests/check_calibration/Makefile
index 2a12134d5f..c4a670712e 100644
--- a/gstlal-calibration/tests/check_calibration/Makefile
+++ b/gstlal-calibration/tests/check_calibration/Makefile
@@ -8,11 +8,11 @@ IFO = H
 # determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
 OBSRUN = O3
 
-START = $(shell echo 1240617600 - 256 | bc)
+START = $(shell echo 1246226240 - 256 | bc)
 #1238997618
-END = $(shell echo 1240617600 + 1024 + 256 | bc)
+END = $(shell echo 1246226240 + 120 + 256 | bc)
 #1240876818
-SHMRUNTIME = 36000
+SHMRUNTIME = 600
 # How much time does the calibration need to settle at the start and end?
 PLOT_WARMUP_TIME = 256
 PLOT_COOLDOWN_TIME = 256
@@ -25,13 +25,13 @@ DCSLINESCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_CleaningLi
 #../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
 DCSFCCCONFIGS = ../../config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini
 GDSTESTCONFIGS = Filters/O3/GDSFilters/H1GDS_1239476409_testAllCorrections.ini
-DCSTESTCONFIGS = Filters/O3/GDSFilters/L1DCS_C01_1239579018_test.ini
-GDSSHMCONFIGS = Filters/ER14/GDSFilters/H1GDS_1234630818_latency_test.ini
+DCSTESTCONFIGS = Filters/O3/GDSFilters/H1DCS_C01_1239472998_AllCorrectionsTest.ini
+GDSSHMCONFIGS = Filters/O3/GDSFilters/H1GDS_lowLatency_1239476409.ini
 GDSOLDCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_old.ini
 GDSBETTERCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_better.ini
 GDSBESTCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_best.ini
 
-all: highpass_filter_ASD_GDS filters_tf_GDS
+all: pcal_DCS_transfer_functions
 
 ###############################################
 ### These commands should change less often ###
@@ -123,8 +123,9 @@ DCS_pcal2darm_plots:
 	#$(IFO)1_easy_raw_frames.cache $(IFO)1_hoft_DCS_frames.cache $(IFO)1_C00_frames.cache
 	python pcal2darm_timeseries.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --raw-frame-cache $(IFO)1_easy_raw_frames.cache --gstlal-frame-cache-list '$(IFO)1_C00_frames.cache,$(IFO)1_hoft_DCS_frames.cache' --config-file '$(DCSCONFIGS)' --pcal-channel-name CAL-PCALY_RX_PD_OUT_DQ --gstlal-channel-list 'GDS-CALIB_STRAIN,DCS-CALIB_STRAIN' --labels 'C00,TEST' --pcal-time-advance 0.00006103515625
 
-pcal_DCS_transfer_functions: $(IFO)1_easy_raw_frames.cache $(IFO)1_hoft_DCS_frames.cache $(IFO)1_hoft_DCS_FCC_frames.cache $(IFO)1_hoft_DCS_TEST_frames.cache
-	python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_TX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_hoft_DCS_frames.cache,$(IFO)1_hoft_DCS_FCC_frames.cache,$(IFO)1_hoft_DCS_TEST_frames.cache --numerator-channel-list DCS-CALIB_STRAIN_CLEAN,DCS-CALIB_STRAIN_CLEAN,DCS-CALIB_STRAIN_CLEAN --config-file $(DCSCONFIGS) --use-median --labels 'Only scalar corrections,Cavity pole corrections,All corrections'
+pcal_DCS_transfer_functions: 
+	#$(IFO)1_hoft_DCS_frames.cache $(IFO)1_hoft_DCS_TEST_frames.cache
+	python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_RX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_hoft_DCS_frames.cache,$(IFO)1_hoft_DCS_TEST_frames.cache --numerator-channel-list DCS-CALIB_STRAIN,DCS-CALIB_STRAIN --config-file $(DCSCONFIGS) --use-median --labels 'Without SRC Corrections,With SRC Corrections'
 
 lines_ratio_DCS: $(IFO)1_hoft_DCS_frames.cache
 	python demod_ratio_timeseries.py --ifo $(IFO)1 --gps-end-time $(PLOT_END) --gps-start-time $(PLOT_START) --denominator-frame-cache $(IFO)1_hoft_DCS_frames.cache --numerator-frame-cache $(IFO)1_hoft_DCS_frames.cache --denominator-channel-name 'DCS-CALIB_STRAIN' --numerator-channel-name 'DCS-CALIB_STRAIN_CLEAN' --frequencies '35.9,36.7,331.9,1083.7;60,120,180' --magnitude-ranges '0.0,0.1;0.0,1.0' --phase-ranges '-180.0,180.0;-180.0,180.0' --plot-titles '$(IFO)1 Calibration Line Subtraction;$(IFO)1 Power Mains Line Subtraction'
diff --git a/gstlal-calibration/tests/check_calibration/plot_transfer_function.py b/gstlal-calibration/tests/check_calibration/plot_transfer_function.py
index be94dfda6c..37c8519fc9 100644
--- a/gstlal-calibration/tests/check_calibration/plot_transfer_function.py
+++ b/gstlal-calibration/tests/check_calibration/plot_transfer_function.py
@@ -69,7 +69,7 @@ parser.add_option("--gps-end-time", metavar = "seconds", type = int, help = "GPS
 parser.add_option("--ifo", metavar = "name", type = str, help = "Name of the interferometer (IFO), e.g., H1, L1")
 parser.add_option("--denominator-frame-cache", metavar = "name", type = str, help = "Name of frame cache file that contains denominator of transfer functions")
 parser.add_option("--denominator-channel-name", metavar = "name", type = str, default = None, help = "Channel-name of denominator")
-parser.add_option("--denominator-name", metavar = "name", type = str, default = 'Pcal(f)', help = "Name of denominator in plot title, in latex math mode")
+parser.add_option("--denominator-name", metavar = "name", type = str, default = '{\\rm Pcal}(f)', help = "Name of denominator in plot title, in latex math mode")
 parser.add_option("--denominator-correction", metavar = "name", type = str, default = None, help = "Name of filters-file parameter needed to apply a correction to the denominator")
 parser.add_option("--numerator-frame-cache-list", metavar = "name", type = str, help = "Comma-separated list of frame cache files that contain numerators of transfer functions")
 parser.add_option("--numerator-channel-list", metavar = "list", type = str, default = None, help = "Comma-separated list of channel-names of numerators")
@@ -177,7 +177,7 @@ if options.config_file is not None:
 num_corr = []
 denom_corr = []
 if options.numerator_correction is not None:
-	if len(filters[options.numerator_correction]) > 1:
+	if numpy.size(filters[options.numerator_correction]) > 1:
 		corr = filters[options.numerator_correction]
 		# Check the frequency spacing of the correction
 		corr_df = corr[0][1] - corr[0][0]
@@ -199,7 +199,7 @@ if options.numerator_correction is not None:
 			num_corr.append(corr[1][before_idx] + 1j * corr[2][before_idx])
 
 if options.denominator_correction is not None:
-	if len(filters[options.denominator_correction]) > 1:
+	if numpy.size(filters[options.denominator_correction]) > 1:
 		corr = filters[options.denominator_correction]
 		# Check the frequency spacing of the correction
 		corr_df = corr[0][1] - corr[0][0]
@@ -243,7 +243,7 @@ def plot_transfer_function(pipeline, name):
 	denominator = calibration_parts.caps_and_progress(pipeline, denominator, "audio/x-raw,format=F64LE", "denominator")
 	denominator = calibration_parts.mkresample(pipeline, denominator, 5, False, int(sample_rate))
 	if options.denominator_correction is not None:
-		if len(filters[options.denominator_correction]) == 1:
+		if numpy.size(filters[options.denominator_correction]) == 1:
 			denominator = pipeparts.mkaudioamplify(pipeline, denominator, float(filters[options.denominator_correction]))
 	denominator = pipeparts.mktee(pipeline, denominator)
 
@@ -251,11 +251,11 @@ def plot_transfer_function(pipeline, name):
 	for i in range(0, len(labels)):
 		numerator = pipeparts.mklalcachesrc(pipeline, location = numerator_frame_cache_list[i], cache_dsc_regex = ifo)
 		numerator = pipeparts.mkframecppchanneldemux(pipeline, numerator, do_file_checksum = False, skip_bad_files = True, channel_list = map("%s:%s".__mod__, channel_list))
-		numerator = calibration_parts.hook_up(pipeline, numerator, numerator_channel_list[i], ifo, 1.0)
+		numerator = calibration_parts.hook_up(pipeline, numerator, numerator_channel_list[i], ifo, 1.0, element_name_suffix = "%d" % i)
 		numerator = calibration_parts.caps_and_progress(pipeline, numerator, "audio/x-raw,format=F64LE", labels[i])
 		numerator = calibration_parts.mkresample(pipeline, numerator, 5, False, int(sample_rate))
 		if options.numerator_correction is not None:
-			if len(filters[options.numerator_correction] == 1):
+			if numpy.size(filters[options.numerator_correction] == 1):
 				numerator = pipeparts.mkaudioamplify(pipeline, numerator, float(filters[options.numerator_correction]))
 		# Interleave the channels to make one stream
 		channels = calibration_parts.mkinterleave(pipeline, [numerator, denominator])
@@ -294,7 +294,7 @@ if options.poles is not None:
 	for i in range(0, len(real_poles) / 2):
 		poles.append(float(real_poles[2 * i]) + 1j * float(real_poles[2 * i + 1]))
 
-colors = ['limegreen', 'g', 'y', 'c', 'm', 'b'] # Hopefully the user will not want to plot more than six datasets on one plot.
+colors = ['blue', 'limegreen', 'y', 'c', 'm', 'b'] # Hopefully the user will not want to plot more than six datasets on one plot.
 for i in range(0, len(labels)):
 	# Remove unwanted lines from file, and re-format wanted lines
 	f = open('%s_%s_over_%s_%d-%d.txt' % (ifo, labels[i].replace(' ', '_').replace('/', 'over'), options.denominator_channel_name, options.gps_start_time, data_duration),"r")
@@ -330,13 +330,13 @@ for i in range(0, len(labels)):
 	if i == 0:
 		plt.figure(figsize = (10, 10))
 	plt.subplot(211)
-	plt.plot(frequency, magnitude, colors[i % 6], linewidth = 0.75, label = labels[i].replace('_', '\_'))
-	#leg = plt.legend(fancybox = True)
-	#leg.get_frame().set_alpha(0.8)
+	plt.plot(frequency, magnitude, colors[i % 6], linewidth = 0.75, label = r'${\rm %s}$' % labels[i].replace('_', '\_').replace(' ', '\ '))
+	leg = plt.legend(fancybox = True)
+	leg.get_frame().set_alpha(0.8)
 	plt.gca().set_xscale(freq_scale)
 	plt.gca().set_yscale(mag_scale)
 	if i == 0:
-		#plt.title(r'%s $%s$ / $%s$' % (ifo, options.numerator_name.replace('_', '\_'), options.denominator_name.replace('_', '\_')))
+		plt.title(r'${\rm %s} \ %s \ / \ %s$' % (ifo, options.numerator_name, options.denominator_name))
 		plt.ylabel(r'${\rm Magnitude}$')
 	plt.xlim(options.frequency_min, options.frequency_max)
 	plt.ylim(options.magnitude_min, options.magnitude_max)
@@ -350,6 +350,6 @@ for i in range(0, len(labels)):
 	plt.xlim(options.frequency_min, options.frequency_max)
 	plt.ylim(options.phase_min, options.phase_max)
 	plt.grid(True, which = "both", linestyle = ':', linewidth = 0.3, color = 'black')
-plt.savefig('%s_transfer_functions_%s_%s%s_%d-%d.png' % (ifo, options.numerator_name, options.denominator_name, options.filename_suffix, options.gps_start_time, data_duration))
-plt.savefig('%s_transfer_functions_%s_%s%s_%d-%d.pdf' % (ifo, options.numerator_name, options.denominator_name, options.filename_suffix, options.gps_start_time, data_duration))
+plt.savefig('%s_%s_over_%s_%d-%d.pdf' % (ifo, numerator_channel_list[-1], options.denominator_channel_name, options.gps_start_time, data_duration))
+plt.savefig('%s_%s_over_%s_%d-%d.png' % (ifo, numerator_channel_list[-1], options.denominator_channel_name, options.gps_start_time, data_duration))
 
-- 
GitLab