diff --git a/gstlal-calibration/configure.ac b/gstlal-calibration/configure.ac index b0185a06b88e4b02df24243daf4ac78cb473f1f9..241c763af3c3f299e9f3671de671fbd1bfa9901b 100644 --- a/gstlal-calibration/configure.ac +++ b/gstlal-calibration/configure.ac @@ -20,6 +20,7 @@ AC_CONFIG_FILES([ \ debian/Makefile \ gst/Makefile \ gst/lal/Makefile \ + gst/cmath/Makefile \ python/Makefile \ tests/Makefile \ ]) diff --git a/gstlal-calibration/gst/Makefile.am b/gstlal-calibration/gst/Makefile.am index 13bae38f0d2ae680a907a2d1bc5099c2f1667c84..740c3abf69ed2473a6eb1ad18c4689cca5f6de91 100644 --- a/gstlal-calibration/gst/Makefile.am +++ b/gstlal-calibration/gst/Makefile.am @@ -1 +1 @@ -SUBDIRS = lal +SUBDIRS = lal cmath diff --git a/gstlal-calibration/gst/cmath/Makefile.am b/gstlal-calibration/gst/cmath/Makefile.am new file mode 100644 index 0000000000000000000000000000000000000000..a2b6c7acb6104621340319de8e2de28ca60205f5 --- /dev/null +++ b/gstlal-calibration/gst/cmath/Makefile.am @@ -0,0 +1,18 @@ +plugin_LTLIBRARIES = libgstcmath.la + +# sources used to compile this plug-in +libgstcmath_la_SOURCES = \ + cmath.c \ + cmath_base.h cmath_base.c \ + cmath_cabs.c \ + cmath_cexp.c \ + cmath_cln.c \ + cmath_clog.c \ + cmath_clog10.c \ + cmath_cpow.c +libgstcmath_la_CPPFLAGS = $(AM_CPPFLAGS) $(PYTHON_CPPFLAGS) +libgstcmath_la_CFLAGS = $(AM_CFLAGS) $(LAL_CFLAGS) $(GSTLAL_CFLAGS) $(gstreamer_CFLAGS) $(gstreamer_audio_CFLAGS) +libgstcmath_la_LDFLAGS = $(AM_LDFLAGS) $(LAL_LIBS) $(GSTLAL_LIBS) $(PYTHON_LIBS) $(gstreamer_LIBS) $(gstreamer_audio_LIBS) $(GSTLAL_PLUGIN_LDFLAGS) + +# headers we need but don't want installed +noinst_HEADERS = cmath_base.h diff --git a/gstlal-calibration/gst/cmath/cmath.c b/gstlal-calibration/gst/cmath/cmath.c new file mode 100644 index 0000000000000000000000000000000000000000..6bbf425bdec5fb056393429b12546ccc1b42ed95 --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath.c @@ -0,0 +1,105 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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 + * + * ============================================================================ + */ + + +/* + * Stuff from GStreamer + */ + + +#include <gst/gst.h> + + +/* + * Our own stuff + */ + + +#include <cmath_base.h> + + +/* + * ============================================================================ + * + * Plugin Entry Point + * + * ============================================================================ + */ + + +GType cmath_cabs_get_type (void); +GType cmath_cexp_get_type (void); +GType cmath_cln_get_type (void); +GType cmath_clog_get_type (void); +GType cmath_clog10_get_type (void); +GType cmath_cpow_get_type (void); + + +static gboolean +plugin_init (GstPlugin *plugin) +{ + struct + { + const gchar *name; + GType type; + } *element, elements[] = { + { + "cmath_base", CMATH_BASE_TYPE}, { + "cabs", cmath_cabs_get_type ()}, { + "cexp", cmath_cexp_get_type ()}, { + "cln", cmath_cln_get_type ()}, { + "clog", cmath_clog_get_type ()}, { + "clog10", cmath_clog10_get_type ()}, { + "cpow", cmath_cpow_get_type ()}, { + NULL, 0},}; + + /* + * Tell GStreamer about the elements. + */ + + for (element = elements; element->name; element++) + if (!gst_element_register (plugin, element->name, GST_RANK_NONE, + element->type)) + return FALSE; + + /* + * Done. + */ + + return TRUE; +} + + +/* + * This is the structure that gst-register looks for. + */ + + +GST_PLUGIN_DEFINE (GST_VERSION_MAJOR, GST_VERSION_MINOR, cmath, + "Complex arithmetic elements", plugin_init, PACKAGE_VERSION, "GPL", + PACKAGE_NAME, "http://www.lsc-group.phys.uwm.edu/daswg") diff --git a/gstlal-calibration/gst/cmath/cmath_base.c b/gstlal-calibration/gst/cmath/cmath_base.c new file mode 100644 index 0000000000000000000000000000000000000000..e943569923aeb2165a52ab9b606cb5b3ed2516c9 --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_base.c @@ -0,0 +1,181 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + +/* + * Start here for explanation of class: + * https://gstreamer.freedesktop.org/data/doc/gstreamer/head/gstreamer-libs/html/GstBaseTransform.html#GstBaseTransform-struct + */ + +#include <cmath_base.h> + + +/* + * ============================================================================ + * + * Type Support + * + * ============================================================================ + */ + + +#define CAPS \ + "audio/x-raw, " \ + "format = (string) { F32LE, F32BE, F64LE, F64BE, Z64LE, Z64BE, Z128LE, Z128BE }, " \ + "rate = (int) [1, MAX], " \ + "channels = (int) [1, MAX], " \ + "layout = (string) {interleaved, non-interleaved}, " \ + "channel-mask = (bitmask) 0" + + +static void +base_init(gpointer klass) +{ + GstElementClass* gstelement_class = GST_ELEMENT_CLASS(klass); + GstBaseTransformClass* gstbasetransform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gst_element_class_set_details_simple(gstelement_class, + "CMath operation base class", + "Filter/Audio", + "Base class for elements that provide complex arithmetic operations", + "Aaron Viets <aaron.viets@ligo.org>, Leo Singer <leo.singer@ligo.org>"); + + gstbasetransform_class -> set_caps = set_caps; + + static GstStaticPadTemplate sink_factory = GST_STATIC_PAD_TEMPLATE("sink", + GST_PAD_SINK, + GST_PAD_ALWAYS, + GST_STATIC_CAPS(CAPS) + ); + + static GstStaticPadTemplate src_factory = GST_STATIC_PAD_TEMPLATE("src", + GST_PAD_SRC, + GST_PAD_ALWAYS, + GST_STATIC_CAPS(CAPS) + ); + + gstbasetransform_class -> set_caps = set_caps; + + gst_element_class_add_pad_template(gstelement_class, gst_static_pad_template_get(&src_factory)); + gst_element_class_add_pad_template(gstelement_class, gst_static_pad_template_get(&sink_factory)); +} + +GType +cmath_base_get_type(void) +{ + static GType type = 0; + + if(!type) { + static const GTypeInfo info = { + .class_size = sizeof(CMathBaseClass), + .base_init = base_init, + .instance_size = sizeof(CMathBase), + }; + type = + g_type_register_static(GST_TYPE_BASE_TRANSFORM, "CMathBase", &info, 0); + } + + return type; +} + + +/* + * ============================================================================ + * + * Caps Negotiation + * + * ============================================================================ + */ + +/* + * set_caps() is called when the caps are re-negotiated. Can return + * false if we do not like what we see. + */ + +gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, + GstCaps *outcaps) +{ + gboolean success = TRUE; + static char *formats[] = {"Z128LE", "Z64LE", "F64LE", "F32LE"}; + static int is_complex[] = {1, 1, 0, 0}; + static int bits[] = {128, 64, 64, 32}; + const gchar *format; + const gchar *interleave; + gint channels; + gint rate; + + CMathBase *element = CMATH_BASE(trans); + + /* Check the incaps to see if it contains e.g. Z128, etc. */ + GstStructure *str = gst_caps_get_structure(incaps, 0); + g_assert(str); + + if(gst_structure_has_field(str, "format")) { + format = gst_structure_get_string(str, "format"); + } else { + GST_ERROR_OBJECT(trans, "No format! Cannot set element caps.\n"); + return 0; + } + + guint i; + for(i = 0; i < sizeof(formats) / sizeof(*formats); i++) { + if(!strcmp(format, formats[i])) { + element->is_complex = is_complex[i]; + element->bits = bits[i]; + goto found_format; + } + } + + GST_ERROR_OBJECT(trans, "Could not find format [%s].\n", format); + element->is_complex = -1; + element->bits = -1; + element->channels = -1; + element->interleave = NULL; + return 0; + + found_format: + + if(gst_structure_has_field(str, "layout")) { + interleave = gst_structure_get_string(str, "layout"); + } else { + GST_ERROR_OBJECT(trans, "No layout! Cannot set element caps.\n"); + return 0; + } + + element->interleave = interleave; + + if(gst_structure_has_field(str, "channels")) { + gst_structure_get_int(str, "channels", &channels); + } else { + GST_ERROR_OBJECT(trans, "No channels! Cannot set element caps.\n"); + return 0; + } + + element->channels = channels; + + if(gst_structure_has_field(str, "rate")) { + gst_structure_get_int(str, "rate", &rate); + } else { + GST_ERROR_OBJECT(trans, "No rate! Cannot set element caps.\n"); + return 0; + } + + element->rate = rate; + + return success; +} diff --git a/gstlal-calibration/gst/cmath/cmath_base.h b/gstlal-calibration/gst/cmath/cmath_base.h new file mode 100644 index 0000000000000000000000000000000000000000..7c2c42a2098a8ffa05d2f820a928e4295586c7cd --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_base.h @@ -0,0 +1,77 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + + +#ifndef __CMATH_BASE_H__ +#define __CMATH_BASE_H__ + + +#include <glib.h> +#include <gst/gst.h> +#include <gst/base/gstbasetransform.h> +#include <stdio.h> +#include <string.h> +#include <complex.h> +#include <stdlib.h> +#include <math.h> + +G_BEGIN_DECLS +#define CMATH_BASE_TYPE \ + (cmath_base_get_type()) +#define CMATH_BASE(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj), CMATH_BASE_TYPE, CMathBase)) +#define CMATH_BASE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass), CMATH_BASE_TYPE, CMathBaseClass)) +#define GST_IS_CMATH_BASE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj), CMATH_BASE_TYPE)) +#define GST_IS_CMATH_BASE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass), CMATH_BASE_TYPE)) + +typedef struct _CMathBase CMathBase; +typedef struct _CMathBaseClass CMathBaseClass; + +struct _CMathBase +{ + GstBaseTransform element; + + int is_complex; + int bits; + int channels; + int rate; + const gchar *interleave; +}; + +struct _CMathBaseClass +{ + GstBaseTransformClass parent_class; +}; + +GType cmath_base_get_type(void); + +/* + * set_caps() is called when the caps are re-negotiated. Can return + * false if we do not like what we see. + */ + +gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, + GstCaps *outcaps); + + +G_END_DECLS +#endif /* __CMATH_BASE_H__ */ diff --git a/gstlal-calibration/gst/cmath/cmath_cabs.c b/gstlal-calibration/gst/cmath/cmath_cabs.c new file mode 100644 index 0000000000000000000000000000000000000000..689a4ea18029c37838bd037db0702baca030ae4d --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_cabs.c @@ -0,0 +1,168 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + + +#include <cmath_base.h> + +#define TYPE_CMATH_CABS \ + (cmath_cabs_get_type()) +#define CMATH_CABS(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj),TYPE_CMATH_CABS,CMathCAbs)) +#define CMATH_CABS_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass),TYPE_CMATH_CABS,CMathCAbsClass)) +#define IS_PLUGIN_TEMPLATE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj),TYPE_CMATH_CABS)) +#define IS_PLUGIN_TEMPLATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass),TYPE_CMATH_CABS)) + +typedef struct _CMathCAbs CMathCAbs; +typedef struct _CMathCAbsClass CMathCAbsClass; + +GType +cmath_cabs_get_type(void); + +struct _CMathCAbs +{ + CMathBase cmath_base; +}; + +struct _CMathCAbsClass +{ + CMathBaseClass parent_class; +}; + + +/* + * ============================================================================ + * + * GstBaseTransform vmethod Implementations + * + * ============================================================================ + */ + +/* An in-place transform really does the same thing as the chain function */ + +static GstFlowReturn +transform_ip(GstBaseTransform *trans, GstBuffer *buf) +{ + CMathCAbs* element = CMATH_CABS(trans); + int bits = element -> cmath_base.bits; + int is_complex = element -> cmath_base.is_complex; + + /* + * Debugging + * + * GstObject* element_gstobj = GST_OBJECT(trans); + * int channels = element -> cmath_base.channels; + * int rate = element -> cmath_base.rate; + * g_print("[%s]: passing GstBuffer: ", element_gstobj->name); + * g_print("%d channels, ", channels); + * g_print("%d bits, ", bits); + * g_print("rate: %d, ", rate); + */ + + GstMapInfo info; + if(!gst_buffer_map(buf, &info, GST_MAP_READWRITE)) { + GST_ERROR_OBJECT(trans, "gst_buffer_map failed\n"); + } + gpointer data = info.data; + gpointer data_end = data + info.size; + + if(is_complex == 1) { + + if(bits == 128) { + /* g_print("COMPLEX FLOAT128\n"); */ + double complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = cabs(*ptr); + } + } else if(bits == 64) { + /* g_print("COMPLEX FLOAT64\n"); */ + float complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = cabsf(*ptr); + } + } else { + g_assert_not_reached(); + } + } else if(is_complex == 0) { + + if(bits == 64) { + /* g_print("REAL FLOAT64\n"); */ + double *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = fabs(*ptr); + } + } else if(bits == 32) { + /* g_print("REAL FLOAT32\n"); */ + float *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = fabsf(*ptr); + } + } else { + g_assert_not_reached(); + } + } else { + g_assert_not_reached(); + } + gst_buffer_unmap(buf, &info); + return GST_FLOW_OK; +} + + +/* + * ============================================================================ + * + * Type Support + * + * ============================================================================ + */ + +/* Initialize the plugin's class */ +static void +cmath_cabs_class_init(gpointer klass, gpointer klass_data) +{ + GstBaseTransformClass *basetransform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gst_element_class_set_details_simple(GST_ELEMENT_CLASS(klass), + "Absolute value", + "Filter/Audio", + "Calculate absolute value, y = |x|", + "Aaron Viets <aaron.viets@ligo.org>, Leo Singer <leo.singer@ligo.org>"); + + basetransform_class -> transform_ip = GST_DEBUG_FUNCPTR(transform_ip); + basetransform_class -> set_caps = GST_DEBUG_FUNCPTR(set_caps); +} + +GType +cmath_cabs_get_type(void) +{ + static GType type = 0; + + if(!type) { + static const GTypeInfo info = { + .class_size = sizeof(CMathBaseClass), + .class_init = cmath_cabs_class_init, + .instance_size = sizeof(CMathBase), + }; + type = g_type_register_static(CMATH_BASE_TYPE, "CMathCAbs", &info, 0); + } + + return type; +} diff --git a/gstlal-calibration/gst/cmath/cmath_cexp.c b/gstlal-calibration/gst/cmath/cmath_cexp.c new file mode 100644 index 0000000000000000000000000000000000000000..6031b7effa882c3939185ee273770904d955b97c --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_cexp.c @@ -0,0 +1,164 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + + +#include <cmath_base.h> + +#define TYPE_CMATH_CEXP \ + (cmath_cexp_get_type()) +#define CMATH_CEXP(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj),TYPE_CMATH_CEXP,CMathCExp)) +#define CMATH_CEXP_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass),TYPE_CMATH_CEXP,CMathCExpClass)) +#define IS_PLUGIN_TEMPLATE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj),TYPE_CMATH_CEXP)) +#define IS_PLUGIN_TEMPLATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass),TYPE_CMATH_CEXP)) + +typedef struct _CMathCExp CMathCExp; +typedef struct _CMathCExpClass CMathCExpClass; + +GType +cmath_cexp_get_type(void); + +struct _CMathCExp +{ + CMathBase cmath_base; +}; + +struct _CMathCExpClass +{ + CMathBaseClass parent_class; +}; + + +/* + * ============================================================================ + * + * GstBaseTransform vmethod Implementations + * + * ============================================================================ + */ + +/* An in-place transform really does the same thing as the chain function */ + +static GstFlowReturn +transform_ip(GstBaseTransform *trans, GstBuffer *buf) +{ + CMathCExp* element = CMATH_CEXP(trans); + int bits = element -> cmath_base.bits; + int is_complex = element -> cmath_base.is_complex; + + /* + * Debugging + * + * GstObject* element_gstobj = GST_OBJECT(trans); + * int channels = element -> cmath_base.channels; + * int rate = element -> cmath_base.rate; + * g_print("[%s]: passing GstBuffer: ", element_gstobj->name); + * g_print("%d channels, ", channels); + * g_print("%d bits, ", bits); + * g_print("rate: %d, ", rate); + */ + + GstMapInfo info; + if(!gst_buffer_map(buf, &info, GST_MAP_READWRITE)) { + GST_ERROR_OBJECT(trans, "gst_buffer_map failed\n"); + } + gpointer data = info.data; + gpointer data_end = data + info.size; + + if(is_complex == 1) { + + if(bits == 128) { + double complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = cexp(*ptr); + } + } else if(bits == 64) { + float complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = cexpf(*ptr); + } + } else { + g_assert_not_reached(); + } + } else if(is_complex == 0) { + + if(bits == 64) { + double *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = exp(*ptr); + } + } else if(bits == 32) { + float *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = expf(*ptr); + } + } else { + g_assert_not_reached(); + } + } else { + g_assert_not_reached(); + } + gst_buffer_unmap(buf, &info); + return GST_FLOW_OK; +} + + +/* + * ============================================================================ + * + * Type Support + * + * ============================================================================ + */ + +/* Initialize the plugin's class */ +static void +cmath_cexp_class_init(gpointer klass, gpointer klass_data) +{ + GstBaseTransformClass *basetransform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gst_element_class_set_details_simple(GST_ELEMENT_CLASS(klass), + "Natural exponent", + "Filter/Audio", + "Calculate natural exponent, y = e^x", + "Aaron Viets <aaron.viets@ligo.org>, Leo Singer <leo.singer@ligo.org>"); + + basetransform_class -> transform_ip = GST_DEBUG_FUNCPTR(transform_ip); + basetransform_class -> set_caps = GST_DEBUG_FUNCPTR(set_caps); +} + +GType +cmath_cexp_get_type(void) +{ + static GType type = 0; + + if(!type) { + static const GTypeInfo info = { + .class_size = sizeof(CMathBaseClass), + .class_init = cmath_cexp_class_init, + .instance_size = sizeof(CMathBase), + }; + type = g_type_register_static(CMATH_BASE_TYPE, "CMathCExp", &info, 0); + } + + return type; +} diff --git a/gstlal-calibration/gst/cmath/cmath_cln.c b/gstlal-calibration/gst/cmath/cmath_cln.c new file mode 100644 index 0000000000000000000000000000000000000000..256ca443d8e500100f912b05eb061a11c829734a --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_cln.c @@ -0,0 +1,164 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + + +#include <cmath_base.h> + +#define TYPE_CMATH_CLN \ + (cmath_cln_get_type()) +#define CMATH_CLN(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj),TYPE_CMATH_CLN,CMathCLn)) +#define CMATH_CLN_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass),TYPE_CMATH_CLN,CMathCLnClass)) +#define IS_PLUGIN_TEMPLATE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj),TYPE_CMATH_CLN)) +#define IS_PLUGIN_TEMPLATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass),TYPE_CMATH_CLN)) + +typedef struct _CMathCLn CMathCLn; +typedef struct _CMathCLnClass CMathCLnClass; + +GType +cmath_cln_get_type(void); + +struct _CMathCLn +{ + CMathBase cmath_base; +}; + +struct _CMathCLnClass +{ + CMathBaseClass parent_class; +}; + + +/* + * ============================================================================ + * + * GstBaseTransform vmethod Implementations + * + * ============================================================================ + */ + +/* An in-place transform really does the same thing as the chain function */ + +static GstFlowReturn +transform_ip(GstBaseTransform *trans, GstBuffer *buf) +{ + CMathCLn* element = CMATH_CLN(trans); + int bits = element -> cmath_base.bits; + int is_complex = element -> cmath_base.is_complex; + + /* + * Debugging + * + * GstObject* element_gstobj = GST_OBJECT(trans); + * int channels = element -> cmath_base.channels; + * int rate = element -> cmath_base.rate; + * g_print("[%s]: passing GstBuffer: ", element_gstobj->name); + * g_print("%d channels, ", channels); + * g_print("%d bits, ", bits); + * g_print("rate: %d, ", rate); + */ + + GstMapInfo info; + if(!gst_buffer_map(buf, &info, GST_MAP_READWRITE)) { + GST_ERROR_OBJECT(trans, "gst_buffer_map failed\n"); + } + gpointer data = info.data; + gpointer data_end = data + info.size; + + if(is_complex == 1) { + + if(bits == 128) { + double complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = clog(*ptr); + } + } else if(bits == 64) { + float complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = clogf(*ptr); + } + } else { + g_assert_not_reached(); + } + } else if(is_complex == 0) { + + if(bits == 64) { + double *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = log(*ptr); + } + } else if(bits == 32) { + float *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = logf(*ptr); + } + } else { + g_assert_not_reached(); + } + } else { + g_assert_not_reached(); + } + gst_buffer_unmap(buf, &info); + return GST_FLOW_OK; +} + + +/* + * ============================================================================ + * + * Type Support + * + * ============================================================================ + */ + +/* Initialize the plugin's class */ +static void +cmath_cln_class_init(gpointer klass, gpointer klass_data) +{ + GstBaseTransformClass *basetransform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gst_element_class_set_details_simple(GST_ELEMENT_CLASS(klass), + "Natural logarithm", + "Filter/Audio", + "Calculate natural logarithm, y = ln x", + "Aaron Viets <aaron.viets@ligo.org>, Leo Singer <leo.singer@ligo.org>"); + + basetransform_class -> transform_ip = GST_DEBUG_FUNCPTR(transform_ip); + basetransform_class -> set_caps = GST_DEBUG_FUNCPTR(set_caps); +} + +GType +cmath_cln_get_type(void) +{ + static GType type = 0; + + if(!type) { + static const GTypeInfo info = { + .class_size = sizeof(CMathBaseClass), + .class_init = cmath_cln_class_init, + .instance_size = sizeof(CMathBase), + }; + type = g_type_register_static(CMATH_BASE_TYPE, "CMathCLn", &info, 0); + } + + return type; +} diff --git a/gstlal-calibration/gst/cmath/cmath_clog.c b/gstlal-calibration/gst/cmath/cmath_clog.c new file mode 100644 index 0000000000000000000000000000000000000000..29e329717102a73b93a38d801f1ff6ae907e77b1 --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_clog.c @@ -0,0 +1,221 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + + +#include <cmath_base.h> + +#define TYPE_CMATH_CLOG \ + (cmath_clog_get_type()) +#define CMATH_CLOG(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj),TYPE_CMATH_CLOG,CMathCLog)) +#define CMATH_CLOG_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass),TYPE_CMATH_CLOG,CMathCLogClass)) +#define IS_PLUGIN_TEMPLATE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj),TYPE_CMATH_CLOG)) +#define IS_PLUGIN_TEMPLATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass),TYPE_CMATH_CLOG)) + +typedef struct _CMathCLog CMathCLog; +typedef struct _CMathCLogClass CMathCLogClass; + +GType +cmath_clog_get_type(void); + +struct _CMathCLog +{ + CMathBase cmath_base; + double base; +}; + +struct _CMathCLogClass +{ + CMathBaseClass parent_class; +}; + + +/* + * ============================================================================ + * + * GstBaseTransform vmethod Implementations + * + * ============================================================================ + */ + +/* An in-place transform really does the same thing as the chain function */ + +static GstFlowReturn +transform_ip(GstBaseTransform *trans, GstBuffer *buf) +{ + CMathCLog* element = CMATH_CLOG(trans); + int bits = element -> cmath_base.bits; + int is_complex = element -> cmath_base.is_complex; + + /* Debugging + * + * GstObject* element_gstobj = GST_OBJECT(trans); + * int channels = element -> cmath_base.channels; + * int rate = element -> cmath_base.rate; + * g_print("[%s]: passing GstBuffer: ", element_gstobj->name); + * g_print("%d channels, ", channels); + * g_print("%d bits, ", bits); + * g_print("rate: %d, ", rate); + */ + + GstMapInfo info; + if(!gst_buffer_map(buf, &info, GST_MAP_READWRITE)) { + GST_ERROR_OBJECT(trans, "gst_buffer_map failed\n"); + } + gpointer data = info.data; + gpointer data_end = data + info.size; + + const double n = element -> base; + const double m = 1. / log(n); + + if(is_complex == 1) { + + if(bits == 128) { + double complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = clog(*ptr) * m; + } + } else if(bits == 64) { + float complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = clogf(*ptr) * m; + } + } else { + g_assert_not_reached(); + } + } else if(is_complex == 0) { + + if(bits == 64) { + double *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = log(*ptr) * m; + } + } else if(bits == 32) { + float *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = logf(*ptr) * m; + } + } else { + g_assert_not_reached(); + } + } else { + g_assert_not_reached(); + } + gst_buffer_unmap(buf, &info); + return GST_FLOW_OK; +} + + +/* + * ============================================================================ + * + * Type Support + * + * ============================================================================ + */ + +/* Set the exponent */ +enum property +{ + PROP_BASE = 1, +}; + + +static void +set_property(GObject * object, enum property id, const GValue * value, + GParamSpec * pspec) +{ + CMathCLog *element = CMATH_CLOG(object); + + GST_OBJECT_LOCK(element); + + switch(id) { + case PROP_BASE: + element->base = g_value_get_double(value); + break; + } + + GST_OBJECT_UNLOCK(element); +} + + +static void +get_property(GObject * object, enum property id, GValue * value, + GParamSpec * pspec) +{ + CMathCLog *element = CMATH_CLOG(object); + + GST_OBJECT_LOCK(element); + + switch(id) { + case PROP_BASE: + g_value_set_double(value, element->base); + break; + } + + GST_OBJECT_UNLOCK(element); +} + +/* Initialize the plugin's class */ +static void +cmath_clog_class_init(gpointer klass, gpointer klass_data) +{ + GObjectClass *gobject_class = G_OBJECT_CLASS(klass); + GstBaseTransformClass *basetransform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gst_element_class_set_details_simple(GST_ELEMENT_CLASS(klass), + "Logarithm with base k", + "Filter/Audio", + "Calculate logarithm, y = log_k x", + "Aaron Viets <aaron.viets@ligo.org>, Leo Singer <leo.singer@ligo.org>"); + + basetransform_class -> transform_ip = GST_DEBUG_FUNCPTR(transform_ip); + basetransform_class -> set_caps = GST_DEBUG_FUNCPTR(set_caps); + + gobject_class->get_property = GST_DEBUG_FUNCPTR(get_property); + gobject_class->set_property = GST_DEBUG_FUNCPTR(set_property); + + g_object_class_install_property(gobject_class, + PROP_BASE, + g_param_spec_double("base", + "Base", + "Base of logarithm", + -G_MAXDOUBLE, G_MAXDOUBLE, 10., + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT)); +} + +GType +cmath_clog_get_type(void) +{ + static GType type = 0; + + if(!type) { + static const GTypeInfo info = { + .class_size = sizeof(CMathBaseClass), + .class_init = cmath_clog_class_init, + .instance_size = sizeof(CMathCLog), + }; + type = g_type_register_static(CMATH_BASE_TYPE, "CMathCLog", &info, 0); + } + + return type; +} diff --git a/gstlal-calibration/gst/cmath/cmath_clog10.c b/gstlal-calibration/gst/cmath/cmath_clog10.c new file mode 100644 index 0000000000000000000000000000000000000000..74afcf36141b43af92f1efe808916fae81d37da5 --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_clog10.c @@ -0,0 +1,163 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + + +#include <cmath_base.h> + +#define TYPE_CMATH_CLOG10 \ + (cmath_clog10_get_type()) +#define CMATH_CLOG10(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj),TYPE_CMATH_CLOG10,CMathCLog10)) +#define CMATH_CLOG10_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass),TYPE_CMATH_CLOG10,CMathCLog10Class)) +#define IS_PLUGIN_TEMPLATE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj),TYPE_CMATH_CLOG10)) +#define IS_PLUGIN_TEMPLATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass),TYPE_CMATH_CLOG10)) + +typedef struct _CMathCLog10 CMathCLog10; +typedef struct _CMathCLog10Class CMathCLog10Class; + +GType +cmath_clog10_get_type(void); + +struct _CMathCLog10 +{ + CMathBase cmath_base; +}; + +struct _CMathCLog10Class +{ + CMathBaseClass parent_class; +}; + + +/* + * ============================================================================ + * + * GstBaseTransform vmethod Implementations + * + * ============================================================================ + */ + +/* An in-place transform really does the same thing as the chain function */ + +static GstFlowReturn +transform_ip(GstBaseTransform *trans, GstBuffer *buf) +{ + CMathCLog10* element = CMATH_CLOG10(trans); + int bits = element -> cmath_base.bits; + int is_complex = element -> cmath_base.is_complex; + + /* Debugging + * + * GstObject* element_gstobj = GST_OBJECT(trans); + * int channels = element -> cmath_base.channels; + * int rate = element -> cmath_base.rate; + * g_print("[%s]: passing GstBuffer: ", element_gstobj->name); + * g_print("%d channels, ", channels); + * g_print("%d bits, ", bits); + * g_print("rate: %d, ", rate); + */ + + GstMapInfo info; + if(!gst_buffer_map(buf, &info, GST_MAP_READWRITE)) { + GST_ERROR_OBJECT(trans, "gst_buffer_map failed\n"); + } + gpointer data = info.data; + gpointer data_end = data + info.size; + + if(is_complex == 1) { + + if(bits == 128) { + double complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = clog(*ptr) / log(10); + } + } else if(bits == 64) { + float complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = clogf(*ptr) / log(10); + } + } else { + g_assert_not_reached(); + } + } else if(is_complex == 0) { + + if(bits == 64) { + double *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = log10(*ptr); + } + } else if(bits == 32) { + float *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = log10f(*ptr); + } + } else { + g_assert_not_reached(); + } + } else { + g_assert_not_reached(); + } + gst_buffer_unmap(buf, &info); + return GST_FLOW_OK; +} + + +/* + * ============================================================================ + * + * Type Support + * + * ============================================================================ + */ + +/* Initialize the plugin's class */ +static void +cmath_clog10_class_init(gpointer klass, gpointer klass_data) +{ + GstBaseTransformClass *basetransform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gst_element_class_set_details_simple(GST_ELEMENT_CLASS(klass), + "Logarithm base 10", + "Filter/Audio", + "Calculate logarithm base 10, y = log_10 x", + "Aaron Viets <aaron.viets@ligo.org>, Leo Singer <leo.singer@ligo.org>"); + + basetransform_class -> transform_ip = GST_DEBUG_FUNCPTR(transform_ip); + basetransform_class -> set_caps = GST_DEBUG_FUNCPTR(set_caps); +} + +GType +cmath_clog10_get_type(void) +{ + static GType type = 0; + + if(!type) { + static const GTypeInfo info = { + .class_size = sizeof(CMathBaseClass), + .class_init = cmath_clog10_class_init, + .instance_size = sizeof(CMathBase), + }; + type = g_type_register_static(CMATH_BASE_TYPE, "CMathCLog10", &info, 0); + } + + return type; +} diff --git a/gstlal-calibration/gst/cmath/cmath_cpow.c b/gstlal-calibration/gst/cmath/cmath_cpow.c new file mode 100644 index 0000000000000000000000000000000000000000..aabd6cce3d6699ba9d57f3bb93b8ac6b94fbf269 --- /dev/null +++ b/gstlal-calibration/gst/cmath/cmath_cpow.c @@ -0,0 +1,218 @@ +/* + * Copyright (C) 2010 Leo Singer + * 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. + */ + + +#include <cmath_base.h> + +#define TYPE_CMATH_CPOW \ + (cmath_cpow_get_type()) +#define CMATH_CPOW(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj),TYPE_CMATH_CPOW,CMathCPow)) +#define CMATH_CPOW_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass),TYPE_CMATH_CPOW,CMathCPowClass)) +#define IS_PLUGIN_TEMPLATE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj),TYPE_CMATH_CPOW)) +#define IS_PLUGIN_TEMPLATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass),TYPE_CMATH_CPOW)) + +typedef struct _CMathCPow CMathCPow; +typedef struct _CMathCPowClass CMathCPowClass; + +GType +cmath_cpow_get_type(void); + +struct _CMathCPow +{ + CMathBase cmath_base; + double exponent; +}; + +struct _CMathCPowClass +{ + CMathBaseClass parent_class; +}; + + +/* + * ============================================================================ + * + * GstBaseTransform vmethod Implementations + * + * ============================================================================ + */ + +/* An in-place transform really does the same thing as the chain function */ + +static GstFlowReturn transform_ip(GstBaseTransform *trans, GstBuffer *buf) +{ + CMathCPow* element = CMATH_CPOW(trans); + int bits = element -> cmath_base.bits; + int is_complex = element -> cmath_base.is_complex; + + /* + * Debugging + * + * GstObject* element_gstobj = GST_OBJECT(trans); + * int channels = element -> cmath_base.channels; + * int rate = element -> cmath_base.rate; + * g_print("[%s]: passing GstBuffer: ", element_gstobj->name); + * g_print("%d channels, ", channels); + * g_print("%d bits, ", bits); + * g_print("rate: %d, ", rate); + */ + + GstMapInfo info; + if(!gst_buffer_map(buf, &info, GST_MAP_READWRITE)) { + GST_ERROR_OBJECT(trans, "gst_buffer_map failed\n"); + } + gpointer data = info.data; + gpointer data_end = data + info.size; + + const double n = element -> exponent; + + if(is_complex == 1) { + + if(bits == 128) { + double complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = cpow(*ptr, n); + } + } else if(bits == 64) { + float complex *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = cpowf(*ptr, n); + } + } else { + g_assert_not_reached(); + } + } else if(is_complex == 0) { + + if(bits == 64) { + double *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = pow(*ptr, n); + } + } else if(bits == 32) { + float *ptr, *end = data_end; + for(ptr = data; ptr < end; ptr++) { + *ptr = powf(*ptr, n); + } + } else { + g_assert_not_reached(); + } + } else { + g_assert_not_reached(); + } + gst_buffer_unmap(buf, &info); + return GST_FLOW_OK; +} + + +/* + * ============================================================================ + * + * Type Support + * + * ============================================================================ + */ + +/* Set the exponent */ +enum property +{ + PROP_EXPONENT = 1, +}; + +static void +set_property(GObject * object, enum property id, const GValue * value, + GParamSpec * pspec) +{ + CMathCPow *element = CMATH_CPOW(object); + + GST_OBJECT_LOCK(element); + + switch(id) { + case PROP_EXPONENT: + element->exponent = g_value_get_double(value); + break; + } + + GST_OBJECT_UNLOCK(element); +} + +static void +get_property(GObject * object, enum property id, GValue * value, + GParamSpec * pspec) +{ + CMathCPow *element = CMATH_CPOW(object); + + GST_OBJECT_LOCK(element); + + switch(id) { + case PROP_EXPONENT: + g_value_set_double(value, element->exponent); + break; + } + + GST_OBJECT_UNLOCK(element); +} + +/* Initialize the plugin's class */ +static void +cmath_cpow_class_init(gpointer klass, gpointer klass_data) +{ + GObjectClass *gobject_class = G_OBJECT_CLASS(klass); + GstBaseTransformClass *basetransform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gst_element_class_set_details_simple(GST_ELEMENT_CLASS(klass), + "Raise input to a power", + "Filter/Audio", + "Calculate input raised to the power n, y = x^n", + "Aaron Viets <aaron.viets@ligo.org>, Leo Singer <leo.singer@ligo.org>"); + + basetransform_class -> transform_ip = GST_DEBUG_FUNCPTR(transform_ip); + basetransform_class -> set_caps = GST_DEBUG_FUNCPTR(set_caps); + + gobject_class->get_property = GST_DEBUG_FUNCPTR(get_property); + gobject_class->set_property = GST_DEBUG_FUNCPTR(set_property); + + g_object_class_install_property(gobject_class, + PROP_EXPONENT, + g_param_spec_double("exponent", + "Exponent", + "Exponent", + -G_MAXDOUBLE, G_MAXDOUBLE, 2., + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT)); +} + +GType +cmath_cpow_get_type(void) +{ + static GType type = 0; + + if(!type) { + static const GTypeInfo info = { + .class_size = sizeof(CMathBaseClass), + .class_init = cmath_cpow_class_init, + .instance_size = sizeof(CMathCPow), + }; + type = g_type_register_static(CMATH_BASE_TYPE, "CMathCPow", &info, 0); + } + + return type; +} diff --git a/gstlal-calibration/gst/lal/Makefile.am b/gstlal-calibration/gst/lal/Makefile.am index 180e52ffbd312e99285184205eb35de6c8a0b021..f21b69651e741e226e9ad9e378c0f3ec5512246b 100644 --- a/gstlal-calibration/gst/lal/Makefile.am +++ b/gstlal-calibration/gst/lal/Makefile.am @@ -14,7 +14,8 @@ lib@GSTPLUGINPREFIX@gstlalcalibration_la_SOURCES = \ gstlal_demodulate.c gstlal_demodulate.h \ gstlal_insertgap.c gstlal_insertgap.h \ gstlal_fccupdate.c gstlal_fccupdate.h \ - gstlal_transferfunction.c gstlal_transferfunction.h + gstlal_transferfunction.c gstlal_transferfunction.h \ + gstlal_trackfrequency.c gstlal_trackfrequency.h lib@GSTPLUGINPREFIX@gstlalcalibration_la_CPPFLAGS = $(AM_CPPFLAGS) $(PYTHON_CPPFLAGS) lib@GSTPLUGINPREFIX@gstlalcalibration_la_CFLAGS = $(AM_CFLAGS) $(LAL_CFLAGS) $(GSTLAL_CFLAGS) $(gstreamer_CFLAGS) $(gstreamer_audio_CFLAGS) lib@GSTPLUGINPREFIX@gstlalcalibration_la_LDFLAGS = $(AM_LDFLAGS) $(LAL_LIBS) $(GSTLAL_LIBS) $(PYTHON_LIBS) $(gstreamer_LIBS) $(gstreamer_audio_LIBS) $(GSTLAL_PLUGIN_LDFLAGS) diff --git a/gstlal-calibration/gst/lal/gstlal_trackfrequency.c b/gstlal-calibration/gst/lal/gstlal_trackfrequency.c new file mode 100644 index 0000000000000000000000000000000000000000..50a341adc2e4620fd2e9ba6dd6778424e7b4faf8 --- /dev/null +++ b/gstlal-calibration/gst/lal/gstlal_trackfrequency.c @@ -0,0 +1,813 @@ +/* + * Copyright (C) 2018 Aaron Viets <aaron.viets@ligo.org> + * + * 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 + * + * ============================================================================ + */ + + +/* + * stuff from the C library + */ + +#include <string.h> +#include <complex.h> +#include <math.h> + + +/* + * stuff from gobject/gstreamer + */ + + +#include <glib.h> +#include <gst/gst.h> +#include <gst/audio/audio.h> +#include <gst/base/gstbasetransform.h> +#include <gst/audio/audio.h> + + +/* + * our own stuff + */ + + +#include <gstlal/gstlal_audio_info.h> +#include <gstlal_trackfrequency.h> + + +/* + * ============================================================================ + * + * Utilities + * + * ============================================================================ + */ + + +static void update_frequency(double *current_frequency, guint64 *crossover_times, guint64 num_halfcycles, guint64 *num_stored, guint64 new_crossover_time) { + + /* First, update the recent times that +/- transitions have occurred */ + if(*num_stored <= num_halfcycles) { + crossover_times[*num_stored] = new_crossover_time; + (*num_stored)++; + } else { + guint64 i; + for(i = 0; i < num_halfcycles; i++) + crossover_times[i] = crossover_times[i + 1]; + crossover_times[num_halfcycles] = new_crossover_time; + } + + /* Now, update the frequency */ + if(*num_stored > 1) { + guint64 time_elapsed = new_crossover_time - *crossover_times; + *current_frequency = 1000000000.0 / (2.0 * (double) time_elapsed / (double) (*num_stored - 1)); + } +} + + +#define DEFINE_TRACK_FREQUENCY(DTYPE) \ +static void trackfrequency_ ## DTYPE(const DTYPE *src, DTYPE *dst, gint64 size, int rate, guint64 pts, guint64 ets, double *current_frequency, guint64 *crossover_times, guint64 num_halfcycles, int *sign, gint64 *check_step, guint64 *num_stored) { \ + \ + /* Check if we have information about previous data. If not, test whether the first sample is positive or negative */ \ + if(*sign == 0) { \ + if(*src < 0) \ + *sign = -1; \ + else \ + *sign = 1; \ + } \ + \ + gint64 i = 0; \ + gint64 j = 0; \ + while(i < size - 1) { \ + \ + while(*check_step) { \ + if(*sign * src[i] < 0) { \ + *check_step = -(*check_step) / 2; \ + *sign = -(*sign); \ + } \ + /* We don't want to fall off the edge of the buffer */ \ + if(*check_step < -i) \ + *check_step = -i / 2; \ + if(*check_step >= size - i) \ + *check_step = (size - i) / 2; \ + i += *check_step; \ + } \ + \ + /* At this point, we could be on either side of a +/- transition, and/or at the end of the buffer */ \ + if(i == 0) { \ + /* There is a transition at the beginning of a buffer */ \ + *dst = (DTYPE) *current_frequency; \ + j++; \ + update_frequency(current_frequency, crossover_times, num_halfcycles, num_stored, pts); \ + } else if(src[i] * src[i - 1] < 0) { \ + /* We are just after a transition (and possibly at the end of the buffer) */ \ + while(j <= i) { \ + dst[j] = (DTYPE) *current_frequency; \ + j++; \ + } \ + update_frequency(current_frequency, crossover_times, num_halfcycles, num_stored, pts + gst_util_uint64_scale_int_round(i, GST_SECOND, rate)); \ + } else if(i != size - 1 && src[i] * src[i + 1] < 0) { \ + /* We are before a transition */ \ + i++; \ + *sign = -(*sign); \ + while(j <= i) { \ + dst[j] = (DTYPE) *current_frequency; \ + j++; \ + } \ + update_frequency(current_frequency, crossover_times, num_halfcycles, num_stored, pts + gst_util_uint64_scale_int_round(i, GST_SECOND, rate)); \ + } else { \ + /* We are at the end of the buffer, and there is no transition here */ \ + while(j <= i) { \ + dst[j] = (DTYPE) *current_frequency; \ + j++; \ + } \ + } \ + \ + if(*current_frequency) \ + *check_step = (int) (0.2 * rate / *current_frequency + 0.5); \ + else \ + *check_step = 1; \ + } \ + /* We should be at the end of the output buffer now */ \ + g_assert_cmpint(size, == , j); \ + \ + /* Set the step size for the start of the next buffer */ \ + if(*current_frequency) { \ + double ETA_from_next_pts = 1.0 / *current_frequency - (ets - crossover_times[*num_stored - 1]) / 1000000000.0; \ + if(ETA_from_next_pts < 0.0) \ + ETA_from_next_pts = 0; \ + *check_step = (int) (0.2 * rate * ETA_from_next_pts + 1.0); \ + } else \ + *check_step = 1; \ +} + +DEFINE_TRACK_FREQUENCY(float) +DEFINE_TRACK_FREQUENCY(double) + + +#define DEFINE_TRACK_FREQUENCY_COMPLEX(DTYPE, F_OR_BLANK) \ +static void trackfrequency_complex_ ## DTYPE(const DTYPE complex *src, DTYPE *dst, gint64 size, int rate, guint64 pts, guint64 ets, double *current_frequency, guint64 *crossover_times, guint64 num_halfcycles, int *sign, gint64 *check_step, guint64 *num_stored) { \ + \ + /* Check if we have information about previous data. If not, test whether the first sample is positive or negative */ \ + if(*sign == 0) { \ + if(creal ## F_OR_BLANK(*src) < 0) \ + *sign = -1; \ + else \ + *sign = 1; \ + } \ + \ + gint64 i = 0; \ + gint64 j = 0; \ + while(i < size - 1) { \ + \ + while(*check_step) { \ + if(*sign * creal ## F_OR_BLANK(src[i]) < 0) { \ + *check_step = -(*check_step) / 2; \ + *sign = -(*sign); \ + } \ + /* We don't want to fall off the edge of the buffer */ \ + if(*check_step < -i) \ + *check_step = -i / 2; \ + if(*check_step >= size - i) \ + *check_step = (size - i) / 2; \ + i += *check_step; \ + } \ + \ + /* At this point, we could be on either side of a +/- transition, and/or at the end of the buffer */ \ + if(i == 0) { \ + /* There is a transition at the beginning of a buffer */ \ + *dst = (DTYPE) *current_frequency; \ + j++; \ + update_frequency(current_frequency, crossover_times, num_halfcycles, num_stored, pts); \ + /* Check if the frequency is negative */ \ + if(creal ## F_OR_BLANK(src[i]) * cimag ## F_OR_BLANK(src[i]) > 0) \ + *current_frequency = -(*current_frequency); \ + } else if(creal ## F_OR_BLANK(src[i]) * creal ## F_OR_BLANK(src[i - 1]) < 0) { \ + /* We are just after a transition (and possibly at the end of the buffer) */ \ + while(j <= i) { \ + dst[j] = (DTYPE) *current_frequency; \ + j++; \ + } \ + update_frequency(current_frequency, crossover_times, num_halfcycles, num_stored, pts + gst_util_uint64_scale_int_round(i, GST_SECOND, rate)); \ + /* Check if the frequency is negative */ \ + if(creal ## F_OR_BLANK(src[i]) * cimag ## F_OR_BLANK(src[i]) > 0) \ + *current_frequency = -(*current_frequency); \ + } else if(i != size - 1 && creal ## F_OR_BLANK(src[i]) * creal ## F_OR_BLANK(src[i + 1]) < 0) { \ + /* We are before a transition */ \ + i++; \ + *sign = -(*sign); \ + while(j <= i) { \ + dst[j] = (DTYPE) *current_frequency; \ + j++; \ + } \ + update_frequency(current_frequency, crossover_times, num_halfcycles, num_stored, pts + gst_util_uint64_scale_int_round(i, GST_SECOND, rate)); \ + /* Check if the frequency is negative */ \ + if(creal ## F_OR_BLANK(src[i]) * cimag ## F_OR_BLANK(src[i]) > 0) \ + *current_frequency = -(*current_frequency); \ + } else { \ + /* We are at the end of the buffer, and there is no transition here */ \ + while(j <= i) { \ + dst[j] = (DTYPE) *current_frequency; \ + j++; \ + } \ + } \ + \ + if(*current_frequency) \ + *check_step = (int) (0.2 * rate / *current_frequency + 0.5); \ + else \ + *check_step = 1; \ + } \ + /* We should be at the end of the output buffer now */ \ + g_assert_cmpint(size, == , j); \ + \ + /* Set the step size for the start of the next buffer */ \ + if(*current_frequency) { \ + double ETA_from_next_pts = 1.0 / *current_frequency - (ets - crossover_times[*num_stored - 1]) / 1000000000.0; \ + if(ETA_from_next_pts < 0.0) \ + ETA_from_next_pts = 0; \ + *check_step = (int) (0.2 * rate * ETA_from_next_pts + 0.5); \ + } else \ + *check_step = 1; \ +} + +DEFINE_TRACK_FREQUENCY_COMPLEX(float, f) +DEFINE_TRACK_FREQUENCY_COMPLEX(double, ) + + +static void trackfrequency(const void *src, void *dst, guint64 src_size, guint64 pts, guint64 ets, GSTLALTrackFrequency *element) { + + switch(element->data_type) { + case GSTLAL_TRACKFREQUENCY_F32: + trackfrequency_float(src, dst, (gint64) src_size / element->unit_size, element->rate, pts, ets, &element->current_frequency, element->crossover_times, element->num_halfcycles, &element->sign, &element->check_step, &element->num_stored); + break; + case GSTLAL_TRACKFREQUENCY_F64: + trackfrequency_double(src, dst, (gint64) src_size / element->unit_size, element->rate, pts, ets, &element->current_frequency, element->crossover_times, element->num_halfcycles, &element->sign, &element->check_step, &element->num_stored); + break; + case GSTLAL_TRACKFREQUENCY_Z64: + trackfrequency_complex_float(src, dst, (gint64) src_size / element->unit_size, element->rate, pts, ets, &element->current_frequency, element->crossover_times, element->num_halfcycles, &element->sign, &element->check_step, &element->num_stored); + break; + case GSTLAL_TRACKFREQUENCY_Z128: + trackfrequency_complex_double(src, dst, (gint64) src_size / element->unit_size, element->rate, pts, ets, &element->current_frequency, element->crossover_times, element->num_halfcycles, &element->sign, &element->check_step, &element->num_stored); + break; + default: + g_assert_not_reached(); + } +} + + +/* + * set the metadata on an output buffer + */ + + +static void set_metadata(GSTLALTrackFrequency *element, GstBuffer *buf, guint64 outsamples, gboolean gap) +{ + if(element->data_type == GSTLAL_TRACKFREQUENCY_Z64 || element->data_type == GSTLAL_TRACKFREQUENCY_Z128) + outsamples *= 2; + GST_BUFFER_OFFSET(buf) = element->next_out_offset; + element->next_out_offset += outsamples; + GST_BUFFER_OFFSET_END(buf) = element->next_out_offset; + GST_BUFFER_TIMESTAMP(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET(buf) - element->offset0, GST_SECOND, element->rate); + GST_BUFFER_DURATION(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, element->rate) - GST_BUFFER_TIMESTAMP(buf); + if(G_UNLIKELY(element->need_discont)) { + GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_DISCONT); + element->need_discont = FALSE; + } + if(gap) + GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_GAP); + else + GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP); +} + + +/* + * ============================================================================ + * + * GStreamer Boiler Plate + * + * ============================================================================ + */ + + +#define INCAPS \ + "audio/x-raw, " \ + "format = (string) {"GST_AUDIO_NE(F32)", "GST_AUDIO_NE(F64)", "GST_AUDIO_NE(Z64)", "GST_AUDIO_NE(Z128)"}, " \ + "rate = " GST_AUDIO_RATE_RANGE ", " \ + "channels = (int) 1, " \ + "layout = (string) interleaved, " \ + "channel-mask = (bitmask) 0" + +#define OUTCAPS \ + "audio/x-raw, " \ + "format = (string) {"GST_AUDIO_NE(F32)", "GST_AUDIO_NE(F64)"}, " \ + "rate = " GST_AUDIO_RATE_RANGE ", " \ + "channels = (int) 1, " \ + "layout = (string) interleaved, " \ + "channel-mask = (bitmask) 0" + + +static GstStaticPadTemplate sink_factory = GST_STATIC_PAD_TEMPLATE( + GST_BASE_TRANSFORM_SINK_NAME, + GST_PAD_SINK, + GST_PAD_ALWAYS, + GST_STATIC_CAPS(INCAPS) +); + + +static GstStaticPadTemplate src_factory = GST_STATIC_PAD_TEMPLATE( + GST_BASE_TRANSFORM_SRC_NAME, + GST_PAD_SRC, + GST_PAD_ALWAYS, + GST_STATIC_CAPS(OUTCAPS) +); + + +G_DEFINE_TYPE( + GSTLALTrackFrequency, + gstlal_trackfrequency, + GST_TYPE_BASE_TRANSFORM +); + + +/* + * ============================================================================ + * + * GstBaseTransform Method Overrides + * + * ============================================================================ + */ + + +/* + * get_unit_size() + */ + + +static gboolean get_unit_size(GstBaseTransform *trans, GstCaps *caps, gsize *size) +{ + GstAudioInfo info; + gboolean success = gstlal_audio_info_from_caps(&info, caps); + if(success) + *size = GST_AUDIO_INFO_BPF(&info); + else + GST_WARNING_OBJECT(trans, "unable to parse caps %" GST_PTR_FORMAT, caps); + return success; +} + + +/* + * transform_caps() + */ + + +static GstCaps *transform_caps(GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps, GstCaps *filter) +{ + guint n; + + caps = gst_caps_normalize(gst_caps_copy(caps)); + GstCaps *othercaps = gst_caps_new_empty(); + + switch(direction) { + case GST_PAD_SRC: + /* There are two possible sink pad formats for each src pad format, so the sink pad has twice as many caps structures */ + for(n = 0; n < gst_caps_get_size(caps); n++) { + gst_caps_append(othercaps, gst_caps_copy_nth(caps, n)); + gst_caps_append(othercaps, gst_caps_copy_nth(caps, n)); + + GstStructure *str = gst_caps_get_structure(othercaps, 2 * n); + const gchar *format = gst_structure_get_string(str, "format"); + + if(!format) { + GST_DEBUG_OBJECT(trans, "unrecognized caps %" GST_PTR_FORMAT, othercaps); + goto error; + } else if(!strcmp(format, GST_AUDIO_NE(F32))) + gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(Z64), NULL); + else if(!strcmp(format, GST_AUDIO_NE(F64))) + gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(Z128), NULL); + else { + GST_DEBUG_OBJECT(trans, "unrecognized format %s in %" GST_PTR_FORMAT, format, othercaps); + goto error; + } + } + break; + + case GST_PAD_SINK: + for(n = 0; n < gst_caps_get_size(caps); n++) { + gst_caps_append(othercaps, gst_caps_copy_nth(caps, n)); + GstStructure *str = gst_caps_get_structure(othercaps, n); + const gchar *format = gst_structure_get_string(str, "format"); + + if(!format) { + GST_DEBUG_OBJECT(trans, "unrecognized caps %" GST_PTR_FORMAT, othercaps); + goto error; + } else if(!strcmp(format, GST_AUDIO_NE(F32)) || !strcmp(format, GST_AUDIO_NE(Z64))) + gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(F32), NULL); + else if(!strcmp(format, GST_AUDIO_NE(F64)) || !strcmp(format, GST_AUDIO_NE(Z128))) + gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(F64), NULL); + else { + GST_DEBUG_OBJECT(trans, "unrecognized format %s in %" GST_PTR_FORMAT, format, othercaps); + goto error; + } + } + break; + + case GST_PAD_UNKNOWN: + GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid direction GST_PAD_UNKNOWN")); + goto error; + } + + if(filter) { + GstCaps *intersection = gst_caps_intersect(othercaps, filter); + gst_caps_unref(othercaps); + othercaps = intersection; + } + gst_caps_unref(caps); + return gst_caps_simplify(othercaps); + +error: + gst_caps_unref(caps); + gst_caps_unref(othercaps); + return GST_CAPS_NONE; +} + + +/* + * set_caps() + */ + + +static gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, GstCaps *outcaps) +{ + GSTLALTrackFrequency *element = GSTLAL_TRACKFREQUENCY(trans); + gint rate_in, rate_out; + gsize unit_size; + const gchar *format = gst_structure_get_string(gst_caps_get_structure(incaps, 0), "format"); + + /* + * parse the caps + */ + + if(!format) { + GST_DEBUG_OBJECT(element, "unable to parse format from %" GST_PTR_FORMAT, incaps); + return FALSE; + } + if(!get_unit_size(trans, incaps, &unit_size)) { + GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed"); + return FALSE; + } + if(!gst_structure_get_int(gst_caps_get_structure(incaps, 0), "rate", &rate_in)) { + GST_DEBUG_OBJECT(element, "unable to parse rate from %" GST_PTR_FORMAT, incaps); + return FALSE; + } + if(!gst_structure_get_int(gst_caps_get_structure(outcaps, 0), "rate", &rate_out)) { + GST_DEBUG_OBJECT(element, "unable to parse rate from %" GST_PTR_FORMAT, outcaps); + return FALSE; + } + + /* + * require the output rate to be equal to the input rate + */ + + if(rate_out != rate_in) { + GST_ERROR_OBJECT(element, "output rate is not equal to input rate. input caps = %" GST_PTR_FORMAT " output caps = %" GST_PTR_FORMAT, incaps, outcaps); + return FALSE; + } + + /* + * record stream parameters + */ + + if(!strcmp(format, GST_AUDIO_NE(F32))) { + element->data_type = GSTLAL_TRACKFREQUENCY_F32; + g_assert_cmpuint(unit_size, ==, 4); + } else if(!strcmp(format, GST_AUDIO_NE(F64))) { + element->data_type = GSTLAL_TRACKFREQUENCY_F64; + g_assert_cmpuint(unit_size, ==, 8); + } else if(!strcmp(format, GST_AUDIO_NE(Z64))) { + element->data_type = GSTLAL_TRACKFREQUENCY_Z64; + g_assert_cmpuint(unit_size, ==, 8); + } else if(!strcmp(format, GST_AUDIO_NE(Z128))) { + element->data_type = GSTLAL_TRACKFREQUENCY_Z128; + g_assert_cmpuint(unit_size, ==, 16); + } else + g_assert_not_reached(); + + element->rate = rate_in; + element->unit_size = unit_size; + + return TRUE; +} + + +/* + * transform_size() + */ + + +static gboolean transform_size(GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps, gsize size, GstCaps *othercaps, gsize *othersize) +{ + GSTLALTrackFrequency *element = GSTLAL_TRACKFREQUENCY(trans); + + gsize unit_size; + + switch(direction) { + case GST_PAD_SRC: + /* + * We have the size of the output buffer, and we set the size of the input buffer, + * which is half as large in bytes. + */ + + if(!element->data_type) { + GST_DEBUG_OBJECT(element, "Data type is not set. Cannot specify incoming buffer size given outgoing buffer size."); + return FALSE; + } + + if(!get_unit_size(trans, caps, &unit_size)) { + GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed"); + return FALSE; + } + + /* buffer size in bytes should be a multiple of unit_size in bytes */ + if(G_UNLIKELY(size % unit_size)) { + GST_DEBUG_OBJECT(element, "buffer size %" G_GSIZE_FORMAT " is not a multiple of %" G_GSIZE_FORMAT, size, unit_size); + return FALSE; + } + + if(element->data_type == GSTLAL_TRACKFREQUENCY_F32 || element->data_type == GSTLAL_TRACKFREQUENCY_F64) + *othersize = size; + else if(element->data_type == GSTLAL_TRACKFREQUENCY_Z64 || element->data_type == GSTLAL_TRACKFREQUENCY_Z128) + *othersize = size * 2; + else + g_assert_not_reached(); + + break; + + case GST_PAD_SINK: + /* + * We have the size of the input buffer, and we set the size of the output buffer, + * which is twice as large in bytes. + */ + + if(!get_unit_size(trans, caps, &unit_size)) { + GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed"); + return FALSE; + } + + /* buffer size in bytes should be a multiple of unit_size in bytes */ + if(G_UNLIKELY(size % unit_size)) { + GST_DEBUG_OBJECT(element, "buffer size %" G_GSIZE_FORMAT " is not a multiple of %" G_GSIZE_FORMAT, size, unit_size); + return FALSE; + } + + GstStructure *str = gst_caps_get_structure(caps, 0); + const gchar *format = gst_structure_get_string(str, "format"); + + if(!format) { + GST_DEBUG_OBJECT(trans, "unrecognized caps %" GST_PTR_FORMAT, caps); + return FALSE; + } else if(!strcmp(format, GST_AUDIO_NE(F32)) || !strcmp(format, GST_AUDIO_NE(F64))) + *othersize = size; + else if(!strcmp(format, GST_AUDIO_NE(Z64)) || !strcmp(format, GST_AUDIO_NE(Z128))) + *othersize = size / 2; + else { + GST_DEBUG_OBJECT(trans, "unrecognized format %s in %" GST_PTR_FORMAT, format, caps); + return FALSE; + } + + if(!get_unit_size(trans, caps, &unit_size)) { + GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed"); + return FALSE; + } + break; + + case GST_PAD_UNKNOWN: + GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid direction GST_PAD_UNKNOWN")); + return FALSE; + } + + return TRUE; +} + + +/* + * start() + */ + + +static gboolean start(GstBaseTransform *trans) +{ + GSTLALTrackFrequency *element = GSTLAL_TRACKFREQUENCY(trans); + + element->t0 = GST_CLOCK_TIME_NONE; + element->offset0 = GST_BUFFER_OFFSET_NONE; + element->next_in_offset = GST_BUFFER_OFFSET_NONE; + element->next_out_offset = GST_BUFFER_OFFSET_NONE; + element->need_discont = TRUE; + + return TRUE; +} + + +/* + * transform() + */ + + +static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuffer *outbuf) +{ + GSTLALTrackFrequency *element = GSTLAL_TRACKFREQUENCY(trans); + GstMapInfo inmap, outmap; + GstFlowReturn result = GST_FLOW_OK; + + /* + * check for discontinuity + */ + + if(G_UNLIKELY(GST_BUFFER_IS_DISCONT(inbuf) || GST_BUFFER_OFFSET(inbuf) != element->next_in_offset || !GST_CLOCK_TIME_IS_VALID(element->t0))) { + element->t0 = GST_BUFFER_PTS(inbuf); + element->offset0 = element->next_out_offset = GST_BUFFER_OFFSET(inbuf); + element->need_discont = TRUE; + element->crossover_times = g_malloc((element->num_halfcycles + 1) * sizeof(guint64)); + element->num_stored = 0; + element->sign = 0; + element->check_step = 1; + } + element->next_in_offset = GST_BUFFER_OFFSET_END(inbuf); + + /* + * process buffer + */ + + if(!GST_BUFFER_FLAG_IS_SET(inbuf, GST_BUFFER_FLAG_GAP)) { + + /* + * input is not gap. + */ + + gst_buffer_map(inbuf, &inmap, GST_MAP_READ); + gst_buffer_map(outbuf, &outmap, GST_MAP_WRITE); + trackfrequency(inmap.data, outmap.data, inmap.size, GST_BUFFER_PTS(inbuf), GST_BUFFER_PTS(inbuf) + GST_BUFFER_DURATION(inbuf), element); + set_metadata(element, outbuf, outmap.size / element->unit_size, FALSE); + gst_buffer_unmap(outbuf, &outmap); + gst_buffer_unmap(inbuf, &inmap); + } else { + + /* + * input is gap. + */ + + element->num_stored = 0; + element->sign = 0; + GST_BUFFER_FLAG_SET(outbuf, GST_BUFFER_FLAG_GAP); + gst_buffer_map(outbuf, &outmap, GST_MAP_WRITE); + memset(outmap.data, 0, outmap.size); + set_metadata(element, outbuf, outmap.size / element->unit_size, TRUE); + gst_buffer_unmap(outbuf, &outmap); + } + + /* + * done + */ + + return result; +} + + +/* + * ============================================================================ + * + * GObject Method Overrides + * + * ============================================================================ + */ + + +/* + * properties + */ + + +enum property { + ARG_NUM_HALFCYCLES = 1, +}; + + +static void set_property(GObject *object, enum property prop_id, const GValue *value, GParamSpec *pspec) +{ + GSTLALTrackFrequency *element = GSTLAL_TRACKFREQUENCY(object); + + GST_OBJECT_LOCK(element); + + switch (prop_id) { + case ARG_NUM_HALFCYCLES: + element->num_halfcycles = g_value_get_uint64(value); + break; + default: + G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec); + break; + } + + GST_OBJECT_UNLOCK(element); +} + + +static void get_property(GObject *object, enum property prop_id, GValue *value, GParamSpec *pspec) +{ + GSTLALTrackFrequency *element = GSTLAL_TRACKFREQUENCY(object); + + GST_OBJECT_LOCK(element); + + switch (prop_id) { + case ARG_NUM_HALFCYCLES: + g_value_set_uint64(value, element->num_halfcycles); + break; + default: + G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec); + break; + } + + GST_OBJECT_UNLOCK(element); +} + +/* + * class_init() + */ + + +static void gstlal_trackfrequency_class_init(GSTLALTrackFrequencyClass *klass) +{ + GObjectClass *gobject_class = G_OBJECT_CLASS(klass); + GstElementClass *element_class = GST_ELEMENT_CLASS(klass); + GstBaseTransformClass *transform_class = GST_BASE_TRANSFORM_CLASS(klass); + + gobject_class->set_property = GST_DEBUG_FUNCPTR(set_property); + gobject_class->get_property = GST_DEBUG_FUNCPTR(get_property); + + transform_class->transform_caps = GST_DEBUG_FUNCPTR(transform_caps); + transform_class->get_unit_size = GST_DEBUG_FUNCPTR(get_unit_size); + transform_class->set_caps = GST_DEBUG_FUNCPTR(set_caps); + transform_class->transform_size = GST_DEBUG_FUNCPTR(transform_size); + transform_class->start = GST_DEBUG_FUNCPTR(start); + transform_class->transform = GST_DEBUG_FUNCPTR(transform); + + gst_element_class_set_details_simple(element_class, + "TrackFrequency", + "Filter/Audio", + "Attempts to measure the loudest frequency of a signal.", + "Aaron Viets <aaron.viets@ligo.org>" + ); + + gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&src_factory)); + gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&sink_factory)); + + g_object_class_install_property( + gobject_class, + ARG_NUM_HALFCYCLES, + g_param_spec_uint64( + "num-halfcycles", + "Number of half-cycles", + "The number of half-periods of a wave to use to compute the frequency.", + 1, G_MAXUINT64, 64, + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT + ) + ); +} + + +/* + * init() + */ + + +static void gstlal_trackfrequency_init(GSTLALTrackFrequency *element) +{ + element->rate = 0; + element->unit_size = 0; + element->current_frequency = 0.0; + element->check_step = 1; + element->num_stored = 0; + element->sign = 0; + gst_base_transform_set_qos_enabled(GST_BASE_TRANSFORM(element), TRUE); + gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(element), TRUE); +} diff --git a/gstlal-calibration/gst/lal/gstlal_trackfrequency.h b/gstlal-calibration/gst/lal/gstlal_trackfrequency.h new file mode 100644 index 0000000000000000000000000000000000000000..abc422b2cc6fea08a61176a6071fcfeadf3f9fa5 --- /dev/null +++ b/gstlal-calibration/gst/lal/gstlal_trackfrequency.h @@ -0,0 +1,109 @@ +/* + * Copyright (C) 2018 Aaron Viets <aaron.viets@ligo.org> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * + * 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., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + + +#ifndef __GST_LAL_TRACKFREQUENCY_H__ +#define __GST_LAL_TRACKFREQUENCY_H__ + + +#include <glib.h> +#include <gst/gst.h> +#include <gst/base/gstbasetransform.h> + + +G_BEGIN_DECLS +#define GSTLAL_TRACKFREQUENCY_TYPE \ + (gstlal_trackfrequency_get_type()) +#define GSTLAL_TRACKFREQUENCY(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj), GSTLAL_TRACKFREQUENCY_TYPE, GSTLALTrackFrequency)) +#define GSTLAL_TRACKFREQUENCY_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass), GSTLAL_TRACKFREQUENCY_TYPE, GSTLALTrackFrequencyClass)) +#define GST_IS_GSTLAL_TRACKFREQUENCY(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj), GSTLAL_TRACKFREQUENCY_TYPE)) +#define GST_IS_GSTLAL_TRACKFREQUENCY_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass), GSTLAL_TRACKFREQUENCY_TYPE)) + + +typedef struct _GSTLALTrackFrequency GSTLALTrackFrequency; +typedef struct _GSTLALTrackFrequencyClass GSTLALTrackFrequencyClass; + + +/* + * GSTLALTrackFrequency: + */ + + +struct _GSTLALTrackFrequency { + GstBaseTransform element; + + /* stream info */ + gint unit_size; + gint rate; + enum gstlal_trackfrequency_data_type { + GSTLAL_TRACKFREQUENCY_F32 = 0, + GSTLAL_TRACKFREQUENCY_F64, + GSTLAL_TRACKFREQUENCY_Z64, + GSTLAL_TRACKFREQUENCY_Z128 + } data_type; + + /* Filter memory */ + guint64 *crossover_times; + gint64 check_step; + double current_frequency; + guint64 num_stored; + int sign; + + /* timestamp book-keeping */ + GstClockTime t0; + guint64 offset0; + guint64 next_in_offset; + guint64 next_out_offset; + gboolean need_discont; + + /* properties */ + guint64 num_halfcycles; +}; + + +/* + * GSTLALTrackFrequencyClass: + * @parent_class: the parent class + */ + + +struct _GSTLALTrackFrequencyClass { + GstBaseTransformClass parent_class; +}; + + +GType gstlal_trackfrequency_get_type(void); + + +G_END_DECLS + + +#endif /* __GST_LAL_TRACKFREQUENCY_H__ */ diff --git a/gstlal-calibration/gst/lal/gstlalcalibration.c b/gstlal-calibration/gst/lal/gstlalcalibration.c index 2fb51e570e5b3cff1969598270fe1fa56f43ba26..f29298198eec58a654c34300e5130c3b9f5fad87 100644 --- a/gstlal-calibration/gst/lal/gstlalcalibration.c +++ b/gstlal-calibration/gst/lal/gstlalcalibration.c @@ -63,6 +63,7 @@ #include <gstlal_insertgap.h> #include <gstlal_fccupdate.h> #include <gstlal_transferfunction.h> +#include <gstlal_trackfrequency.h> /* @@ -93,6 +94,7 @@ static gboolean plugin_init(GstPlugin *plugin) {"lal_insertgap", GSTLAL_INSERTGAP_TYPE}, {"lal_fcc_update", GSTLAL_FCC_UPDATE_TYPE}, {"lal_transferfunction", GSTLAL_TRANSFERFUNCTION_TYPE}, + {"lal_trackfrequency", GSTLAL_TRACKFREQUENCY_TYPE}, {NULL, 0}, }; diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py index 8532ff8d49a8231f61f8a15685cffedcb6d9947a..75de467aece7a72b7d60988b525e985e0ed82d05 100644 --- a/gstlal-calibration/python/calibration_parts.py +++ b/gstlal-calibration/python/calibration_parts.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright (C) 2015 Madeline Wade +# Copyright (C) 2015 Madeline Wade, 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 @@ -176,13 +176,14 @@ def remove_harmonics(pipeline, signal, f0, num_harmonics, f0_var, filter_latency # remove any line(s) from a spectrum. filter length for demodulation (given in seconds) is adjustable # function argument caps must be complex caps + filter_param = 0.0625 head = pipeparts.mktee(pipeline, head) elem = pipeparts.mkgeneric(pipeline, None, "lal_adder", sync = True) mkqueue(pipeline, head).link(elem) for i in range(1, num_harmonics + 1): line = pipeparts.mkgeneric(pipeline, head, "lal_demodulate", line_frequency = i * f0) line = mkresample(pipeline, line, 5, filter_latency == 0, compute_rate) - line_in_witness = lowpass(pipeline, line_in_witness, compute_rate, length = 0.0625 / (f0_var * i), fcut = 0, filter_latency = filter_latency) + line_in_witness = lowpass(pipeline, line_in_witness, compute_rate, length = filter_param / (f0_var * i), fcut = 0, filter_latency = filter_latency) line = mkresample(pipeline, line, 3, filter_latency == 0.0, rate_out) line = pipeparts.mkgeneric(pipeline, line, "lal_demodulate", line_frequency = -1.0 * i * f0, prefactor_real = -2.0) real, imag = split_into_real(pipeline, line) @@ -195,22 +196,53 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, # remove line(s) from a spectrum. filter length for demodulation (given in seconds) is adjustable # function argument caps must be complex caps - if num_harmonics > 1: - witness = pipeparts.mktee(pipeline, witness) + filter_param = 0.0625 + downsample_quality = 5 + upsample_quality = 4 + resample_shift = 96.0 + 16.5 + zero_latency = filter_latency == 0.0 + witness = pipeparts.mktee(pipeline, witness) signal = pipeparts.mktee(pipeline, signal) 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", frequency = 0.0) + f0_beat_frequency = pipeparts.mktee(pipeline, f0_beat_frequency) + for i in range(1, num_harmonics + 1): - # Find amplitude and phase of line in witness channel - line_in_witness = pipeparts.mkgeneric(pipeline, witness, "lal_demodulate", line_frequency = i * f0) - line_in_witness = mkresample(pipeline, line_in_witness, 5, filter_latency == 0.0, compute_rate) - line_in_witness = lowpass(pipeline, line_in_witness, compute_rate, length = 0.0625 / (f0_var * i), fcut = 0, filter_latency = filter_latency) - line_in_witness = pipeparts.mktee(pipeline, line_in_witness) + # Length of low-pass filter + filter_length = filter_param / (f0_var * i) + 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 + + # Find phase shift due to timestamp shift for each harmonic + phase_shift = complex_audioamplify(pipeline, f0_beat_frequency, 0, 2 * i * numpy.pi * time_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 line in signal line_in_signal = pipeparts.mkgeneric(pipeline, signal, "lal_demodulate", line_frequency = i * f0) - line_in_signal = mkresample(pipeline, line_in_signal, 5, filter_latency == 0.0, compute_rate) - line_in_signal = lowpass(pipeline, line_in_signal, compute_rate, length = 0.0625 / (f0_var * i), fcut = 0, filter_latency = filter_latency) + 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) # Find transfer function between witness channel and signal at this frequency tf_at_f = complex_division(pipeline, line_in_signal, line_in_witness) @@ -221,8 +253,8 @@ 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)) - reconstructed_line_in_signal = mkresample(pipeline, reconstructed_line_in_signal, 3, filter_latency == 0.0, rate_out) + 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, imag = split_into_real(pipeline, reconstructed_line_in_signal) pipeparts.mkfakesink(pipeline, imag) @@ -253,7 +285,7 @@ def lowpass(pipeline, head, rate, length = 1.0, fcut = 500, filter_latency = 0.5 lowpass /= numpy.sum(lowpass) # Now apply the filter - return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * filter_latency), fir_matrix = [lowpass], time_domain = td) + return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * filter_latency + 0.5), fir_matrix = [lowpass], time_domain = td) def highpass(pipeline, head, rate, length = 1.0, fcut = 9.0, filter_latency = 0.5, td = True): length = int(length * rate) @@ -270,7 +302,7 @@ def highpass(pipeline, head, rate, length = 1.0, fcut = 9.0, filter_latency = 0. highpass[int((length - 1) / 2)] += 1 # Now apply the filter - return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * filter_latency), fir_matrix = [highpass], time_domain = td) + return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * filter_latency + 0.5), fir_matrix = [highpass], time_domain = td) def bandpass(pipeline, head, rate, length = 1.0, f_low = 100, f_high = 400, filter_latency = 0.5, td = True): length = int(length * rate / 2) @@ -295,7 +327,7 @@ def bandpass(pipeline, head, rate, length = 1.0, f_low = 100, f_high = 400, filt bandpass = numpy.convolve(highpass, lowpass) # Now apply the filter - return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * 2 * filter_latency), fir_matrix = [bandpass], time_domain = td) + return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * 2 * filter_latency + 0.5), fir_matrix = [bandpass], time_domain = td) def bandstop(pipeline, head, rate, length = 1.0, f_low = 100, f_high = 400, filter_latency = 0.5, td = True): length = int(length * rate / 2) @@ -324,7 +356,7 @@ def bandstop(pipeline, head, rate, length = 1.0, f_low = 100, f_high = 400, filt bandstop[length - 1] += 1 # Now apply the filter - return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * 2 * filter_latency), fir_matrix = [bandstop], time_domain = td) + return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * 2 * filter_latency + 0.5), fir_matrix = [bandstop], time_domain = td) def compute_rms(pipeline, head, rate, average_time, f_min = None, f_max = None, filter_latency = 0.5, rate_out = 16, td = True): # Find the root mean square amplitude of a signal between two frequencies