From bc0a17c7fbdf99fd5dabb620cf344d0323d5aaf1 Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Wed, 27 Oct 2021 12:49:23 -0700
Subject: [PATCH] lal_transferfunction: Added ability to take sparse TFs and
shift frequencies
---
.../gst/lal/gstlal_transferfunction.c | 282 +++++++++++-------
.../gst/lal/gstlal_transferfunction.h | 2 +
gstlal-calibration/lib/gstlal_firtools.c | 9 +-
.../tests/check_calibration/Makefile | 13 +-
.../check_calibration/pcal2darm_timeseries.py | 10 +-
.../plot_transfer_function.py | 2 +-
6 files changed, 195 insertions(+), 123 deletions(-)
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.c b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
index 772c4f32be..ff11416a28 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.c
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
@@ -227,7 +227,7 @@ DEFINE_MAXIMUM(32);
DEFINE_MAXIMUM(64);
-static void write_transfer_functions(complex double *tfs, char *element_name, double df, gint64 rows, int columns, double t_start, double t_finish, gboolean write_to_screen, char *filename, gboolean free_name) {
+static void write_transfer_functions(complex double *tfs, char *element_name, double f0, double df, gint64 rows, int columns, double t_start, double t_finish, gboolean write_to_screen, char *filename, gboolean free_name) {
gint64 i;
int j, j_stop;
if(write_to_screen) {
@@ -238,17 +238,17 @@ static void write_transfer_functions(complex double *tfs, char *element_name, do
j_stop = columns - 1;
for(i = 0; i < rows; i++) {
- g_print("%10f\t", i * df);
+ g_print("%10f\t", f0 + i * df);
for(j = 0; j < j_stop; j++) {
if(cimag(tfs[i + j * rows]) < 0.0)
- g_print("%10e - %10ei\t\t", creal(tfs[i + j * rows]), -cimag(tfs[i + j * rows]));
+ g_print("%1.15e - %1.15ei\t\t", creal(tfs[i + j * rows]), -cimag(tfs[i + j * rows]));
else
- g_print("%10e + %10ei\t\t", creal(tfs[i + j * rows]), cimag(tfs[i + j * rows]));
+ g_print("%1.15e + %1.15ei\t\t", creal(tfs[i + j * rows]), cimag(tfs[i + j * rows]));
}
if(cimag(tfs[i + j_stop * rows]) < 0.0)
- g_print("%10e - %10ei\n", creal(tfs[i + j_stop * rows]), -cimag(tfs[i + j_stop * rows]));
+ g_print("%1.15e - %1.15ei\n", creal(tfs[i + j_stop * rows]), -cimag(tfs[i + j_stop * rows]));
else
- g_print("%10e + %10ei\n", creal(tfs[i + j_stop * rows]), cimag(tfs[i + j_stop * rows]));
+ g_print("%1.15e + %1.15ei\n", creal(tfs[i + j_stop * rows]), cimag(tfs[i + j_stop * rows]));
}
g_print("\n\n");
}
@@ -263,17 +263,17 @@ static void write_transfer_functions(complex double *tfs, char *element_name, do
j_stop = columns - 1;
for(i = 0; i < rows; i++) {
- g_fprintf(fp, "%10f\t", i * df);
+ g_fprintf(fp, "%10f\t", f0 + i * df);
for(j = 0; j < j_stop; j++) {
if(cimag(tfs[i + j * rows]) < 0.0)
- g_fprintf(fp, "%10e - %10ei\t\t", creal(tfs[i + j * rows]), -cimag(tfs[i + j * rows]));
+ g_fprintf(fp, "%1.15e - %1.15ei\t\t", creal(tfs[i + j * rows]), -cimag(tfs[i + j * rows]));
else
- g_fprintf(fp, "%10e + %10ei\t\t", creal(tfs[i + j * rows]), cimag(tfs[i + j * rows]));
+ g_fprintf(fp, "%1.15e + %1.15ei\t\t", creal(tfs[i + j * rows]), cimag(tfs[i + j * rows]));
}
if(cimag(tfs[i + j_stop * rows]) < 0.0)
- g_fprintf(fp, "%10e - %10ei\n", creal(tfs[i + j_stop * rows]), -cimag(tfs[i + j_stop * rows]));
+ g_fprintf(fp, "%1.15e - %1.15ei\n", creal(tfs[i + j_stop * rows]), -cimag(tfs[i + j_stop * rows]));
else
- g_fprintf(fp, "%10e + %10ei\n", creal(tfs[i + j_stop * rows]), cimag(tfs[i + j_stop * rows]));
+ g_fprintf(fp, "%1.15e + %1.15ei\n", creal(tfs[i + j_stop * rows]), cimag(tfs[i + j_stop * rows]));
}
g_fprintf(fp, "\n\n");
fclose(fp);
@@ -295,8 +295,8 @@ static void write_fir_filters(double *filters, char *element_name, gint64 rows,
j_stop = columns - 1;
for(i = 0; i < rows; i++) {
for(j = 0; j < j_stop; j++)
- g_print("%10e\t", filters[i + j * rows]);
- g_print("%10e\n", filters[i + j_stop * rows]);
+ g_print("%1.15e\t", filters[i + j * rows]);
+ g_print("%1.15e\n", filters[i + j_stop * rows]);
}
g_print("\n\n");
}
@@ -312,8 +312,8 @@ static void write_fir_filters(double *filters, char *element_name, gint64 rows,
j_stop = columns - 1;
for(i = 0; i < rows; i++) {
for(j = 0; j < j_stop; j++)
- g_fprintf(fp, "%10e\t", filters[i + j * rows]);
- g_fprintf(fp, "%10e\n", filters[i + j_stop * rows]);
+ g_fprintf(fp, "%1.15e\t", filters[i + j * rows]);
+ g_fprintf(fp, "%1.15e\n", filters[i + j_stop * rows]);
}
g_fprintf(fp, "\n\n");
fclose(fp);
@@ -402,9 +402,9 @@ static gboolean update_transfer_functions_ ## DTYPE(complex DTYPE *autocorrelati
* First, apply the sinc filter to higher frequencies. We could hit the edge
* of the sinc table or the Nyquist frequency of the transfer function.
*/ \
- sinc0 = (sinc_taps_per_input - i * sinc_taps_per_output % sinc_taps_per_input) % sinc_taps_per_input; \
+ sinc0 = sinc_length / 2 + (sinc_taps_per_input - i * sinc_taps_per_output % sinc_taps_per_input) % sinc_taps_per_input; \
input0 = j + (i * (fd_fft_length - 1) + fd_tf_length - 2) / (fd_tf_length - 1) * elements_per_freq; \
- k_stop = minimum64((sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input, fd_fft_length - (i * (fd_fft_length - 1) + fd_tf_length - 2) / (fd_tf_length - 1)); \
+ k_stop = minimum64((sinc_taps_per_input - 1 + sinc_length - sinc0) / sinc_taps_per_input, fd_fft_length - (i * (fd_fft_length - 1) + fd_tf_length - 2) / (fd_tf_length - 1)); \
for(k = 0; k < k_stop; k++) \
z += sinc_table[sinc0 + k * sinc_taps_per_input] * autocorrelation_matrix[input0 + k * elements_per_freq]; \
/*
@@ -413,27 +413,27 @@ static gboolean update_transfer_functions_ ## DTYPE(complex DTYPE *autocorrelati
*/ \
sinc0 += k_stop * sinc_taps_per_input; \
input0 += (k_stop - 2) * elements_per_freq; \
- k_stop = (sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input; \
+ k_stop = (sinc_taps_per_input - 1 + sinc_length - sinc0) / sinc_taps_per_input; \
for(k = 0; k < k_stop; k++) \
z += sinc_table[sinc0 + k * sinc_taps_per_input] * conj ## F_OR_BLANK(autocorrelation_matrix[input0 - k * elements_per_freq]); \
/*
* Now, go back and apply the sinc filter to the lower frequencies. We could hit the edge
* of the sinc table or the DC component of the transfer function.
*/ \
- sinc0 = 1 + (sinc_taps_per_input + i * sinc_taps_per_output - 1) % sinc_taps_per_input; \
+ sinc0 = sinc_length / 2 + (sinc_taps_per_input - i * sinc_taps_per_output % sinc_taps_per_input) % sinc_taps_per_input - sinc_taps_per_input; \
input0 = j + (i * (fd_fft_length - 1) - 1) / (fd_tf_length - 1) * elements_per_freq; \
- k_stop = minimum64((sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input, (i * (fd_fft_length - 1) - 1) / (fd_tf_length - 1)); \
+ k_stop = minimum64((sinc0 + sinc_taps_per_input) / sinc_taps_per_input, (i * (fd_fft_length - 1) - 1) / (fd_tf_length - 1)); \
for(k = 0; k < k_stop; k++) \
- z += sinc_table[sinc0 + k * sinc_taps_per_input] * autocorrelation_matrix[input0 - k * elements_per_freq]; \
+ z += sinc_table[sinc0 - k * sinc_taps_per_input] * autocorrelation_matrix[input0 - k * elements_per_freq]; \
/*
* If we hit the DC component of the transfer function but not the edge of the
* sinc table, turn around and keep going until we hit the edge of the sinc table.
*/ \
- sinc0 += k_stop * sinc_taps_per_input; \
+ sinc0 -= k_stop * sinc_taps_per_input; \
input0 = j; \
- k_stop = (sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input; \
+ k_stop = (sinc0 + sinc_taps_per_input) / sinc_taps_per_input; \
for(k = 0; k < k_stop; k++) \
- z += sinc_table[sinc0 + k * sinc_taps_per_input] * conj ## F_OR_BLANK(autocorrelation_matrix[input0 + k * elements_per_freq]); \
+ z += sinc_table[sinc0 - k * sinc_taps_per_input] * conj ## F_OR_BLANK(autocorrelation_matrix[input0 + k * elements_per_freq]); \
\
/* Set an element of the GSL vector transfer_functions_at_f */ \
gsl_vector_complex_set(transfer_functions_at_f, j, gsl_complex_rect(creal((complex double) z) / num_avg, cimag((complex double) z) / num_avg)); \
@@ -447,9 +447,9 @@ static gboolean update_transfer_functions_ ## DTYPE(complex DTYPE *autocorrelati
* First, apply the sinc filter to higher frequencies. We could hit the edge
* of the sinc table or the Nyquist frequency of the transfer function.
*/ \
- sinc0 = (sinc_taps_per_input - i * sinc_taps_per_output % sinc_taps_per_input) % sinc_taps_per_input; \
+ sinc0 = sinc_length / 2 + (sinc_taps_per_input - i * sinc_taps_per_output % sinc_taps_per_input) % sinc_taps_per_input; \
input0 = j + num_tfs + (i * (fd_fft_length - 1) + fd_tf_length - 2) / (fd_tf_length - 1) * elements_per_freq; \
- k_stop = minimum64((sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input, fd_fft_length - (i * (fd_fft_length - 1) + fd_tf_length - 2) / (fd_tf_length - 1)); \
+ k_stop = minimum64((sinc_taps_per_input - 1 + sinc_length - sinc0) / sinc_taps_per_input, fd_fft_length - (i * (fd_fft_length - 1) + fd_tf_length - 2) / (fd_tf_length - 1)); \
for(k = 0; k < k_stop; k++) \
z += sinc_table[sinc0 + k * sinc_taps_per_input] * autocorrelation_matrix[input0 + k * elements_per_freq]; \
/*
@@ -458,27 +458,27 @@ static gboolean update_transfer_functions_ ## DTYPE(complex DTYPE *autocorrelati
*/ \
sinc0 += k_stop * sinc_taps_per_input; \
input0 += (k_stop - 2) * elements_per_freq; \
- k_stop = (sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input; \
+ k_stop = (sinc_taps_per_input - 1 + sinc_length - sinc0) / sinc_taps_per_input; \
for(k = 0; k < k_stop; k++) \
z += sinc_table[sinc0 + k * sinc_taps_per_input] * conj ## F_OR_BLANK(autocorrelation_matrix[input0 - k * elements_per_freq]); \
/*
* Now, go back and apply the sinc filter to the lower frequencies. We could hit the edge
* of the sinc table or the DC component of the transfer function.
*/ \
- sinc0 = 1 + (sinc_taps_per_input + i * sinc_taps_per_output - 1) % sinc_taps_per_input; \
+ sinc0 = sinc_length / 2 + (sinc_taps_per_input - i * sinc_taps_per_output % sinc_taps_per_input) % sinc_taps_per_input - sinc_taps_per_input; \
input0 = j + num_tfs + (i * (fd_fft_length - 1) - 1) / (fd_tf_length - 1) * elements_per_freq; \
- k_stop = minimum64((sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input, (i * (fd_fft_length - 1) - 1) / (fd_tf_length - 1)); \
+ k_stop = minimum64((sinc0 + sinc_taps_per_input) / sinc_taps_per_input, (i * (fd_fft_length - 1) - 1) / (fd_tf_length - 1)); \
for(k = 0; k < k_stop; k++) \
- z += sinc_table[sinc0 + k * sinc_taps_per_input] * autocorrelation_matrix[input0 - k * elements_per_freq]; \
+ z += sinc_table[sinc0 - k * sinc_taps_per_input] * autocorrelation_matrix[input0 - k * elements_per_freq]; \
/*
* If we hit the DC component of the transfer function but not the edge of the
* sinc table, turn around and keep going until we hit the edge of the sinc table.
*/ \
- sinc0 += k_stop * sinc_taps_per_input; \
+ sinc0 -= k_stop * sinc_taps_per_input; \
input0 = j + num_tfs; \
- k_stop = (sinc_taps_per_input + sinc_length / 2 - sinc0) / sinc_taps_per_input; \
+ k_stop = (sinc0 + sinc_taps_per_input) / sinc_taps_per_input; \
for(k = 0; k < k_stop; k++) \
- z += sinc_table[sinc0 + k * sinc_taps_per_input] * conj ## F_OR_BLANK(autocorrelation_matrix[input0 + k * elements_per_freq]); \
+ z += sinc_table[sinc0 - k * sinc_taps_per_input] * conj ## F_OR_BLANK(autocorrelation_matrix[input0 + k * elements_per_freq]); \
\
/* Set an element of the GSL matrix autocorrelation_matrix_at_f */ \
gsl_matrix_complex_set(autocorrelation_matrix_at_f, j / num_tfs, j % num_tfs, gsl_complex_rect(creal((complex double) z) / num_avg, cimag((complex double) z) / num_avg)); \
@@ -566,7 +566,10 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
src_size /= element->unit_size; \
\
gint64 i, j, k, m, num_ffts, num_ffts_in_avg_if_nogap, k_start, k_stop, first_index, first_index2, fd_fft_length, fd_tf_length, stride, num_tfs, notch_start, notch_end; \
- fd_fft_length = element->fft_length / 2 + 1; \
+ if(element->df <= (double) element->rate / element->fft_length) \
+ fd_fft_length = element->fft_length / 2 + 1; \
+ else \
+ fd_fft_length = (gint64) (element->b / 2 / element->a + 1); \
fd_tf_length = element->fir_length / 2 + 1; \
stride = element->fft_length - element->fft_overlap; \
num_tfs = element->channels - 1; \
@@ -600,7 +603,17 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
real_fft[k] = element->workspace.w ## S_OR_D ## pf.fft_window[k] * src[j + element->channels * (k - k_start)]; \
\
/* Take an FFT */ \
- if(element->use_fir_fft) { \
+ if(element->df > (double) element->rate / element->fft_length) { \
+ /* Use a sparse FFT algorithm to process a "comb" */ \
+ firfft = gstlal_rfft_comb_ ## DTYPE(real_fft, (guint) element->fft_length, element->a, element->b, (DTYPE) element->f0, FALSE, NULL); \
+ /* Copy FFT to the proper location */ \
+ first_index = j * fd_fft_length; \
+ for(k = 0; k < fd_fft_length; k++) \
+ element->workspace.w ## S_OR_D ## pf.ffts[first_index + k] = firfft[k]; \
+ \
+ g_free(firfft); \
+ \
+ } else if(element->use_fir_fft) { \
firfft = gstlal_rfft_ ## DTYPE(real_fft, (guint) element->fft_length, NULL, 0, NULL, FALSE, 0, NULL, 0, NULL, NULL, NULL, FALSE); \
\
/* Copy FFT to the proper location */ \
@@ -735,7 +748,17 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
real_fft[k] = element->workspace.w ## S_OR_D ## pf.fft_window[k] * ptr[first_index + k * element->channels]; \
\
/* Take an FFT */ \
- if(element->use_fir_fft) { \
+ if(element->df > (double) element->rate / element->fft_length) { \
+ /* Use a sparse FFT algorithm to process a "comb" */ \
+ firfft = gstlal_rfft_comb_ ## DTYPE(real_fft, (guint) element->fft_length, element->a, element->b, (DTYPE) element->f0, FALSE, NULL); \
+ /* Copy FFT to the proper location */ \
+ first_index = j * fd_fft_length; \
+ for(k = 0; k < fd_fft_length; k++) \
+ element->workspace.w ## S_OR_D ## pf.ffts[first_index + k] = firfft[k]; \
+ \
+ g_free(firfft); \
+ \
+ } else if(element->use_fir_fft) { \
firfft = gstlal_rfft_ ## DTYPE(real_fft, (guint) element->fft_length, NULL, 0, NULL, FALSE, 0, NULL, 0, NULL, NULL, NULL, FALSE); \
\
/* Copy FFT to the proper location */ \
@@ -927,8 +950,11 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
/* Let other elements know about the update */ \
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]); \
/* Write transfer functions to the screen or a file if we want */ \
- if(element->write_to_screen || element->filename) \
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (num_ffts_in_avg_if_nogap * stride + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE); \
+ if(element->write_to_screen || element->filename) { \
+ gint64 rows = fd_tf_length; \
+ double df = element->rate / 2.0 / (fd_tf_length - 1.0); \
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->make_fir_filters ? 0 : element->f0 * element->df, df, rows, num_tfs, element->t_start_tf, element->t_start_tf + (double) (num_ffts_in_avg_if_nogap * stride + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE); \
+ } \
/* If this is this first transfer function after a gap, we may wish to store it */ \
if(element->use_first_after_gap && !element->num_tfs_since_gap) { \
for(i = 0; i < num_tfs * fd_tf_length; i++) \
@@ -1031,8 +1057,8 @@ static gboolean start(GstBaseSink *sink) {
element->fft_overlap = element->fft_length - 1;
}
if(element->fir_length % 2) {
- GST_WARNING_OBJECT(element, "The chosen fir-length must be even. Adding 1 to make it even.");
- element->fir_length += 1;
+ GST_WARNING_OBJECT(element, "The chosen fir-length must be even. Removing 1 sample to make it even.");
+ element->fir_length -= 1;
}
if(element->make_fir_filters && element->frequency_resolution < (double) element->rate / element->fir_length)
GST_WARNING_OBJECT(element, "The specified frequency resolution is finer than 1/fir_length, which cannot be achieved. The actual frequency resolution will be reset to 1/MIN(fft_length, fir_length).");
@@ -1064,9 +1090,13 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
if(element->workspace.wspf.num_ffts_in_avg - element->workspace.wspf.num_ffts_dropped > element->min_ffts && !element->computed_full_tfs) {
/* Then we should use the transfer functions we have now, since they won't get any better */
- gint64 fd_fft_length = element->fft_length / 2 + 1;
- gint64 fd_tf_length = element->fir_length / 2 + 1;
- gint64 num_tfs = element->channels - 1;
+ gint64 fd_fft_length, fd_tf_length, num_tfs;
+ if(element->df <= (double) element->rate / element->fft_length)
+ fd_fft_length = element->fft_length / 2 + 1;
+ else
+ fd_fft_length = (gint64) (element->b / 2 / element->a + 1);
+ fd_tf_length = element->fir_length / 2 + 1;
+ num_tfs = element->channels - 1;
if(element->use_median) {
/* Then we still need to fill the autocorrelation matrix with median values */
gint64 median_length = element->num_ffts / 2 + 1;
@@ -1119,7 +1149,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (element->workspace.wspf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->make_fir_filters ? 0 : element->f0 * element->df, element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (element->workspace.wspf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
} else
GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. No transfer functions will be produced.");
/* Update FIR filters if we want */
@@ -1139,9 +1169,13 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
} else {
if(element->workspace.wdpf.num_ffts_in_avg - element->workspace.wdpf.num_ffts_dropped > element->min_ffts && !element->computed_full_tfs) {
/* Then we should use the transfer functions we have now, since they won't get any better */
- gint64 fd_fft_length = element->fft_length / 2 + 1;
- gint64 fd_tf_length = element->fir_length / 2 + 1;
- gint64 num_tfs = element->channels - 1;
+ gint64 fd_fft_length, fd_tf_length, num_tfs;
+ if(element->df <= (double) element->rate / element->fft_length)
+ fd_fft_length = element->fft_length / 2 + 1;
+ else
+ fd_fft_length = (gint64) (element->b / 2 / element->a + 1);
+ fd_tf_length = element->fir_length / 2 + 1;
+ num_tfs = element->channels - 1;
if(element->use_median) {
/* Then we still need to fill the autocorrelation matrix with median values */
gint64 median_length = element->num_ffts / 2 + 1;
@@ -1194,7 +1228,7 @@ static gboolean event(GstBaseSink *sink, GstEvent *event) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (element->workspace.wdpf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->make_fir_filters ? 0 : element->f0 * element->df, element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->t_start_tf, element->t_start_tf + (double) (element->workspace.wdpf.num_ffts_in_avg * (element->fft_length - element->fft_overlap) + element->fft_overlap) / element->rate, element->write_to_screen, element->filename, TRUE);
} else
GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. No transfer functions will be produced.");
/* Update FIR filters if we want */
@@ -1427,13 +1461,52 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
goto done;
}
+ /*
+ * We can compute a "sparse" transfer function if the desired frequency spacing of the
+ * FFTs is greater than 1 / fft_length. This is useful for measuring comb injections.
+ */
+
+ gint64 fd_fft_length;
+
+ if(element->df <= (double) element->rate / element->fft_length) {
+ element->df = (double) element->rate / element->fft_length;
+ fd_fft_length = element->fft_length / 2 + 1;
+ element->f0 = 0.0;
+ } else {
+ /* Find integers a and b such that a/b = element->df / element->rate */
+ guint a, best_a;
+ double b, best_b, error;
+ best_a = a = 1;
+ best_b = b = a * element->rate / element->df;
+ error = fabs((int) (b + 0.5) - b);
+ while(error > 1e-10 && a < 1000) {
+ a += 1;
+ b = a * element->rate / element->df;
+ if(fabs((int) (b + 0.5) - b) < error) {
+ best_a = a;
+ best_b = b;
+ error = fabs((int) (b + 0.5) - b);
+ }
+ a++;
+ }
+ element->a = best_a;
+ element->b = (guint) best_b;
+ element->df = best_a * element->rate / best_b;
+ /* Unit conversion of f0 into units of df */
+ element->f0 /= element->df;
+ element->f0 -= (int) element->f0;
+ if(element->f0 < 0)
+ element->f0 += 1; \
+
+ fd_fft_length = (gint64) (element->b / 2 / element->a + 1); \
+ }
+
/*
* Allocate any memory that depends on stream parameters
*/
- gint64 fd_fft_length = element->fft_length / 2 + 1;
if(!element->fir_length)
- element->fir_length = element->fft_length;
+ element->fir_length = 2 * fd_fft_length - 2;
gint64 fd_fir_length = element->fir_length / 2 + 1;
element->transfer_functions = g_malloc((element->channels - 1) * fd_fir_length * sizeof(*element->transfer_functions));
memset(element->transfer_functions, 0, (element->channels - 1) * fd_fir_length * sizeof(*element->transfer_functions));
@@ -1513,7 +1586,7 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
* Make a sinc table to resample and/or low-pass filter the transfer functions when we make FIR filters
*/
- if(element->fir_length == element->fft_length) {
+ if(fd_fir_length == fd_fft_length && !(element->f0 && element->make_fir_filters)) {
element->workspace.wspf.sinc_length = 1;
element->workspace.wspf.sinc_table = g_malloc(sizeof(*element->workspace.wspf.sinc_table));
*element->workspace.wspf.sinc_table = 1.0;
@@ -1530,47 +1603,41 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
while(common_denominator % short_length)
common_denominator += long_length;
element->workspace.wspf.sinc_taps_per_df = common_denominator / long_length;
+ /* Number of taps per frequency bin of sinc table at fft frequency resolution */
+ gint64 taps_per_fftdf = common_denominator / (fd_fft_length - 1);
/* taps_per_osc is the number of taps per half-oscillation in the sinc table */
- gint64 taps_per_osc = element->workspace.wspf.sinc_taps_per_df * (gint64) (maximum(maximum(element->frequency_resolution * element->fir_length / element->rate, element->frequency_resolution * element->fft_length / element->rate), maximum((double) element->fft_length / element->fir_length, (double) element->fir_length / element->fft_length)) + 0.5);
+ gint64 taps_per_osc = element->workspace.wspf.sinc_taps_per_df * (gint64) (maximum(maximum(element->frequency_resolution * element->fir_length / element->rate, element->frequency_resolution * element->fft_length / element->rate), maximum((double) (fd_fft_length - 1) / (fd_fir_length - 1), (double) (fd_fir_length - 1) / (fd_fft_length - 1))) + 0.5);
element->workspace.wspf.sinc_length = minimum64(element->workspace.wspf.sinc_taps_per_df * maximum64(fd_fir_length / 2, fd_fft_length / 2) - 1, 1 + SINC_LENGTH * taps_per_osc);
+ element->workspace.wspf.sinc_length += (1 + element->workspace.wspf.sinc_length) % 2;
- /* To save memory, we use symmetry and record only half of the sinc table */
- element->workspace.wspf.sinc_table = g_malloc((1 + element->workspace.wspf.sinc_length / 2) * sizeof(*element->workspace.wspf.sinc_table));
- *element->workspace.wspf.sinc_table = 1.0;
+ element->workspace.wspf.sinc_table = g_malloc(element->workspace.wspf.sinc_length * sizeof(*element->workspace.wspf.sinc_table));
gint64 j;
- float sin_arg, normalization;
+ float sin_arg_up, sin_arg_down, normalization;
+ sin_arg_up = sin_arg_down = element->f0 * M_PI * taps_per_fftdf / taps_per_osc;
+ if(sin_arg_up)
+ element->workspace.wspf.sinc_table[element->workspace.wspf.sinc_length / 2] = sinf(sin_arg_up) / sin_arg_up;
+ else
+ element->workspace.wspf.sinc_table[element->workspace.wspf.sinc_length / 2] = 1.0;
for(i = 1; i <= element->workspace.wspf.sinc_length / 2; i++) {
- sin_arg = M_PI * i / taps_per_osc;
- element->workspace.wspf.sinc_table[i] = sinf(sin_arg) / sin_arg;
+ sin_arg_up += M_PI / taps_per_osc;
+ sin_arg_down -= M_PI / taps_per_osc;
+ element->workspace.wspf.sinc_table[element->workspace.wspf.sinc_length / 2 + i] = sinf(sin_arg_up) / sin_arg_up;
+ element->workspace.wspf.sinc_table[element->workspace.wspf.sinc_length / 2 - i] = sinf(sin_arg_down) / sin_arg_down;
}
/* Window the sinc table */
- kaiser_float(element->workspace.wspf.sinc_length, 10.0, element->workspace.wspf.sinc_table, TRUE);
+ kaiser_float(element->workspace.wspf.sinc_length, 10.0, element->workspace.wspf.sinc_table, FALSE);
/*
* Normalize the sinc table to make the DC gain exactly 1. We need to account for the fact
* that the density of taps in the filter could be higher than the density of input samples.
*/
gint64 taps_per_input = fd_fir_length > fd_fft_length ? common_denominator / short_length : element->workspace.wspf.sinc_taps_per_df;
- for(i = 0; i < (taps_per_input + 1) / 2; i++) {
+ for(i = 0; i < taps_per_input; i++) {
normalization = 0.0;
- for(j = i; j <= element->workspace.wspf.sinc_length / 2; j += taps_per_input)
- normalization += element->workspace.wspf.sinc_table[j];
- for(j = taps_per_input - i; j <= element->workspace.wspf.sinc_length / 2; j += taps_per_input)
+ for(j = i; j < element->workspace.wspf.sinc_length; j += taps_per_input)
normalization += element->workspace.wspf.sinc_table[j];
- for(j = i; j <= element->workspace.wspf.sinc_length / 2; j += taps_per_input)
- element->workspace.wspf.sinc_table[j] /= normalization;
- if(i) {
- for(j = taps_per_input - i; j <= element->workspace.wspf.sinc_length / 2; j += taps_per_input)
- element->workspace.wspf.sinc_table[j] /= normalization;
- }
- }
- /* If taps_per_input is even, we need to account for one more normalization without "over-normalizing." */
- if(!((taps_per_input) % 2)) {
- normalization = 0.0;
- for(j = taps_per_input / 2; j <= element->workspace.wspf.sinc_length / 2; j += taps_per_input)
- normalization += 2 * element->workspace.wspf.sinc_table[j];
- for(j = taps_per_input / 2; j <= element->workspace.wspf.sinc_length / 2; j += taps_per_input)
+ for(j = i; j < element->workspace.wspf.sinc_length; j += taps_per_input)
element->workspace.wspf.sinc_table[j] /= normalization;
}
}
@@ -1683,7 +1750,7 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
GST_LOG_OBJECT(element, "starting FFTWF planning");
/* data that will be Fourier transformed into frequency domain */
- element->workspace.wspf.fft = (complex float *) fftwf_malloc(fd_fft_length * sizeof(*element->workspace.wspf.fft));
+ element->workspace.wspf.fft = (complex float *) fftwf_malloc((element->fft_length / 2 + 1) * sizeof(*element->workspace.wspf.fft));
element->workspace.wspf.plan = fftwf_plan_dft_r2c_1d(element->fft_length, (float *) element->workspace.wspf.fft, element->workspace.wspf.fft, FFTW_ESTIMATE);
if(element->make_fir_filters && !element->workspace.wspf.fir_filter) {
@@ -1738,7 +1805,7 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
* Make a sinc table to resample and/or low-pass filter the transfer functions when we make FIR filters
*/
- if(element->fir_length == element->fft_length) {
+ if(fd_fir_length == fd_fft_length && !(element->f0 && element->make_fir_filters)) {
element->workspace.wdpf.sinc_length = 1;
element->workspace.wdpf.sinc_table = g_malloc(sizeof(*element->workspace.wdpf.sinc_table));
*element->workspace.wdpf.sinc_table = 1.0;
@@ -1755,47 +1822,41 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
while(common_denominator % short_length)
common_denominator += long_length;
element->workspace.wdpf.sinc_taps_per_df = common_denominator / long_length;
+ /* Number of taps per frequency bin of sinc table at fft frequency resolution */
+ gint64 taps_per_fftdf = common_denominator / (fd_fft_length - 1);
/* taps_per_osc is the number of taps per half-oscillation in the sinc table */
- gint64 taps_per_osc = element->workspace.wdpf.sinc_taps_per_df * (gint64) (maximum(maximum(element->frequency_resolution * element->fir_length / element->rate, element->frequency_resolution * element->fft_length / element->rate), maximum((double) element->fft_length / element->fir_length, (double) element->fir_length / element->fft_length)) + 0.5);
+ gint64 taps_per_osc = element->workspace.wdpf.sinc_taps_per_df * (gint64) (maximum(maximum(element->frequency_resolution * element->fir_length / element->rate, element->frequency_resolution * element->fft_length / element->rate), maximum((double) (fd_fft_length - 1) / (fd_fir_length - 1), (double) (fd_fir_length - 1) / (fd_fft_length - 1))) + 0.5);
element->workspace.wdpf.sinc_length = minimum64(element->workspace.wdpf.sinc_taps_per_df * maximum64(fd_fir_length / 2, fd_fft_length / 2) - 1, 1 + SINC_LENGTH * taps_per_osc);
+ element->workspace.wdpf.sinc_length += (1 + element->workspace.wdpf.sinc_length) % 2;
- /* To save memory, we use symmetry and record only half of the sinc table */
- element->workspace.wdpf.sinc_table = g_malloc((1 + element->workspace.wdpf.sinc_length / 2) * sizeof(*element->workspace.wdpf.sinc_table));
- *element->workspace.wdpf.sinc_table = 1.0;
+ element->workspace.wdpf.sinc_table = g_malloc(element->workspace.wdpf.sinc_length * sizeof(*element->workspace.wdpf.sinc_table));
gint64 j;
- double sin_arg, normalization;
+ double sin_arg_up, sin_arg_down, normalization;
+ sin_arg_up = sin_arg_down = element->f0 * M_PI * taps_per_fftdf / taps_per_osc;
+ if(sin_arg_up)
+ element->workspace.wdpf.sinc_table[element->workspace.wdpf.sinc_length / 2] = sinf(sin_arg_up) / sin_arg_up;
+ else
+ element->workspace.wdpf.sinc_table[element->workspace.wdpf.sinc_length / 2] = 1.0;
for(i = 1; i <= element->workspace.wdpf.sinc_length / 2; i++) {
- sin_arg = M_PI * i / taps_per_osc;
- element->workspace.wdpf.sinc_table[i] = sin(sin_arg) / sin_arg;
+ sin_arg_up += M_PI / taps_per_osc;
+ sin_arg_down -= M_PI / taps_per_osc;
+ element->workspace.wdpf.sinc_table[element->workspace.wdpf.sinc_length / 2 + i] = sin(sin_arg_up) / sin_arg_up;
+ element->workspace.wdpf.sinc_table[element->workspace.wdpf.sinc_length / 2 - i] = sin(sin_arg_down) / sin_arg_down;
}
/* Window the sinc table */
- kaiser_double(element->workspace.wdpf.sinc_length, 10.0, element->workspace.wdpf.sinc_table, TRUE);
+ kaiser_double(element->workspace.wdpf.sinc_length, 10.0, element->workspace.wdpf.sinc_table, FALSE);
/*
* Normalize the sinc table to make the DC gain exactly 1. We need to account for the fact
* that the density of taps in the filter could be higher than the density of input samples.
*/
gint64 taps_per_input = fd_fir_length > fd_fft_length ? common_denominator / short_length : element->workspace.wdpf.sinc_taps_per_df;
- for(i = 0; i < (taps_per_input + 1) / 2; i++) {
+ for(i = 0; i < taps_per_input; i++) {
normalization = 0.0;
- for(j = i; j <= element->workspace.wdpf.sinc_length / 2; j += taps_per_input)
- normalization += element->workspace.wdpf.sinc_table[j];
- for(j = taps_per_input - i; j <= element->workspace.wdpf.sinc_length / 2; j += taps_per_input)
+ for(j = i; j < element->workspace.wdpf.sinc_length; j += taps_per_input)
normalization += element->workspace.wdpf.sinc_table[j];
- for(j = i; j <= element->workspace.wdpf.sinc_length / 2; j += taps_per_input)
- element->workspace.wdpf.sinc_table[j] /= normalization;
- if(i) {
- for(j = taps_per_input - i; j <= element->workspace.wdpf.sinc_length / 2; j += taps_per_input)
- element->workspace.wdpf.sinc_table[j] /= normalization;
- }
- }
- /* If taps_per_input is even, we need to account for one more normalization without "over-normalizing." */
- if(!((taps_per_input) % 2)) {
- normalization = 0.0;
- for(j = taps_per_input / 2; j <= element->workspace.wdpf.sinc_length / 2; j += taps_per_input)
- normalization += 2 * element->workspace.wdpf.sinc_table[j];
- for(j = taps_per_input / 2; j <= element->workspace.wdpf.sinc_length / 2; j += taps_per_input)
+ for(j = i; j < element->workspace.wdpf.sinc_length; j += taps_per_input)
element->workspace.wdpf.sinc_table[j] /= normalization;
}
}
@@ -1908,7 +1969,7 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
GST_LOG_OBJECT(element, "starting FFTW planning");
/* data that will be Fourier transformed into frequency domain */
- element->workspace.wdpf.fft = (complex double *) fftw_malloc(fd_fft_length * sizeof(*element->workspace.wdpf.fft));
+ element->workspace.wdpf.fft = (complex double *) fftw_malloc((element->fft_length / 2 + 1) * sizeof(*element->workspace.wdpf.fft));
element->workspace.wdpf.plan = fftw_plan_dft_r2c_1d(element->fft_length, (double *) element->workspace.wdpf.fft, element->workspace.wdpf.fft, FFTW_ESTIMATE);
if(element->make_fir_filters && !element->workspace.wdpf.fir_filter) {
@@ -1995,7 +2056,7 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
/* Write transfer functions to the screen or a file if we want */
if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), (double) element->rate / element->fir_length, element->fir_length / 2 + 1, element->channels - 1, 0, 0, element->write_to_screen, element->filename, TRUE);
+ write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->make_fir_filters ? 0 : element->f0 * element->df, (double) element->rate / element->fir_length, element->fir_length / 2 + 1, element->channels - 1, 0, 0, element->write_to_screen, element->filename, TRUE);
if(element->make_fir_filters) {
for(i = 0; i < (element->channels - 1) * element->fir_length; i++)
element->fir_filters[i] = element->post_gap_fir_filters[i];
@@ -2060,11 +2121,16 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
/* Check whether we need to do anything with this data */
if(element->sample_count > element->update_samples && mapinfo.size) {
+ gint64 fd_fft_length;
+ if(element->df <= (double) element->rate / element->fft_length) \
+ fd_fft_length = element->fft_length / 2 + 1; \
+ else \
+ fd_fft_length = (gint64) (element->b / 2 / element->a + 1); \
gboolean success;
if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
/* If we are just beginning to compute new transfer functions with this data, initialize memory that we will fill to zero */
if(!element->workspace.wspf.num_ffts_in_avg) {
- memset(element->workspace.wspf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * (element->fft_length / 2 + 1) * sizeof(*element->workspace.wspf.autocorrelation_matrix));
+ memset(element->workspace.wspf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * fd_fft_length * sizeof(*element->workspace.wspf.autocorrelation_matrix));
element->t_start_tf = (double) (gst_util_uint64_scale_int_round(GST_BUFFER_PTS(buffer) + GST_BUFFER_DURATION(buffer), element->rate, GST_SECOND) - element->sample_count + element->update_samples) / element->rate;
}
/* Send the data to a function to compute fft's and transfer functions */
@@ -2072,7 +2138,7 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
} else {
/* If we are just beginning to compute new transfer functions with this data, initialize memory that we will fill to zero */
if(!element->workspace.wdpf.num_ffts_in_avg) {
- memset(element->workspace.wdpf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * (element->fft_length / 2 + 1) * sizeof(*element->workspace.wdpf.autocorrelation_matrix));
+ memset(element->workspace.wdpf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * fd_fft_length * sizeof(*element->workspace.wdpf.autocorrelation_matrix));
element->t_start_tf = (double) (gst_util_uint64_scale_int_round(GST_BUFFER_PTS(buffer) + GST_BUFFER_DURATION(buffer), element->rate, GST_SECOND) - element->sample_count + element->update_samples) / element->rate;
}
/* Send the data to a function to compute fft's and transfer functions */
@@ -3032,6 +3098,8 @@ static void gstlal_transferfunction_init(GSTLALTransferFunction *element) {
element->gap_samples = 0;
element->num_tfs_since_gap = 0;
element->computed_full_tfs = FALSE;
+ element->a = 0;
+ element->b = 0;
gst_base_sink_set_sync(GST_BASE_SINK(element), FALSE);
gst_base_sink_set_async_enabled(GST_BASE_SINK(element), FALSE);
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.h b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
index d3c62dd867..cb9531e2dd 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.h
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
@@ -202,6 +202,8 @@ struct _GSTLALTransferFunction {
gboolean use_fir_fft;
double df;
double f0;
+ guint a;
+ guint b;
};
diff --git a/gstlal-calibration/lib/gstlal_firtools.c b/gstlal-calibration/lib/gstlal_firtools.c
index b610eb7cf6..0195aa8911 100644
--- a/gstlal-calibration/lib/gstlal_firtools.c
+++ b/gstlal-calibration/lib/gstlal_firtools.c
@@ -1139,6 +1139,8 @@ GSTLAL_IRFFT(long, double, l, );
LONG complex DTYPE *gstlal_rfft_comb_ ## LONG ## DTYPE(LONG DTYPE *td_data, guint N, guint a, guint b, LONG DTYPE f0, gboolean free_input, guint *N_out) { \
\
guint factor, N_blocks, *prime_factors, num_factors, i, j, M, M_min, M_num_factors, *M_prime_factors, M_out, M_in, rfft_N_out; \
+ if(N_out == NULL) \
+ N_out = g_malloc(sizeof(guint)); \
gboolean M_prime_factors_need_freed, exp_array_need_freed; \
long complex double *rot_array, *rot_td_data, *M_exp_array, *M_exp_array2, *exp_array, *long_fd_data, *mini_fd_data, *fd_data; \
LONG DTYPE rot_pow; \
@@ -1192,9 +1194,9 @@ LONG complex DTYPE *gstlal_rfft_comb_ ## LONG ## DTYPE(LONG DTYPE *td_data, guin
/* Find prime factors */ \
prime_factors = find_prime_factors(b, &num_factors); \
\
+ rfft_N_out = b / 2 + 1; \
+ *N_out = (rfft_N_out + a - 1) / a; \
if(f0 > 0) { \
- rfft_N_out = b / 2; \
- *N_out = (guint) (rfft_N_out + a - 1 - f0 * a) / a; \
/* Apply phase rotation to shift the frequencies being measured. */ \
rot_array = find_exp_array(N, FALSE); \
rot_td_data = g_malloc(N * sizeof(long complex double)); \
@@ -1267,9 +1269,6 @@ LONG complex DTYPE *gstlal_rfft_comb_ ## LONG ## DTYPE(LONG DTYPE *td_data, guin
} \
\
} else { \
- rfft_N_out = b / 2 + 1; \
- *N_out = (rfft_N_out + a - 1) / a; \
- \
/* Check if we will need to use gstlal_prime_rfft_() for this */ \
if(prime_factors[num_factors - 2] >= 607) { \
/* Find the first member greater than or equal to 607 */ \
diff --git a/gstlal-calibration/tests/check_calibration/Makefile b/gstlal-calibration/tests/check_calibration/Makefile
index 3f04a6026c..bccfc16add 100644
--- a/gstlal-calibration/tests/check_calibration/Makefile
+++ b/gstlal-calibration/tests/check_calibration/Makefile
@@ -8,12 +8,13 @@ IFO = H
# determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
OBSRUN = O3
-START = $(shell echo 1269090018 | bc)
+START = $(shell echo 1269114368 | bc)
# 1237831461 start of O3A
# 1256655618 start of O3B
# 1238288418
# 1269090018
-END = $(shell echo 1269091018 | bc)
+# 1269114645-21449
+END = $(shell echo 1269136384 | bc)
# 1253977218 end of O3A
# 1269367218 end of O3B
# 1238338818
@@ -48,7 +49,9 @@ DCSTESTCONFIGS = Filters/O3/GDSFilters/H1DCS_C01_1252083618_AllCorrectionsTest.i
GDSSHMCONFIGS = Filters/O3/GDSFilters/H1GDS_1258216456_testing.ini
GDSOLDCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_old.ini
-all: $(IFO)1_hoft_GDS_SHM_frames
+all: TDCFs_pcal2darm_plots
+# TDCFs_pcal2darm_plots
+# pcal_DCS_transfer_functions
#approxKappasErrorPlot
#$(IFO)1_hoft_DCS_STATSUB_frames.cache $(IFO)1_hoft_DCS_NONSTATSUB_frames.cache
@@ -206,7 +209,7 @@ approxKappasErrorPlot: $(IFO)1_C01_frames.cache
#1239472998,1251057623,1256655618,1278428752
TDCFs_pcal2darm_plots: $(IFO)1_easy_raw_frames.cache $(IFO)1_SCALAR_CORR_frames.cache $(IFO)1_FCC_CORR_frames.cache $(IFO)1_FCCFS_CORR_frames.cache $(IFO)1_ALL_CORR_frames.cache
- python3 pcal2darm_timeseries.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --raw-frame-cache $(IFO)1_easy_raw_frames.cache --gstlal-frame-cache-list '$(IFO)1_SCALAR_CORR_frames.cache,$(IFO)1_FCC_CORR_frames.cache,$(IFO)1_FCCFS_CORR_frames.cache,$(IFO)1_ALL_CORR_frames.cache' --config-file '$(ALL_CORR_CONFIGS)' --pcal-channel-name CAL-PCALY_RX_PD_OUT_DQ --gstlal-channel-list 'DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN' --labels '{\rm Scalars},+f_{\rm cc},+f_{\rm s}+Q,+\tau_i' --latex-labels --pcal-time-advance 0.00006103515625 --show-stats
+ python3 pcal2darm_timeseries.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --raw-frame-cache $(IFO)1_easy_raw_frames.cache --gstlal-frame-cache-list '$(IFO)1_SCALAR_CORR_frames.cache,$(IFO)1_FCC_CORR_frames.cache,$(IFO)1_FCCFS_CORR_frames.cache,$(IFO)1_ALL_CORR_frames.cache' --config-file-list '$(SCALAR_CORR_CONFIGS),$(FCC_CORR_CONFIGS),$(FCCFS_CORR_CONFIGS),$(ALL_CORR_CONFIGS)' --pcal-channel-name CAL-PCALY_RX_PD_OUT_DQ --gstlal-channel-list 'DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN' --labels '{\rm Multipliers},+f_{\rm cc},+f_{\rm s}+Q,+\tau_i' --latex-labels --pcal-time-advance 0.00006103515625 --show-stats
detect_actuation_timing_change: $(IFO)1_C01_frames.cache
python3 detect_kappa_change.py --gps-start-time $(START) --gps-end-time $(END) --ifo $(IFO)1 --frame-cache $(IFO)1_C01_frames.cache --statevector-channel DCS-CALIB_STATE_VECTOR_C01 filename $(IFO)_actuation_timing_changes_$(START)-$(END).txt
@@ -219,7 +222,7 @@ pcal_GDS_transfer_functions: $(IFO)1_easy_raw_frames.cache $(IFO)1_hoft_GDS_fram
python3 plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_RX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_hoft_GDS_frames.cache --numerator-channel-list GDS-CALIB_STRAIN --config-file $(GDSCONFIGS) --use-median --labels GDS_NO_TDCFs
pcal_DCS_transfer_functions: $(IFO)1_SCALAR_CORR_frames.cache $(IFO)1_FCC_CORR_frames.cache $(IFO)1_FCCFS_CORR_frames.cache $(IFO)1_ALL_CORR_frames.cache
- python3 plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_RX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_SCALAR_CORR_frames.cache,$(IFO)1_FCC_CORR_frames.cache,$(IFO)1_FCCFS_CORR_frames.cache,$(IFO)1_ALL_CORR_frames.cache --numerator-channel-list DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN --config-file $(ALL_CORR_CONFIGS) --use-median --labels '{\rm Scalars},+f_{\rm cc},+f_{\rm s}+Q,+\tau_i' --latex-labels
+ python3 plot_transfer_function.py --gps-start-time 1269021699 --gps-end-time 1269021877 --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_RX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_SCALAR_CORR_frames.cache,$(IFO)1_FCC_CORR_frames.cache,$(IFO)1_FCCFS_CORR_frames.cache,$(IFO)1_ALL_CORR_frames.cache --numerator-channel-list DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN,DCS-CALIB_STRAIN --config-file $(ALL_CORR_CONFIGS) --use-median --labels '{\rm Multipliers},+f_{\rm cc},+f_{\rm s}+Q,+\tau_i' --latex-labels
pcal_DCS_EXACTKAPPAS_transfer_functions: $(IFO)1_easy_raw_frames.cache $(IFO)1_hoft_DCS_APPROXKAPPAS_frames.cache $(IFO)1_hoft_DCS_EXACTKAPPAS_frames.cache $(IFO)1_C01_frames.cache
python3 plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_RX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_C01_frames.cache,$(IFO)1_hoft_DCS_APPROXKAPPAS_frames.cache,$(IFO)1_hoft_DCS_EXACTKAPPAS_frames.cache --numerator-channel-list DCS-CALIB_STRAIN_C01,DCS-CALIB_STRAINAPPROXKAPPAS,DCS-CALIB_STRAINEXACTKAPPAS --config-file $(DCSEXACTKAPPASCONFIGS) --labels 'C01,Approx,Exact' --use-median
diff --git a/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py b/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
index d33790cfc2..2c84195df2 100644
--- a/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
+++ b/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
@@ -331,7 +331,7 @@ else:
plot_labels.append("{\\rm %s}" % label.replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_').replace(' ', '\ '))
# Read data from files and plot it
-colors = ['red', 'limegreen', 'mediumblue', 'gold', 'b', 'm'] # Hopefully the user will not want to plot more than six datasets on one plot.
+colors = ['tomato', 'green', 'mediumblue', 'gold', 'b', 'm'] # Hopefully the user will not want to plot more than six datasets on one plot.
channels = calcs_channels
channels.extend(gstlal_channels)
for i in range(0, len(frequencies)):
@@ -374,14 +374,14 @@ for i in range(0, len(frequencies)):
plt.ylabel(r'${\rm Magnitude}$')
magnitude_range = options.magnitude_ranges.split(';')[i]
ticks_and_grid(plt.gca(), ymin = float(magnitude_range.split(',')[0]), ymax = float(magnitude_range.split(',')[1]))
- leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'upper right')
+ leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'lower right')
leg.get_frame().set_alpha(0.8)
plt.subplot(2, len(frequencies), len(frequencies) + i + 1)
if options.show_stats:
plt.plot(times[0], phases[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'$%s \ [\mu_{1/2} = %0.3f^{\circ}, \sigma = %0.3f^{\circ}]$' % (plot_labels[0], numpy.median(phases[0]), stdmed(phases[0])))
else:
plt.plot(times[0], phases[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'$%s$' % (plot_labels[0]))
- leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'upper right')
+ leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'lower right')
leg.get_frame().set_alpha(0.8)
if i == 0:
plt.ylabel(r'${\rm Phase \ [deg]}$')
@@ -404,14 +404,14 @@ for i in range(0, len(frequencies)):
plt.plot(times[j], magnitudes[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'$%s \ [\mu_{1/2} = %0.4f, \sigma = %0.4f]$' % (plot_labels[j], numpy.median(magnitudes[j]), stdmed(magnitudes[j])))
else:
plt.plot(times[j], magnitudes[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'$%s$' % (plot_labels[j]))
- leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'upper right')
+ leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'lower right')
leg.get_frame().set_alpha(0.8)
plt.subplot(2, len(frequencies), len(frequencies) + i + 1)
if options.show_stats:
plt.plot(times[j], phases[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'$%s \ [\mu_{1/2} = %0.3f^{\circ}, \sigma = %0.3f^{\circ}]$' % (plot_labels[j], numpy.median(phases[j]), stdmed(phases[j])))
else:
plt.plot(times[j], phases[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'$%s$' % (plot_labels[j]))
- leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'upper right')
+ leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'lower right')
leg.get_frame().set_alpha(0.8)
plt.savefig("%s_deltal_over_pcal%s_%d-%d.png" % (ifo, options.file_name_suffix, int(t_start), int(dur_in_seconds)))
plt.savefig("%s_deltal_over_pcal%s_%d-%d.pdf" % (ifo, options.file_name_suffix, int(t_start), int(dur_in_seconds)))
diff --git a/gstlal-calibration/tests/check_calibration/plot_transfer_function.py b/gstlal-calibration/tests/check_calibration/plot_transfer_function.py
index 7ebf8be10a..c05e1fc0d5 100644
--- a/gstlal-calibration/tests/check_calibration/plot_transfer_function.py
+++ b/gstlal-calibration/tests/check_calibration/plot_transfer_function.py
@@ -340,7 +340,7 @@ for i in range(len(real_poles)):
colors = []
# Start with defaults to use if we don't recognize what we are plotting.
-default_colors = ['tomato', 'red', 'limegreen', 'gold', 'orange', 'aqua', 'magenta', 'blue']
+default_colors = ['tomato', 'green', 'blue', 'gold', 'orange', 'aqua', 'magenta', 'blue']
# Use specific colors for known version of calibration
C02_labels = ['C02', 'c02']
--
GitLab