From cc396eddb2c771b62a3d238877a7d10c1ee6f648 Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Wed, 27 Jun 2018 14:01:13 -0700
Subject: [PATCH] gstlal_compute_strain: Options for data cleaning. Also added
a test of the behavior of lal_gate.
---
gstlal-calibration/bin/gstlal_compute_strain | 12 ++-
.../python/calibration_parts.py | 4 +-
gstlal-calibration/tests/lal_gate_test.py | 88 +++++++++++++++++++
3 files changed, 101 insertions(+), 3 deletions(-)
create mode 100755 gstlal-calibration/tests/lal_gate_test.py
diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain
index f277bdc4e4..6c10636766 100755
--- a/gstlal-calibration/bin/gstlal_compute_strain
+++ b/gstlal-calibration/bin/gstlal_compute_strain
@@ -283,6 +283,11 @@ parser.add_option("--powerlines-channel-sr", metavar = "Hz", type = int, default
parser.add_option("--powerlines-freq-var", metavar = "Hz", type = float, default = 0.02, help = "Amount by which freqency of the power mains varies with time. Supposedly they are intentionally varied by +- 0.02 Hz. (Default = 0.02)")
parser.add_option("--witness-channel-list", metavar = "name", default = None, help = "Comma-separated list of witness channels to use to subtract noise from h(t). (Default is to not use any witness channels.)")
parser.add_option("--witness-channel-sr", metavar = "Hz", type = int, default = 2048, help = "Sample rate at which transfer functions will be computed and witness channels will be filtered. (Default = 2048 Hz)")
+parser.add_option("--witness-channel-fft-length", metavar = "seconds", type = float, default = 1.0, help = "The length in seconds of the fast Fourier transforms used to compute transfer functions between witness channels and h(t). The fft's are windowed with Hann windows and overlapped. (Default = 1.0)")
+parser.add_option("--num-witness-ffts", type = int, default = 128, help = "The number of fft's to average for the witness -> h(t) transfer functions calculation. The average is taken after the ratio h(f) / witness(f). (Default = 128)")
+parser.add_option("--witness-fir-length", metavar = "seconds", type = float, default = 0.5, help = "The length in seconds of the filters applied to the witness channels before subtracting from h(t) (Default = 0.5)")
+parser.add_option("--witness-frequency-resolution", metavar = "Hz", type = float, default = 5.0, help = "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. (Default = 5.0)")
+parser.add_option("--witness-tf-update-time", metavar = "seconds", type = int, default = 3600, help = "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 (Default = 3600)")
parser.add_option("--cleaning-check-rms-time", metavar = "seconds", type = float, default = 1.0, help = "The amount of data from h(t) and cleaned h(t) that is used to compute and compare the rms. This comparison between cleaned and uncleaned h(t) determines whether the HOFT_CLEAN bits of the calibration state vector are on or off.")
parser.add_option("--cleaning-check-range-low-min", metavar = "Hz", type = int, default = 15, help = "Minimum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_LOWFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range. (Default = 15 Hz)")
parser.add_option("--cleaning-check-range-low-max", metavar = "Hz", type = int, default = 40, help = "Maximum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_LOWFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range. (Default = 40 Hz)")
@@ -1939,6 +1944,11 @@ if options.remove_powerlines:
# Remove excess noise using any provided witness channels
if options.witness_channel_list is not None:
+ witness_fft_samples = int(options.witness_channel_fft_length * options.witness_channel_sr)
+ witness_fft_overlap = int(witness_fft_samples / 2)
+ witness_tf_update_samples = int(options.witness_channel_sr) * int(options.witness_tf_update_time)
+ witness_fir_samples = int(options.witness_fir_length * options.witness_channel_sr)
+
if not (options.remove_callines or options.remove_powerlines):
clean_strain = straintee
if options.no_dq_vector:
@@ -1949,7 +1959,7 @@ if options.witness_channel_list is not None:
witnesses.append(calibration_parts.caps_and_progress(pipeline, head_dict[key], "audio/x-raw, format=F64LE, channels=1, channel-mask=(bitmask)0x0", key))
if len(witnesses) != len(witness_channel_list):
print("WARNING: Not all requested witness channels are being used")
- clean_strain = calibration_parts.clean_data(pipeline, clean_strain, hoftsr, witnesses, int(options.witness_channel_sr), int(options.witness_channel_sr), int(options.witness_channel_sr) / 2, 128, int(options.witness_channel_sr) * 3600, int(options.witness_channel_sr) / 2, 10, obsready = None if options.no_dq_vector else obsreadytee)
+ clean_strain = calibration_parts.clean_data(pipeline, clean_strain, hoftsr, witnesses, options.witness_channel_sr, witness_fft_samples, witness_fft_overlap, options.num_witness_ffts, witness_tf_update_samples, witness_fir_samples, options.witness_frequency_resolution, obsready = None if options.no_dq_vector else obsreadytee, attack_length = -filter_settle_time * options.witness_channel_sr, filename = "transfer_functions.txt")
if options.remove_callines or options.remove_powerlines or options.witness_channel_list is not None:
clean_strain = pipeparts.mkprogressreport(pipeline, clean_strain, "progress_hoft_cleaned_%s" % instrument)
diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py
index 989812163e..edf9c46d21 100644
--- a/gstlal-calibration/python/calibration_parts.py
+++ b/gstlal-calibration/python/calibration_parts.py
@@ -754,7 +754,7 @@ def update_filters(filter_maker, arg, filter_taker, maker_prop_name, taker_prop_
firfilter = filter_maker.get_property(maker_prop_name)[filter_number][::-1]
filter_taker.set_property(taker_prop_name, firfilter)
-def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, update_samples, fir_length, frequency_resolution, obsready = None, filename = None):
+def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, update_samples, fir_length, frequency_resolution, obsready = None, attack_length = 0, filename = None):
#
# Use witness channels that monitor the environment to remove environmental noise
@@ -775,7 +775,7 @@ def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_lengt
resampled_signal = mkresample(pipeline, signal_tee, 5, False, witness_rate)
transfer_functions = mkinterleave(pipeline, numpy.insert(witness_tees, 0, resampled_signal, axis = 0))
if obsready is not None:
- transfer_functions = mkgate(pipeline, transfer_functions, obsready, 1)
+ transfer_functions = mkgate(pipeline, transfer_functions, obsready, 1, attack_length = attack_length)
transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, update_samples = update_samples, make_fir_filters = -1, fir_length = fir_length, frequency_resolution = frequency_resolution, high_pass = 9, update_after_gap = True, filename = filename)
signal_minus_noise = [signal_tee]
for i in range(0, len(witnesses)):
diff --git a/gstlal-calibration/tests/lal_gate_test.py b/gstlal-calibration/tests/lal_gate_test.py
new file mode 100755
index 0000000000..215834ecf1
--- /dev/null
+++ b/gstlal-calibration/tests/lal_gate_test.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+# Copyright (C) 2016 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
+import test_common
+
+
+#
+# =============================================================================
+#
+# Pipelines
+#
+# =============================================================================
+#
+
+
+def lal_gate_01(pipeline, name):
+ #
+ # This pipeline tests the behavior of lal_gate with respect to attack_length, start-of-stream, etc.
+ #
+
+ rate = 512 # Hz
+ buffer_length = 1.0 # seconds
+ test_duration = 100 # seconds
+ frequency = 0.1 # Hz
+ attack_length = -3 # seconds
+ hold_length = 0 # seconds
+ threshold = 1.0
+ DC_offset = 1.0
+
+ #
+ # build pipeline
+ #
+
+ # Make a sine wave
+ src = test_common.test_src(pipeline, buffer_length = buffer_length, rate = rate, test_duration = test_duration, wave = 0, volume = 1.0, freq = frequency, width = 64)
+ # Add a DC offset
+ head = pipeparts.mkgeneric(pipeline, src, "lal_add_constant", value = DC_offset)
+ head = pipeparts.mktee(pipeline, head)
+ pipeparts.mknxydumpsink(pipeline, head, "%s_in.txt" % name)
+
+ # Gate it
+ head = calibration_parts.mkgate(pipeline, head, head, threshold, attack_length = int(attack_length * rate), hold_length = int(hold_length * rate))
+ pipeparts.mknxydumpsink(pipeline, head, "%s_out.txt" % name)
+
+ #
+ # done
+ #
+
+ return pipeline
+
+#
+# =============================================================================
+#
+# Main
+#
+# =============================================================================
+#
+
+
+test_common.build_and_run(lal_gate_01, "lal_gate_01")
+
--
GitLab