Optimize line subtraction median array length updates for high latency

In order to mitigate "oversubtraction" of spectral lines and the unwanted narrow-band impact on nearby frequencies in h(t), we use a rather long (4096 seconds) running median to estimate the TFs used to reconstruct the subtraction signal from the witness channels. However, at times when those TFs undergo significant secular change (e.g., thermalizations at the beginnings of lock stretches), this could lead to ineffective subtraction. Therefore, we use a method that detects secular changes in the TFs and adjusts the length of the median accordingly. During stable times, the full 4096s are used. When significant secular change is detected, this is temporarily reduced to 128s. Once the TFs become stable again, the median grows until it again reaches 4096s.

It occurred to me that, when I first implemented this method, I assumed that it would only be used in low latency. This means that, when the running median is shortened to 128s, it uses only the most recent 128s of data, dropping older data. In high latency, since we use future data to estimate TFs, this is not desirable behavior. In this case, we want to retain the oldest data in the running median; that is, the data that is temporally closest to the data we are subtracting. This MR should ensure that, regardless of whether we are doing low-latency or high-latency subtraction, the data that is retained in the running median during these times is the appropriate data.

Here is an example of a time at H1 when there were several thermalizations during which the line subtraction running medians were temporarily shortened due to secular changes. Here is a spectrogram from before the fix: Here is a spectrogram from the same stretch of data calibrated (and line subtracted) with this update:

If carefully compared, an improvement can be seen in the 15-20 Hz region.

gst/lal/gstlal_smoothkappas.c

@@ -161,7 +161,7 @@ static void set_metadata(GSTLALSmoothKappas *element, GstBuffer *buf, guint64 ou
 }
 
 
-static void get_new_median(double new_element, double *fifo_array, double *current_median, gint array_size, int *index_re, int *index_im, gboolean array_is_imaginary, int samples_in_filter, gboolean no_default, int *future_unused) {
+static void get_new_median(double new_element, double *fifo_array, double *current_median, gint array_size, int *index_re, int *index_im, gboolean array_is_imaginary, int samples_in_filter, gboolean no_default, int future_unused) {
 	if(array_is_imaginary) {
 		fifo_array[*index_im] = new_element;
 		(*index_im)++;
@@ -180,7 +180,7 @@ static void get_new_median(double new_element, double *fifo_array, double *curre
 	greater = Greater = G_MAXDOUBLE;
 	less = Less = -G_MAXDOUBLE;
 	num_samples = (no_default && samples_in_filter < array_size) ? samples_in_filter : array_size;
-	j_start = (array_is_imaginary ? *index_im : *index_re) + array_size - num_samples - *future_unused;
+	j_start = (array_is_imaginary ? *index_im : *index_re) + array_size - num_samples - future_unused;
 	for(j = j_start; j < j_start + num_samples; j++) {
 		jj = j % array_size;
 		if(fifo_array[jj] < *current_median) {
@@ -204,14 +204,15 @@ static void get_new_median(double new_element, double *fifo_array, double *curre
 		else
 			g_assert_not_reached();
 	}
-	if(*future_unused > 0)
-		(*future_unused)--;
-
 	g_assert_cmpint(number_less + number_equal + number_greater, ==, num_samples);
 
 	if(num_samples % 2) {
-		if(number_greater > num_samples / 2)
+		if(number_greater > 1 + num_samples / 2)
+			*current_median = Greater;
+		else if(number_greater > num_samples / 2)
 			*current_median = greater;
+		else if(number_less > 1 + num_samples / 2)
+			*current_median = Less;
 		else if(number_less > num_samples / 2)
 			*current_median = less;
 	} else {
@@ -272,8 +273,9 @@ static void reset(GSTLALSmoothKappas *element) {
 	g_assert_cmpint(element->future_unused, >=, 0);
 	/* Drop samples to make truncated array */
 	element->samples_in_filter = element->min_array_size;
-	int i, num_default, ii, i_start, i_stop, num_less, num_greater, num_equal;
+	int i, num_default, ii, i_start, i_stop, num_less, num_greater, num_equal, num_samples;
 	double greater, less, replace_value_re, replace_value_im;
+	num_samples = (element->no_default && element->samples_in_filter < element->array_size) ? element->samples_in_filter : element->array_size;
 
 	/* Update median array and bookkeeping */
 	if(element->array_size > element->samples_in_filter) {
@@ -292,10 +294,12 @@ static void reset(GSTLALSmoothKappas *element) {
 		greater = G_MAXDOUBLE;
 		less = -G_MAXDOUBLE;
 		/* Real part */
-		while(num_less + num_equal + num_greater == 0 || num_greater > element->array_size / 2 || num_less > element->array_size / 2) {
+		while(num_less + num_equal + num_greater == 0 || num_greater > num_samples / 2 || num_less > num_samples / 2) {
 			num_less = 0;
 			num_equal = 0;
 			num_greater = 0;
+			greater = G_MAXDOUBLE;
+			less = -G_MAXDOUBLE;
 			i_start = (element->index_re + 2 * element->array_size - element->samples_in_filter - element->future_unused) % element->array_size;
 			i_stop = i_start + (element->no_default ? element->samples_in_filter : element->array_size);
 			for(i = i_start; i < i_stop; i++) {
@@ -312,15 +316,15 @@ static void reset(GSTLALSmoothKappas *element) {
 						greater = element->fifo_array_re[ii];
 				}
 			}
-			if(num_less == element->array_size / 2 && num_greater == element->array_size / 2 && !(element->array_size % 2))
+			if(num_less == num_samples / 2 && num_greater == num_samples / 2 && !(num_samples % 2))
 				element->current_median_re = (less + greater) / 2;
-			else if(num_less == element->array_size / 2 && !(element->array_size % 2))
+			else if(num_less == num_samples / 2 && !(num_samples % 2))
 				element->current_median_re = (less + element->current_median_re) / 2;
-			else if(num_greater == element->array_size / 2 && !(element->array_size % 2))
+			else if(num_greater == num_samples / 2 && !(num_samples % 2))
 				element->current_median_re = (element->current_median_re + greater) / 2;
-			else if(num_greater > element->array_size / 2)
+			else if(num_greater > num_samples / 2)
 				element->current_median_re = greater;
-			else if(num_less > element->array_size / 2)
+			else if(num_less > num_samples / 2)
 				element->current_median_re = less;
 			/* In all other cases, do nothing */
 		}
@@ -330,10 +334,12 @@ static void reset(GSTLALSmoothKappas *element) {
 		greater = G_MAXDOUBLE;
 		less = -G_MAXDOUBLE;
 		/* Imaginary part */
-		while(num_less + num_equal + num_greater == 0 || num_greater > element->array_size / 2 || num_less > element->array_size / 2) {
+		while(num_less + num_equal + num_greater == 0 || num_greater > num_samples / 2 || num_less > num_samples / 2) {
 			num_less = 0;
 			num_equal = 0;
 			num_greater = 0;
+			greater = G_MAXDOUBLE;
+			less = -G_MAXDOUBLE;
 			i_start = (element->index_im + 2 * element->array_size - element->samples_in_filter - element->future_unused) % element->array_size;
 			i_stop = i_start + (element->no_default ? element->samples_in_filter : element->array_size);
 			for(i = i_start; i < i_stop; i++) {
@@ -350,15 +356,15 @@ static void reset(GSTLALSmoothKappas *element) {
 						greater = element->fifo_array_im[ii];
 				}
 			}
-			if(num_less == element->array_size / 2 && num_greater == element->array_size / 2 && !(element->array_size % 2))
+			if(num_less == num_samples / 2 && num_greater == num_samples / 2 && !(num_samples % 2))
 				element->current_median_im = (less + greater) / 2;
-			else if(num_less == element->array_size / 2 && !(element->array_size % 2))
+			else if(num_less == num_samples / 2 && !(num_samples % 2))
 				element->current_median_im = (less + element->current_median_im) / 2;
-			else if(num_greater == element->array_size / 2 && !(element->array_size % 2))
+			else if(num_greater == num_samples / 2 && !(num_samples % 2))
 				element->current_median_im = (element->current_median_im + greater) / 2;
-			else if(num_greater > element->array_size / 2)
+			else if(num_greater > num_samples / 2)
 				element->current_median_im = greater;
-			else if(num_less > element->array_size / 2)
+			else if(num_less > num_samples / 2)
 				element->current_median_im = less;
 			/* In all other cases, do nothing */
 		}
@@ -370,6 +376,9 @@ static void reset(GSTLALSmoothKappas *element) {
 				element->fifo_array_im[i % element->array_size] = element->current_median_im;
 		}
 	}
+	if(element->future_unused > 0)
+		(element->future_unused)--;
+
 	/* Update average array and bookkeeping */
 	if(element->avg_array_size > element->samples_in_filter) {
 		num_default = element->avg_array_size - element->min_array_size;
@@ -405,7 +414,9 @@ static GstFlowReturn smooth_buffer_ ## DTYPE(const DTYPE *src, guint64 src_size,
 		} \
  \
 		/* Compute new median */ \
-		get_new_median(new_element, fifo_array, current_median, array_size, index_re, index_im, FALSE, *samples_in_filter + (int) i + 1, no_default, future_unused); \
+		get_new_median(new_element, fifo_array, current_median, array_size, index_re, index_im, FALSE, *samples_in_filter + (int) i + 1, no_default, *future_unused); \
+		if(*future_unused > 0) \
+			(*future_unused)--; \
 		/* Compute new average */ \
 		new_avg = (DTYPE) get_average(*current_median, avg_array, avg_array_size, avg_index_re, avg_index_im, FALSE, *samples_in_filter + (int) i + 1, no_default); \
 		if (track_bad_kappa) { \
@@ -463,8 +474,10 @@ static GstFlowReturn smooth_complex_buffer_ ## DTYPE(const DTYPE complex *src, g
 			*num_bad_in_avg_im = 0; \
 		} \
 		/* Compute new median */ \
-		get_new_median(new_element_re, fifo_array_re, current_median_re, array_size, index_re, index_im, FALSE, *samples_in_filter + (int) i + 1, no_default, future_unused); \
-		get_new_median(new_element_im, fifo_array_im, current_median_im, array_size, index_re, index_im, TRUE, *samples_in_filter + (int) i + 1, no_default, future_unused); \
+		get_new_median(new_element_re, fifo_array_re, current_median_re, array_size, index_re, index_im, FALSE, *samples_in_filter + (int) i + 1, no_default, *future_unused); \
+		get_new_median(new_element_im, fifo_array_im, current_median_im, array_size, index_re, index_im, TRUE, *samples_in_filter + (int) i + 1, no_default, *future_unused); \
+		if(*future_unused > 0) \
+			(*future_unused)--; \
 		/* Compute new average */ \
 		new_avg = (DTYPE) get_average(*current_median_re, avg_array_re, avg_array_size, avg_index_re, avg_index_im, FALSE, *samples_in_filter + (int) i + 1, no_default) + I * (DTYPE) get_average(*current_median_im, avg_array_im, avg_array_size, avg_index_re, avg_index_im, TRUE, *samples_in_filter + (int) i + 1, no_default); \
 		if(track_bad_kappa) { \
@@ -791,7 +804,7 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
 	g_mutex_lock(&element->reset_mutex);
 	if(element->num_reset_times > 0) {
 		if(GST_BUFFER_PTS(inbuf) >= *element->reset_times) {
-			if(element->min_array_size < element->samples_in_filter)
+			if(element->min_array_size < element->samples_in_filter && element->min_array_size > 0)
 				reset(element);
 			element->num_reset_times--;
 			for(i = 0; (guint) i < element->num_reset_times; i++)
@@ -817,9 +830,8 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
 		result = smooth_complex_buffer_float((const float complex *) inmap.data, inmap.size / element->unit_size, (float complex *) outmap.data, outmap.size / element->unit_size, element->fifo_array_re, element->fifo_array_im, element->avg_array_re, element->avg_array_im, element->default_kappa_re, element->default_kappa_im, &element->current_median_re, &element->current_median_im, element->maximum_offset_re, element->maximum_offset_im, element->array_size, element->avg_array_size, &element->index_re, &element->index_im, &element->avg_index_re, &element->avg_index_im, &element->num_bad_in_avg_re, &element->num_bad_in_avg_im, gap, element->default_to_median, element->track_bad_kappa, &element->samples_in_filter, element->no_default, &element->total_samples, element->reject_zeros, &element->future_unused);
 	} else if(element->data_type == GSTLAL_SMOOTHKAPPAS_Z128) {
 		result = smooth_complex_buffer_double((const double complex *) inmap.data, inmap.size / element->unit_size, (double complex *) outmap.data, outmap.size / element->unit_size, element->fifo_array_re, element->fifo_array_im, element->avg_array_re, element->avg_array_im, element->default_kappa_re, element->default_kappa_im, &element->current_median_re, &element->current_median_im, element->maximum_offset_re, element->maximum_offset_im, element->array_size, element->avg_array_size, &element->index_re, &element->index_im, &element->avg_index_re, &element->avg_index_im, &element->num_bad_in_avg_re, &element->num_bad_in_avg_im, gap, element->default_to_median, element->track_bad_kappa, &element->samples_in_filter, element->no_default, &element->total_samples, element->reject_zeros, &element->future_unused);
-	} else {
+	} else
 		g_assert_not_reached();
-	}
 
 	set_metadata(element, outbuf, outmap.size / element->unit_size);
 
@@ -936,6 +948,8 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
 				old_reset_time = new_reset_time;
 			}
 			g_mutex_unlock(&element->reset_mutex);
+			if(element->min_array_size == 0)
+				GST_WARNING_OBJECT(element, "received a reset time, but min_array_size is not set, so reset will not occur");
 		}
 		break;
 	case ARG_REJECT_ZEROS:
@@ -1079,7 +1093,7 @@ static void gstlal_smoothkappas_class_init(GSTLALSmoothKappasClass *klass)
 			"array-size",
 			"Median array size",
 			"Size of the array of values from which the median is calculated.",
-			G_MININT, G_MAXINT, 2049,
+			0, G_MAXINT, 2049,
 			G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
 		)
 	);
@@ -1090,7 +1104,7 @@ static void gstlal_smoothkappas_class_init(GSTLALSmoothKappasClass *klass)
 			"min-array-size",
 			"Minimum Median Array Size",
 			"Number of samples to be retained in median array when element is reset",
-			G_MININT, G_MAXINT, 0,
+			0, G_MAXINT, 0,
 			G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
 		)
 	);
@@ -1102,7 +1116,7 @@ static void gstlal_smoothkappas_class_init(GSTLALSmoothKappasClass *klass)
 			"Average array size",
 			"Size of the array of values from which the average is calculated\n\t\t\t"
 			"from the median values. By default, no average is taken.",
-			G_MININT, G_MAXINT, 1,
+			0, G_MAXINT, 1,
 			G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
 		)
 	);