From 318562006c942478fe1320351bc1ca12e8e747ae Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Wed, 18 Jul 2018 15:12:01 -0700
Subject: [PATCH] gstlal_compute_strain: added ability to process witness
channels separately for noise subtraction, if desired.
---
gstlal-calibration/bin/gstlal_compute_strain | 29 +++++++++++--------
.../python/calibration_parts.py | 4 +++
.../tests/lal_transferfunction_test.py | 13 +++++----
3 files changed, 28 insertions(+), 18 deletions(-)
diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain
index ce52aed05c..b8c5bf484e 100755
--- a/gstlal-calibration/bin/gstlal_compute_strain
+++ b/gstlal-calibration/bin/gstlal_compute_strain
@@ -279,7 +279,7 @@ parser.add_option("--remove-powerlines", action = "store_true", help = "Set this
parser.add_option("--powerlines-channel-name", metavar = "name", default = "PEM-EY_MAINSMON_EBAY_1_DQ", help = "Set the name of the channel used as input for 60 Hz power lines to be removed. (Default = PEM-EY_MAINSMON_EBAY_1_DQ)")
parser.add_option("--powerlines-channel-sr", metavar = "Hz", type = int, default = 1024, help = "Sample rate of the powerlines channel (Default = 1024 Hz)")
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-list", metavar = "name", default = None, help = "List of witness channels to use to subtract noise from h(t), given as a string. Use semicolons (;) to separate channels/groups of channels that should be processed separately. Use commas to separate channels that should be processed together. For example, --witness-channel-list=\"chan1,chan2,chan3;chan4,chan5\". This means the pipeline will process the first 3 channels together and the last 2 channels together. (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)")
@@ -692,10 +692,14 @@ if options.remove_powerlines:
# If we are using witness channels to clean h(t), add those to the channel list
if options.witness_channel_list is not None:
- witness_channel_list = options.witness_channel_list.split(',')
- for i in range(0, len(witness_channel_list)):
- channel_list.append((instrument, witness_channel_list[i]))
- headkeys.append(witness_channel_list[i])
+ witness_channel_list = []
+ split_witness_list = options.witness_channel_list.split(';')
+ for i in range(0, len(split_witness_list)):
+ witness_channels = split_witness_list[i].split(',')
+ witness_channel_list.append(witness_channels)
+ for j in range(0, len(witness_channels)):
+ channel_list.append((instrument, witness_channels[j]))
+ headkeys.append(witness_channels[j])
####################################################################################################
@@ -1886,13 +1890,14 @@ if options.witness_channel_list is not None:
clean_strain = straintee
if options.no_dq_vector:
obsreadytee = None
- witnesses = []
- for key in headkeys:
- if key in witness_channel_list:
- 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, 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")
+ for i in range(0, len(witness_channel_list)):
+ witnesses = []
+ for key in headkeys:
+ if key in witness_channel_list[i]:
+ 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[i]):
+ print("WARNING: Not all requested witness channels are being used")
+ 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 edf9c46d21..c979058dff 100644
--- a/gstlal-calibration/python/calibration_parts.py
+++ b/gstlal-calibration/python/calibration_parts.py
@@ -252,6 +252,8 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics,
# Find transfer function between witness channel and signal at this frequency
tf_at_f = complex_division(pipeline, line_in_signal, line_in_witness)
+ # It may be necessary to remove the first few samples since denominator samples may have arrived before numerator samples, in which case the adder assumes the numerator is one.
+ tf_at_f = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_insertgap", chop_length = 1000000000) # Removing one second
# Remove worthless data from computation of transfer function if we can
if obsready is not None:
@@ -263,6 +265,8 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics,
reconstructed_line_in_signal = mkmultiplier(pipeline, list_srcs(pipeline, tf_at_f, line_in_witness))
else:
reconstructed_line_in_signal = mkmultiplier(pipeline, list_srcs(pipeline, tf_at_f, line_in_witness, phase_factor))
+ # It may be necessary to remove the first few samples since line_in_witness may have arrived first, in which case the result of the above multiplication would be line_in_witness
+ reconstructed_line_in_signal = pipeparts.mkgeneric(pipeline, reconstructed_line_in_signal, "lal_insertgap", chop_length = 1000000000)
reconstructed_line_in_signal = mkresample(pipeline, reconstructed_line_in_signal, upsample_quality, zero_latency, rate_out)
reconstructed_line_in_signal = pipeparts.mkgeneric(pipeline, reconstructed_line_in_signal, "lal_demodulate", line_frequency = -1.0 * n * f0, prefactor_real = -2.0)
reconstructed_line_in_signal = pipeparts.mkgeneric(pipeline, reconstructed_line_in_signal, "creal")
diff --git a/gstlal-calibration/tests/lal_transferfunction_test.py b/gstlal-calibration/tests/lal_transferfunction_test.py
index 43b7672910..31bd49b6aa 100755
--- a/gstlal-calibration/tests/lal_transferfunction_test.py
+++ b/gstlal-calibration/tests/lal_transferfunction_test.py
@@ -142,20 +142,21 @@ def lal_transferfunction_03(pipeline, name):
buffer_length = 1.0 # seconds
test_duration = 50.0 # seconds
width = 64 # bits
- channels = 1
- freq = 512
+ freq = 512 # Hz
+ num_witnesses = 2
+ witness_factor = 0.5 # Ratio of amplitudes of witness / signal
#
# build pipeline
#
- hoft = test_common.test_src(pipeline, buffer_length = buffer_length, wave = 5, volume = 1, freq = freq, channels = channels, rate = rate, test_duration = test_duration, width = width, verbose = False)
+ hoft = test_common.test_src(pipeline, buffer_length = buffer_length, wave = 5, volume = 1, freq = freq, channels = 1, rate = rate, test_duration = test_duration, width = width, verbose = False)
hoft = pipeparts.mktee(pipeline, hoft)
noise = []
- for i in range(0, 10):
- difference = test_common.test_src(pipeline, buffer_length = buffer_length, wave = 5, volume = 0.001, channels = channels, rate = rate, test_duration = test_duration, width = width, verbose = False)
- noisy_hoft = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, hoft, difference))
+ for i in range(0, num_witnesses):
+ difference = test_common.test_src(pipeline, buffer_length = buffer_length, wave = 5, volume = 0.001, channels = 1, rate = rate, test_duration = test_duration, width = width, verbose = False)
+ noisy_hoft = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, pipeparts.mkaudioamplify(pipeline, hoft, witness_factor), difference))
noisy_hoft = pipeparts.mkprogressreport(pipeline, noisy_hoft, "noisy_hoft_%d" % i)
noise.append(noisy_hoft)
--
GitLab