From 1934c8f3b280e0fef24534aa0b22e8bf6d00ab38 Mon Sep 17 00:00:00 2001 From: Aaron Viets <aaron.viets@ligo.org> Date: Wed, 23 May 2018 14:12:18 -0700 Subject: [PATCH] calibration_parts.py: update to 60 Hz line removal function --- .../python/calibration_parts.py | 71 +++++++++---------- 1 file changed, 33 insertions(+), 38 deletions(-) diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py index 40cae75f40..7c0e4a4baf 100644 --- a/gstlal-calibration/python/calibration_parts.py +++ b/gstlal-calibration/python/calibration_parts.py @@ -203,52 +203,47 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, zero_latency = filter_latency == 0.0 witness = pipeparts.mktee(pipeline, witness) - pipeparts.mknxydumpsink(pipeline, witness, "witness.txt") - f0_measured = pipeparts.mkgeneric(pipeline, witness, "lal_trackfrequency", num_halfcycles = 1024) - pipeparts.mknxydumpsink(pipeline, f0_measured, "f0_measured.txt") signal = pipeparts.mktee(pipeline, signal) - signal_at_60 = bandpass(pipeline, signal, 16384, f_low = 58, f_high = 62) - pipeparts.mknxydumpsink(pipeline, signal_at_60, "signal_at_60.txt") - signal_minus_lines = [signal] - - # Find amplitude and phase of first harmonic in witness channel - line_in_witness0 = pipeparts.mkgeneric(pipeline, witness, "lal_demodulate", line_frequency = f0) - line_in_witness0 = mkresample(pipeline, line_in_witness0, downsample_quality, zero_latency, compute_rate) - line_in_witness0 = lowpass(pipeline, line_in_witness0, compute_rate, length = filter_param / f0_var, fcut = 0, filter_latency = filter_latency) - line_in_witness0 = pipeparts.mktee(pipeline, line_in_witness0) # If f0 strays from its nominal value and there is a timestamp shift in the signal - # (e.g., to achieve zero latency), we need to correct the phase in the reconstructed - # signal. Start by finding the beat frequency between f0 and the actual frequency. - f0_beat_frequency = pipeparts.mkgeneric(pipeline, line_in_witness0, "lal_trackfrequency", num_halfcycles = 10) - f0_beat_frequency = pipeparts.mktee(pipeline, f0_beat_frequency) + # (e.g., to achieve zero latency), we need to correct the phase in the reconstructed + # signal. To do this, we measure the frequency in the witness and find the beat + # frequency between that and the nominal frequency f0. + if filter_latency != 0.5: + # The low-pass and resampling filters are not centered in time + f0_measured = pipeparts.mkgeneric(pipeline, witness, "lal_trackfrequency", num_halfcycles = int(2.0 * f0)) + f0_measured = mkresample(pipeline, f0_measured, 3, zero_latency, compute_rate) + f0_measured = pipeparts.mkgeneric(pipeline, f0_measured, "splitcounter") + f0_measured = pipeparts.mkgeneric(pipeline, f0_measured, "lal_smoothkappas", array_size = 1, avg_array_size = compute_rate, default_kappa_re = 0, default_to_median = True, filter_latency = filter_latency) + f0_beat_frequency = pipeparts.mkgeneric(pipeline, f0_measured, "lal_add_constant", value = -f0) + f0_beat_frequency = pipeparts.mktee(pipeline, f0_beat_frequency) - for i in range(1, num_harmonics + 1): + signal_minus_lines = [signal] + + for n in range(1, num_harmonics + 1): # Length of low-pass filter - filter_length = filter_param / (f0_var * i) + filter_length = filter_param / (f0_var * n) filter_samples = int(filter_length * compute_rate) + (1 - int(filter_length * compute_rate) % 2) sample_shift = filter_samples / 2 - int((filter_samples - 1) * filter_latency + 0.5) # shift of timestamp relative to data time_shift = float(sample_shift) / compute_rate + zero_latency * resample_shift / compute_rate - phase_angle = -2 * i * numpy.pi * time_shift - print "phase_angle = %f" % phase_angle + two_n_pi_delta_t = 2 * n * numpy.pi * time_shift - # Find phase shift due to timestamp shift for each harmonic - phase_shift = pipeparts.mkmatrixmixer(pipeline, f0_beat_frequency, matrix=[[0, phase_angle]]) - phase_shift = pipeparts.mktogglecomplex(pipeline, phase_shift) - phase_factor = pipeparts.mkgeneric(pipeline, phase_shift, "cexp") + # Only do this if we have to + if filter_latency != 0.5: + # Find phase shift due to timestamp shift for each harmonic + phase_shift = pipeparts.mkmatrixmixer(pipeline, f0_beat_frequency, matrix=[[0, two_n_pi_delta_t]]) + phase_shift = pipeparts.mktogglecomplex(pipeline, phase_shift) + phase_factor = pipeparts.mkgeneric(pipeline, phase_shift, "cexp") - if i == 1: - line_in_witness = line_in_witness0 - else: - # Find amplitude and phase of higher harmonics in witness channel - line_in_witness = pipeparts.mkgeneric(pipeline, witness, "lal_demodulate", line_frequency = i * f0) - line_in_witness = mkresample(pipeline, line_in_witness, downsample_quality, zero_latency, compute_rate) - line_in_witness = lowpass(pipeline, line_in_witness, compute_rate, length = filter_length, fcut = 0, filter_latency = filter_latency) - line_in_witness = pipeparts.mktee(pipeline, line_in_witness) + # Find amplitude and phase of each harmonic in the witness channel + line_in_witness = pipeparts.mkgeneric(pipeline, witness, "lal_demodulate", line_frequency = n * f0) + line_in_witness = mkresample(pipeline, line_in_witness, downsample_quality, zero_latency, compute_rate) + line_in_witness = lowpass(pipeline, line_in_witness, compute_rate, length = filter_length, fcut = 0, filter_latency = filter_latency) + line_in_witness = pipeparts.mktee(pipeline, line_in_witness) # Find amplitude and phase of line in signal - line_in_signal = pipeparts.mkgeneric(pipeline, signal, "lal_demodulate", line_frequency = i * f0) + line_in_signal = pipeparts.mkgeneric(pipeline, signal, "lal_demodulate", line_frequency = n * f0) line_in_signal = mkresample(pipeline, line_in_signal, downsample_quality, zero_latency, compute_rate) line_in_signal = lowpass(pipeline, line_in_signal, compute_rate, length = filter_length, fcut = 0, filter_latency = filter_latency) @@ -261,14 +256,14 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, tf_at_f = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0, array_size = 1, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency) # Use gated, averaged transfer function to reconstruct the sinusoid as it appears in the signal from the witness channel - reconstructed_line_in_signal = mkmultiplier(pipeline, list_srcs(pipeline, tf_at_f, line_in_witness, phase_factor)) + if filter_latency == 0.5: + 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)) 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 * i * f0, prefactor_real = -2.0) + 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, imag = split_into_real(pipeline, reconstructed_line_in_signal) pipeparts.mkfakesink(pipeline, imag) - if i == 1: - reconstructed_line_in_signal = pipeparts.mktee(pipeline, reconstructed_line_in_signal) - pipeparts.mknxydumpsink(pipeline, reconstructed_line_in_signal, "reconstructed_line_in_signal.txt") signal_minus_lines.append(reconstructed_line_in_signal) -- GitLab