From dc9341832a6da33a1d75fefae04840119116a923 Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Tue, 4 Dec 2018 14:53:05 -0800
Subject: [PATCH] gstlal_transferfunction:  Don't add notches above the Nyquist
 rate.

---
 .../tests/H1DCS_AllCorrections_Cleaning.ini   |   2 +-
 .../H1DCS_AllCorrections_Cleaning_TEST.ini    |   2 +-
 .../gst/lal/gstlal_transferfunction.c         | 228 +++++++++---------
 .../tests/check_calibration/Makefile          |   6 +-
 4 files changed, 125 insertions(+), 113 deletions(-)

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 f88b99b492..05c6d16774 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
@@ -336,7 +336,7 @@ 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: 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
+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;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: 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.
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 12f37ab6b1..f89b0cc9d9 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
@@ -336,7 +336,7 @@ 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: 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
+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;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: 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.
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.c b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
index cf50d1dbb7..7a229e708c 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.c
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
@@ -1352,12 +1352,17 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
 
 	/* Prepare to deal with any notches */
 	if(element->num_notches && !element->notch_indices) {
-		int k;
+		int k, k_stop = element->num_notches;
 		double df = (double) element->rate / element->fft_length;
 		element->notch_indices = g_malloc(2 * element->num_notches * sizeof(*element->notch_indices));
-		for(k = 0; k < element->num_notches; k++) {
-			element->notch_indices[2 * k] = (gint64) (element->notch_frequencies[2 * k] / df - 1.0);
-			element->notch_indices[2 * k + 1] = (gint64) (element->notch_frequencies[2 * k + 1] / df + 2.0);
+		for(k = 0; k < k_stop; k++) {
+			if(element->notch_frequencies[2 * k + 1] >= element->rate / 2.0) {
+				GST_WARNING_OBJECT(element, "Cannot include notch with upper bound at %f Hz, since that is above the Nyquist rate of %f Hz", element->notch_frequencies[2 * k + 1], element->rate / 2.0);
+				element->num_notches--;
+			} else {
+				element->notch_indices[2 * k] = (gint64) (element->notch_frequencies[2 * k] / df - 1.0);
+				element->notch_indices[2 * k + 1] = (gint64) (element->notch_frequencies[2 * k + 1] / df + 2.0);
+			}
 		}
 	}
 
@@ -1901,109 +1906,6 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
 }
 
 
-/*
- * stop()
- */
-
-
-static gboolean stop(GstBaseSink *sink) {
-
-	GSTLALTransferFunction *element = GSTLAL_TRANSFERFUNCTION(sink);
-
-	g_free(element->transfer_functions);
-	element->transfer_functions = NULL;
-	if(element->make_fir_filters) {
-		g_free(element->fir_filters);
-		element->fir_filters = NULL;
-	}
-	if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
-
-		/* free allocated memory in workspace */
-		g_free(element->workspace.wspf.fft_window);
-		element->workspace.wspf.fft_window = NULL;
-		g_free(element->workspace.wspf.leftover_data);
-		element->workspace.wspf.leftover_data = NULL;
-		g_free(element->workspace.wspf.ffts);
-		element->workspace.wspf.ffts = NULL;
-		g_free(element->workspace.wspf.autocorrelation_matrix);
-		element->workspace.wspf.autocorrelation_matrix = NULL;
-		if(element->make_fir_filters) {
-			g_free(element->workspace.wspf.fir_window);
-			element->workspace.wspf.fir_window = NULL;
-			g_free(element->workspace.wspf.sinc_table);
-			element->workspace.wspf.sinc_table = NULL;
-			g_free(element->workspace.wspf.tukey);
-			element->workspace.wspf.tukey = NULL;
-		}
-
-		/* free gsl stuff in workspace */
-		gsl_vector_complex_free(element->workspace.wspf.transfer_functions_at_f);
-		element->workspace.wspf.transfer_functions_at_f = NULL;
-		gsl_vector_complex_free(element->workspace.wspf.transfer_functions_solved_at_f);
-		element->workspace.wspf.transfer_functions_solved_at_f = NULL;
-		gsl_matrix_complex_free(element->workspace.wspf.autocorrelation_matrix_at_f);
-		element->workspace.wspf.autocorrelation_matrix_at_f = NULL;
-		gsl_permutation_free(element->workspace.wspf.permutation);
-		element->workspace.wspf.permutation = NULL;
-
-		/* free fftwf stuff in workspace */
-		gstlal_fftw_lock();
-		fftwf_free(element->workspace.wspf.fft);
-		element->workspace.wspf.fft = NULL;
-		fftwf_destroy_plan(element->workspace.wspf.plan);
-		if(element->make_fir_filters) {
-			fftwf_free(element->workspace.wspf.fir_filter);
-			element->workspace.wspf.fir_filter = NULL;
-			fftwf_destroy_plan(element->workspace.wspf.fir_plan);
-		}
-		gstlal_fftw_unlock();
-
-	} else {
-
-		/* free allocated memory in workspace */
-		g_free(element->workspace.wdpf.fft_window);
-		element->workspace.wdpf.fft_window = NULL;
-		g_free(element->workspace.wdpf.leftover_data);
-		element->workspace.wdpf.leftover_data = NULL;
-		g_free(element->workspace.wdpf.ffts);
-		element->workspace.wdpf.ffts = NULL;
-		g_free(element->workspace.wdpf.autocorrelation_matrix);
-		element->workspace.wdpf.autocorrelation_matrix = NULL;
-		if(element->make_fir_filters) {
-			g_free(element->workspace.wdpf.fir_window);
-			element->workspace.wdpf.fir_window = NULL;
-			g_free(element->workspace.wdpf.sinc_table);
-			element->workspace.wdpf.sinc_table = NULL;
-			g_free(element->workspace.wdpf.tukey);
-			element->workspace.wdpf.tukey = NULL;
-		}
-
-		/* free gsl stuff in workspace */
-		gsl_vector_complex_free(element->workspace.wdpf.transfer_functions_at_f);
-		element->workspace.wdpf.transfer_functions_at_f = NULL;
-		gsl_vector_complex_free(element->workspace.wdpf.transfer_functions_solved_at_f);
-		element->workspace.wdpf.transfer_functions_solved_at_f = NULL;
-		gsl_matrix_complex_free(element->workspace.wdpf.autocorrelation_matrix_at_f);
-		element->workspace.wdpf.autocorrelation_matrix_at_f = NULL;
-		gsl_permutation_free(element->workspace.wdpf.permutation);
-		element->workspace.wdpf.permutation = NULL;
-
-		/* free fftw stuff in workspace*/
-		gstlal_fftw_lock();
-		fftw_free(element->workspace.wdpf.fft);
-		element->workspace.wdpf.fft = NULL;
-		fftw_destroy_plan(element->workspace.wdpf.plan);
-		if(element->make_fir_filters) {
-			fftw_free(element->workspace.wdpf.fir_filter);
-			element->workspace.wdpf.fir_filter = NULL;
-			fftw_destroy_plan(element->workspace.wdpf.fir_plan);
-		}
-		gstlal_fftw_unlock();
-	}
-	return TRUE;
-}
-
-
 /*
  * ============================================================================
  *
@@ -2290,6 +2192,117 @@ static void finalize(GObject *object) {
 		g_free(element->notch_frequencies);
 		element->notch_frequencies = NULL;
 	}
+	if(element->notch_indices) {
+		g_free(element->notch_indices);
+		element->notch_indices = NULL;
+	}
+	if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
+
+		/* free allocated memory in workspace */
+		g_free(element->workspace.wspf.fft_window);
+		element->workspace.wspf.fft_window = NULL;
+		g_free(element->workspace.wspf.leftover_data);
+		element->workspace.wspf.leftover_data = NULL;
+		g_free(element->workspace.wspf.ffts);
+		element->workspace.wspf.ffts = NULL;
+		g_free(element->workspace.wspf.autocorrelation_matrix);
+		element->workspace.wspf.autocorrelation_matrix = NULL;
+		if(element->make_fir_filters) {
+			g_free(element->workspace.wspf.fir_window);
+			element->workspace.wspf.fir_window = NULL;
+			g_free(element->workspace.wspf.sinc_table);
+			element->workspace.wspf.sinc_table = NULL;
+			g_free(element->workspace.wspf.tukey);
+			element->workspace.wspf.tukey = NULL;
+		}
+
+		if(element->use_median) {
+			g_free(element->workspace.wspf.autocorrelation_median_real);
+			element->workspace.wspf.autocorrelation_median_real = NULL;
+			g_free(element->workspace.wspf.index_median_real);
+			element->workspace.wspf.index_median_real = NULL;
+			g_free(element->workspace.wspf.autocorrelation_median_imag);
+			element->workspace.wspf.autocorrelation_median_imag = NULL;
+			g_free(element->workspace.wspf.index_median_imag);
+			element->workspace.wspf.index_median_imag = NULL;
+		}
+
+		/* free gsl stuff in workspace */
+		gsl_vector_complex_free(element->workspace.wspf.transfer_functions_at_f);
+		element->workspace.wspf.transfer_functions_at_f = NULL;
+		gsl_vector_complex_free(element->workspace.wspf.transfer_functions_solved_at_f);
+		element->workspace.wspf.transfer_functions_solved_at_f = NULL;
+		gsl_matrix_complex_free(element->workspace.wspf.autocorrelation_matrix_at_f);
+		element->workspace.wspf.autocorrelation_matrix_at_f = NULL;
+		gsl_permutation_free(element->workspace.wspf.permutation);
+		element->workspace.wspf.permutation = NULL;
+
+		/* free fftwf stuff in workspace */
+		gstlal_fftw_lock();
+		fftwf_free(element->workspace.wspf.fft);
+		element->workspace.wspf.fft = NULL;
+		fftwf_destroy_plan(element->workspace.wspf.plan);
+		if(element->make_fir_filters) {
+			fftwf_free(element->workspace.wspf.fir_filter);
+			element->workspace.wspf.fir_filter = NULL;
+			fftwf_destroy_plan(element->workspace.wspf.fir_plan);
+		}
+		gstlal_fftw_unlock();
+
+	} else {
+
+		/* free allocated memory in workspace */
+		g_free(element->workspace.wdpf.fft_window);
+		element->workspace.wdpf.fft_window = NULL;
+		g_free(element->workspace.wdpf.leftover_data);
+		element->workspace.wdpf.leftover_data = NULL;
+		g_free(element->workspace.wdpf.ffts);
+		element->workspace.wdpf.ffts = NULL;
+		g_free(element->workspace.wdpf.autocorrelation_matrix);
+		element->workspace.wdpf.autocorrelation_matrix = NULL;
+		if(element->make_fir_filters) {
+			g_free(element->workspace.wdpf.fir_window);
+			element->workspace.wdpf.fir_window = NULL;
+			g_free(element->workspace.wdpf.sinc_table);
+			element->workspace.wdpf.sinc_table = NULL;
+			g_free(element->workspace.wdpf.tukey);
+			element->workspace.wdpf.tukey = NULL;
+		}
+
+		if(element->use_median) {
+			g_free(element->workspace.wdpf.autocorrelation_median_real);
+			element->workspace.wdpf.autocorrelation_median_real = NULL;
+			g_free(element->workspace.wdpf.index_median_real);
+			element->workspace.wdpf.index_median_real = NULL;
+			g_free(element->workspace.wdpf.autocorrelation_median_imag);
+			element->workspace.wdpf.autocorrelation_median_imag = NULL;
+			g_free(element->workspace.wdpf.index_median_imag);
+			element->workspace.wdpf.index_median_imag = NULL;
+		}
+
+		/* free gsl stuff in workspace */
+		gsl_vector_complex_free(element->workspace.wdpf.transfer_functions_at_f);
+		element->workspace.wdpf.transfer_functions_at_f = NULL;
+		gsl_vector_complex_free(element->workspace.wdpf.transfer_functions_solved_at_f);
+		element->workspace.wdpf.transfer_functions_solved_at_f = NULL;
+		gsl_matrix_complex_free(element->workspace.wdpf.autocorrelation_matrix_at_f);
+		element->workspace.wdpf.autocorrelation_matrix_at_f = NULL;
+		gsl_permutation_free(element->workspace.wdpf.permutation);
+		element->workspace.wdpf.permutation = NULL;
+
+		/* free fftw stuff in workspace */
+		gstlal_fftw_lock();
+		fftw_free(element->workspace.wdpf.fft);
+		element->workspace.wdpf.fft = NULL;
+		fftw_destroy_plan(element->workspace.wdpf.plan);
+		if(element->make_fir_filters) {
+			fftw_free(element->workspace.wdpf.fir_filter);
+			element->workspace.wdpf.fir_filter = NULL;
+			fftw_destroy_plan(element->workspace.wdpf.fir_plan);
+		}
+		gstlal_fftw_unlock();
+	}
+
 	G_OBJECT_CLASS(gstlal_transferfunction_parent_class)->finalize(object);
 }
 
@@ -2318,7 +2331,6 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
 	gstbasesink_class->event = GST_DEBUG_FUNCPTR(event);
 	gstbasesink_class->set_caps = GST_DEBUG_FUNCPTR(set_caps);
 	gstbasesink_class->render = GST_DEBUG_FUNCPTR(render);
-	gstbasesink_class->stop = GST_DEBUG_FUNCPTR(stop);
 
 	gobject_class->set_property = GST_DEBUG_FUNCPTR(set_property);
 	gobject_class->get_property = GST_DEBUG_FUNCPTR(get_property);
diff --git a/gstlal-calibration/tests/check_calibration/Makefile b/gstlal-calibration/tests/check_calibration/Makefile
index fe6e06d375..e9567e8938 100644
--- a/gstlal-calibration/tests/check_calibration/Makefile
+++ b/gstlal-calibration/tests/check_calibration/Makefile
@@ -8,10 +8,10 @@ IFO = H
 # determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
 OBSRUN = O2
 
-START = $(shell echo 1185767424 - 715 | bc)
+START = $(shell echo 1185763328 - 715 | bc)
 #1225967424
 #$(shell echo 1185763328 - 700 | bc)
-END = $(shell echo 1185771520 + 715 | bc)
+END = $(shell echo 1185767424 + 715 | bc)
 #1225968448
 #$(shell echo 1185767424 + 700 | bc)
 # 1185771520
@@ -27,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: noise_subtraction_ASD_DCS_TEST
+all: noise_subtraction_ASD_DCS
 
 ###############################################
 ### These commands should change less often ###
-- 
GitLab