interp.patch

From 5cd668c98c9e1b48a04109a1c36356ec324a7e62 Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Wed, 11 Feb 2015 17:56:36 -0500
Subject: [PATCH 1/8] gstlal_interpolator: first commit

---
 gstlal/gst/lal/Makefile.am           |   1 +
 gstlal/gst/lal/gstlal.c              |   2 +
 gstlal/gst/lal/gstlal_interpolator.c | 631 +++++++++++++++++++++++++++++++++++
 gstlal/gst/lal/gstlal_interpolator.h | 127 +++++++
 4 files changed, 761 insertions(+)
 create mode 100644 gstlal/gst/lal/gstlal_interpolator.c
 create mode 100644 gstlal/gst/lal/gstlal_interpolator.h

diff --git a/gstlal/gst/lal/Makefile.am b/gstlal/gst/lal/Makefile.am
index 81d8ffb..ad98219 100644
--- a/gstlal/gst/lal/Makefile.am
+++ b/gstlal/gst/lal/Makefile.am
@@ -10,6 +10,7 @@ libgstlal_la_SOURCES = \
    gstlal_drop.h gstlal_drop.c \
    gstlal_firbank.h gstlal_firbank.c \
    gstlal_gate.h gstlal_gate.c \
+   gstlal_interpolator.h gstlal_interpolator.c \
    gstlal_matrixmixer.h gstlal_matrixmixer.c \
    gstlal_nofakedisconts.h gstlal_nofakedisconts.c \
    gstlal_nxydump.h gstlal_nxydump.c \
diff --git a/gstlal/gst/lal/gstlal.c b/gstlal/gst/lal/gstlal.c
index 8832ed7..d492df9 100644
--- a/gstlal/gst/lal/gstlal.c
+++ b/gstlal/gst/lal/gstlal.c
@@ -57,6 +57,7 @@
 #include <gstlal_drop.h>
 #include <gstlal_firbank.h>
 #include <gstlal_gate.h>
+#include <gstlal_interpolator.h>
 #include <gstlal_matrixmixer.h>
 #include <gstlal_nofakedisconts.h>
 #include <gstlal_nxydump.h>
@@ -92,6 +93,7 @@ static gboolean plugin_init(GstPlugin *plugin)
        {"lal_drop", GSTLAL_DROP_TYPE},
        {"lal_firbank", GSTLAL_FIRBANK_TYPE},
        {"lal_gate", GSTLAL_GATE_TYPE},
+       {"lal_interpolator", GSTLAL_INTERPOLATOR_TYPE},
        {"lal_matrixmixer", GSTLAL_MATRIXMIXER_TYPE},
        {"lal_nofakedisconts", GSTLAL_NOFAKEDISCONTS_TYPE},
        {"lal_nxydump", GST_TSVENC_TYPE},
diff --git a/gstlal/gst/lal/gstlal_interpolator.c b/gstlal/gst/lal/gstlal_interpolator.c
new file mode 100644
index 0000000..731d161
--- /dev/null
+++ b/gstlal/gst/lal/gstlal_interpolator.c
@@ -0,0 +1,631 @@
+/*
+ * An interpolator element
+ *
+ * Copyright (C) 2011  Chad Hanna, Kipp Cannon
+ *
+ * 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
+ *
+ * ========================================================================
+ */
+
+
+/*
+ * struff from the C library
+ */
+
+
+/*
+ * stuff from glib/gstreamer
+ */
+
+
+#include <glib.h>
+#include <glib/gprintf.h>
+#include <gst/gst.h>
+#include <gst/base/gstbasetransform.h>
+#include <gst/base/gstadapter.h>
+#include <math.h>
+#include <string.h>
+#include <fftw3.h>
+
+/*
+ * our own stuff
+ */
+
+#include <gstlal/gstlal.h>
+#include <gstlal/gstlal_debug.h>
+#include <gstlal/gstaudioadapter.h>
+#include <gstlal_interpolator.h>
+
+/*
+ * Utility functions for tapering the FFT buffers
+ */
+
+#define PI 3.141592653589793
+
+static float* taperup(int samps) {
+   float x;
+   float *out = (float *) calloc(samps, sizeof(float));
+   for (int i = 0; i < samps; i++) {
+       x = cos(PI / 2. * (float) i / samps);
+       out[i] = 1. - x * x;
+   }
+   return out;
+}
+
+static float* taperdown(int samps) {
+   float x;
+   float *out = (float *) calloc(samps, sizeof(float));
+   for (int i = 0; i < samps; i++) {
+       x = cos(PI / 2. * (float) i / samps);
+       out[i] = x * x;
+   }
+   return out;
+}
+
+static int applytaper(float *in, int end, float *taper) {
+   for (int i = 0; i < end; i++) {
+       in[i] *= taper[i];
+   }
+   return 0;
+}
+
+static int blend(float *in1, float *in2, int start, int end) {
+   for (int i = start; i < end; i++)
+       in1[i] += in2[i];
+   return 0;
+}
+
+
+/*
+ * gstreamer boiler plate
+ */
+
+#define GST_CAT_DEFAULT gstlal_interpolator_debug
+GST_DEBUG_CATEGORY_STATIC(GST_CAT_DEFAULT);
+
+static void additional_initializations(GType type)
+{
+   GST_DEBUG_CATEGORY_INIT(GST_CAT_DEFAULT, "lal_interpolator", 0, "lal_interpolator element");
+}
+
+GST_BOILERPLATE_FULL(
+        GSTLALInterpolator,
+        gstlal_interpolator,
+        GstBaseTransform,
+        GST_TYPE_BASE_TRANSFORM,
+        additional_initializations
+);
+
+static void gstlal_interpolator_base_init(gpointer klass){}
+
+/* Pads */
+
+static GstStaticPadTemplate sink_template =
+   GST_STATIC_PAD_TEMPLATE ("sink",
+   GST_PAD_SINK,
+   GST_PAD_ALWAYS,
+   GST_STATIC_CAPS ("audio/x-raw-float, "
+       "endianness = (int) BYTE_ORDER, "
+       "width = (int) 32, "
+       "rate = (int) {4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}, "
+       "channels = (int) [1, MAX]")
+   );
+
+static GstStaticPadTemplate src_template =
+   GST_STATIC_PAD_TEMPLATE ("src",
+   GST_PAD_SRC,
+   GST_PAD_ALWAYS,
+   GST_STATIC_CAPS ("audio/x-raw-float, "
+       "endianness = (int) BYTE_ORDER, "
+       "width = (int) 32, "
+       "rate = (int) {4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}, "
+       "channels = (int) [1, MAX]")
+   );
+
+/*
+ * Virtual method protototypes
+ */
+
+static void finalize(GObject *object);
+static gboolean get_unit_size(GstBaseTransform *trans, GstCaps *caps, guint *size);
+static gboolean set_caps (GstBaseTransform * base, GstCaps * incaps, GstCaps * outcaps);
+static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuffer *outbuf);
+static GstCaps* transform_caps (GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps);
+static gboolean transform_size(GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps, guint size, GstCaps *othercaps, guint *othersize);
+static gboolean start(GstBaseTransform *trans);
+static gboolean stop(GstBaseTransform *trans);
+
+/*
+ * class_init()
+ */
+
+static void gstlal_interpolator_class_init(GSTLALInterpolatorClass *klass)
+{
+
+        GObjectClass *gobject_class = G_OBJECT_CLASS(klass);
+        GstElementClass *element_class = GST_ELEMENT_CLASS(klass);
+        GstBaseTransformClass *transform_class = GST_BASE_TRANSFORM_CLASS(klass);
+
+        gst_element_class_set_details_simple(element_class, "Interpolator", "Filter/Audio", "Interpolates multichannel audio data using FFTs", "Chad Hanna <chad.hanna@ligo.org>");
+
+   gobject_class->finalize = GST_DEBUG_FUNCPTR(finalize);
+   transform_class->get_unit_size = GST_DEBUG_FUNCPTR(get_unit_size);
+   transform_class->set_caps = GST_DEBUG_FUNCPTR(set_caps);
+   transform_class->transform = GST_DEBUG_FUNCPTR(transform);
+   transform_class->transform_caps = GST_DEBUG_FUNCPTR(transform_caps);
+   transform_class->transform_size = GST_DEBUG_FUNCPTR(transform_size);
+   transform_class->start = GST_DEBUG_FUNCPTR(start);
+   transform_class->stop = GST_DEBUG_FUNCPTR(stop);
+
+        gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&src_template));
+        gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&sink_template));
+}
+
+
+static void gstlal_interpolator_init(GSTLALInterpolator *element, GSTLALInterpolatorClass *klass)
+{
+   gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(element), TRUE);
+
+   /* internal data */
+   element->inrate = 0;
+   element->outrate = 0;
+
+   element->nrin = 0; // size of real input to FFT
+   element->ncin = 0; // size of complex input to FFFT
+   element->nrout = 0; // size of real output to FFT
+   element->ncout = 0; // size of complex output to FFT
+   element->tapersampsin = 0;
+   element->tapersampsout = 0;
+   element->up = NULL;
+   element->down = NULL;
+   element->last = NULL;
+   element->rin = NULL;
+   element->cin = NULL;
+   element->rout = NULL;
+   element->cout = NULL;
+   element->data = NULL;
+
+   element->adapter = g_object_new(GST_TYPE_AUDIOADAPTER, NULL);
+}
+
+static GstCaps* transform_caps (GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps) {
+
+   /* 
+         * FIXME actually pull out the allowed rates so that we can prevent
+    *  downsampling at the negotiation stage
+    */
+   GstStructure *capsstruct;
+   gint channels;
+   capsstruct = gst_caps_get_structure (caps, 0);
+   char capsstr[256] = {0};
+
+   if (direction == GST_PAD_SINK && gst_structure_get_int (capsstruct, "channels", &channels)) {
+       sprintf(capsstr, "audio/x-raw-float, endianness = (int) BYTE_ORDER, width = (int) 32, rate = (int) {4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}, channels = (int) %d", channels);
+       return gst_caps_from_string(capsstr);
+   }
+
+   return gst_caps_from_string("audio/x-raw-float, endianness = (int) BYTE_ORDER, width = (int) 32, rate = (int) {4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}, channels = (int) [1, MAX]");
+
+}
+
+static gboolean set_caps (GstBaseTransform * base, GstCaps * incaps, GstCaps * outcaps) {
+   GSTLALInterpolator *element = GSTLAL_INTERPOLATOR (base);
+   GstStructure *instruct, *outstruct;
+   gint inchannels, inrate, outchannels, outrate;
+
+   instruct = gst_caps_get_structure (incaps, 0);
+   outstruct = gst_caps_get_structure (outcaps, 0);
+   g_return_val_if_fail(gst_structure_get_int (instruct, "channels", &inchannels), FALSE);
+   g_return_val_if_fail(gst_structure_get_int (instruct, "rate", &inrate), FALSE);
+   g_return_val_if_fail(gst_structure_get_int (outstruct, "channels", &outchannels), FALSE);
+   g_return_val_if_fail(gst_structure_get_int (outstruct, "rate", &outrate), FALSE);
+
+   g_return_val_if_fail(inchannels == outchannels, FALSE);
+   g_return_val_if_fail(outrate >= inrate, FALSE);
+   g_return_val_if_fail(outrate % inrate == 0, FALSE);
+
+   element->inrate = inrate;
+   element->outrate = outrate;
+   element->channels = inchannels;
+
+   /* Timestamp and offset bookeeping */
+
+   element->t0 = GST_CLOCK_TIME_NONE;
+        element->offset0 = GST_BUFFER_OFFSET_NONE;
+        element->next_input_offset = GST_BUFFER_OFFSET_NONE;
+        element->next_output_offset = GST_BUFFER_OFFSET_NONE;
+   element->need_discont = TRUE;
+
+   /* 
+    * NOTE: This element hardcodes transform sizes to be 1 second so we
+    * can set up a lot of stuff right now. Might as well.
+    */
+
+   element->nrin = inrate;
+   element->ncin = element->nrin / 2 + 1; // FFTW documentation for complex size
+   element->nrout = outrate;
+   element->ncout = element->nrout / 2 + 1; // FFTW documentation for complex size
+
+   element->tapersampsin = element->nrin / 4;
+   element->tapersampsout = element->nrout / 4;
+   element->blocksampsin = element->nrin - element->tapersampsin;
+   element->blocksampsout = element->nrout - element->tapersampsout;
+
+   if (element->up)
+       free(element->up);
+   element->up = taperup(element->tapersampsin);
+
+   if (element->down)
+       free(element->down);
+   element->down = taperdown(element->tapersampsin);
+
+   if (element->last)
+       free(element->last);
+   element->last = (float *) fftw_alloc_real(element->tapersampsout * element->channels);
+   memset(element->last, 0, sizeof(float) * element->tapersampsout * element->channels);
+       
+   if (!element->data)
+       free(element->data);
+   element->data = (float*) fftwf_alloc_real(element->nrin * element->channels);
+
+   if (element->rin)
+       fftwf_free(element->rin);
+   element->rin = (float*) fftwf_alloc_real(element->nrin);
+
+   if (element->cin)
+       fftwf_free(element->cin);
+   element->cin = (fftwf_complex*) fftwf_alloc_complex(element->ncin);
+
+   if (element->rout)
+       fftwf_free(element->rout);
+   element->rout = (float*) fftwf_alloc_real(element->nrout);
+
+   if (element->cout)
+       fftwf_free(element->cout);
+   element->cout = (fftwf_complex*) fftwf_alloc_complex(element->ncout);
+
+   gstlal_fftw_lock();
+   element->fwdplan_in = fftwf_plan_dft_r2c_1d(element->nrin, element->rin, element->cin, FFTW_PATIENT);
+   element->revplan_out = fftwf_plan_dft_c2r_1d(element->nrout, element->cout, element->rout, FFTW_PATIENT);
+   gstlal_fftw_unlock();
+
+   return TRUE;
+}
+
+
+/*
+ * ============================================================================
+ *
+ *                     GstBaseTransform Method Overrides
+ *
+ * ============================================================================
+ */
+
+
+/*
+ * get_unit_size()
+ */
+
+
+static gboolean get_unit_size(GstBaseTransform *trans, GstCaps *caps, guint *size)
+{
+   GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
+   GstStructure *str;
+   gint width, channels;
+   gboolean success = TRUE;
+
+   str = gst_caps_get_structure(caps, 0);
+   success &= gst_structure_get_int(str, "channels", &channels);
+   success &= gst_structure_get_int(str, "width", &width);
+       
+
+   if(success) {
+       *size = width / 8 * channels;
+       element->unitsize = *size;
+       g_object_set(element->adapter, "unit-size", *size, NULL);
+       GST_INFO_OBJECT(element, "get_unit_size(): channels %d, width %d", channels, width);
+   }
+   else
+       GST_WARNING_OBJECT(trans, "unable to parse channels from %" GST_PTR_FORMAT, caps);
+
+   return success;
+}
+
+static guint minimum_input_length(GSTLALInterpolator *element, guint samps) {
+   return ceil(samps / element->blocksampsin) * element->blocksampsin + element->tapersampsin;
+}
+
+static guint minimum_input_size(GSTLALInterpolator *element, guint size) {
+   return minimum_input_length(element, size / element->unitsize) * element->unitsize;
+}
+
+static guint64 get_available_samples(GSTLALInterpolator *element)
+{
+   //FIXME is size here really samples, I guess so??
+   guint size;
+   g_object_get(element->adapter, "size", &size, NULL);
+   return size;
+}
+
+static guint get_output_length(GSTLALInterpolator *element, guint samps) {
+   guint remainder;
+   guint numinsamps = get_available_samples(element) + samps;
+   if (numinsamps == 0)
+       return 0;
+   guint numoutsamps = numinsamps * element->outrate / element->inrate;
+   guint numblocks = numoutsamps / element->blocksampsout; //truncation
+   if (numblocks != 0)
+       remainder = numoutsamps % numblocks; // what doesn't fit into a multiple of block sizes
+   else
+       remainder = numoutsamps;
+   if ((remainder < element->tapersampsout) && (numblocks > 0)) // we can't actually produce output for the last tapersampsout worth of samples since those will be blended
+       numblocks -= 1;
+   return numblocks * element->blocksampsout; // Could be zero
+}
+
+static guint get_output_size(GSTLALInterpolator *element, guint size) {
+   return get_output_length(element, size / element->unitsize) * element->unitsize;
+}
+
+static gboolean transform_size(GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps, guint size, GstCaps *othercaps, guint *othersize)
+{
+   GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
+   guint unit_size;
+   guint other_unit_size;
+   gboolean success = TRUE;
+
+   if(!get_unit_size(trans, caps, &unit_size))
+       return FALSE;
+   if(size % unit_size) {
+       GST_ERROR_OBJECT(element, "size not a multiple of %u", unit_size);
+       return FALSE;
+   }
+   if(!get_unit_size(trans, othercaps, &other_unit_size))
+       return FALSE;
+
+   switch(direction) {
+   case GST_PAD_SRC:
+       /*
+        * number of input bytes required to produce an output
+        * buffer of (at least) the requested size
+        */
+       *othersize = minimum_input_size(element, size);
+       GST_INFO_OBJECT(element, "transform_size() producing %d (bytes) buffer for request on src pad", *othersize);
+       break;
+
+   case GST_PAD_SINK:
+       /*
+        * number of output bytes to be generated by the receipt of
+        * an input buffer of the given size.
+        */
+
+       *othersize = get_output_size(element, size);
+       GST_INFO_OBJECT(element, "transform_size() SINK pad buffer of size %d (bytes) %d (samples) provided. Transforming to size %d (bytes) %d (samples).", size, size / element->unitsize,  *othersize, *othersize / element->unitsize);
+       break;
+
+   case GST_PAD_UNKNOWN:
+       GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid direction GST_PAD_UNKNOWN"));
+       success = FALSE;
+       break;
+   }
+
+   return success;
+}
+
+static gboolean start(GstBaseTransform *trans)
+{
+   GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
+   element->t0 = GST_CLOCK_TIME_NONE;
+   element->offset0 = GST_BUFFER_OFFSET_NONE;
+   element->next_input_offset = GST_BUFFER_OFFSET_NONE;
+   element->next_output_offset = GST_BUFFER_OFFSET_NONE;
+   element->need_discont = TRUE;
+   //element->need_new_segment = TRUE;
+   return TRUE;
+}
+
+static gboolean stop(GstBaseTransform *trans)
+{
+        GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
+        g_object_unref(element->adapter);
+        element->adapter = NULL;
+        return TRUE;
+}
+
+static void flush_history(GSTLALInterpolator *element) {
+   GST_INFO_OBJECT(element, "flushing adapter contents");
+   gst_audioadapter_clear(element->adapter);
+}
+
+static void set_metadata(GSTLALInterpolator *element, GstBuffer *buf, guint64 outsamples, gboolean gap)
+{
+   GST_BUFFER_OFFSET(buf) = element->next_output_offset;
+   element->next_output_offset += outsamples;
+   GST_BUFFER_OFFSET_END(buf) = element->next_output_offset;
+   GST_BUFFER_TIMESTAMP(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET(buf) - element->offset0, GST_SECOND, element->outrate);
+   GST_BUFFER_DURATION(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, element->outrate) - 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);
+   GST_INFO_OBJECT(element, "set_metadata() %s%s output buffer %p spans %" GST_BUFFER_BOUNDARIES_FORMAT, gap ? "gap" : "nongap", GST_BUFFER_FLAG_IS_SET(buf, GST_BUFFER_FLAG_DISCONT) ? "+discont" : "", buf, GST_BUFFER_BOUNDARIES_ARGS(buf));
+}
+
+
+static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuffer *outbuf)
+{
+// FIXME, finish this. It is just (partially) copied from firbank
+   GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
+   guint expected_output_length, expected_output_size;
+   GstFlowReturn result = GST_FLOW_OK;
+
+   g_assert(GST_BUFFER_TIMESTAMP_IS_VALID(inbuf));
+   g_assert(GST_BUFFER_DURATION_IS_VALID(inbuf));
+   g_assert(GST_BUFFER_OFFSET_IS_VALID(inbuf));
+   g_assert(GST_BUFFER_OFFSET_END_IS_VALID(inbuf));
+
+   if(G_UNLIKELY(GST_BUFFER_IS_DISCONT(inbuf) || GST_BUFFER_OFFSET(inbuf) != element->next_input_offset || !GST_CLOCK_TIME_IS_VALID(element->t0))) {
+       /*
+        * flush any previous history and clear the adapter
+        */
+
+       flush_history(element);
+
+       /*
+        * (re)sync timestamp and offset book-keeping
+        */
+
+       element->t0 = GST_BUFFER_TIMESTAMP(inbuf);
+       element->offset0 = GST_BUFFER_OFFSET(inbuf);
+       element->next_output_offset = element->offset0 + get_output_length(element, GST_BUFFER_SIZE(inbuf) / element->unitsize);
+
+       /*
+        * be sure to flag the next output buffer as a discontinuity
+        */
+
+       element->need_discont = TRUE;
+       GST_INFO_OBJECT(element, "transform() A discontinuity was detected. The offset has been reset to %" G_GUINT64_FORMAT " and the timestamp has been reset to %" GST_TIME_SECONDS_FORMAT, element->offset0, element->t0);
+
+   } else
+       g_assert_cmpuint(GST_BUFFER_TIMESTAMP(inbuf), ==, gst_audioadapter_expected_timestamp(element->adapter));
+
+   element->next_input_offset = GST_BUFFER_OFFSET_END(inbuf);
+
+   /*
+    * Put the input buffer into the adapter first
+    * Then check the output buffer size that was expected.
+    * Note that the transform size function tells you what you can produce
+    * *after* receiving the next input buffer so this order is important
+    */
+
+   GST_INFO_OBJECT(element, "transform() pushing %d (bytes) %d (samples) sample buffer into adapter with size %d (samples)", GST_BUFFER_SIZE(inbuf), GST_BUFFER_SIZE(inbuf) / element->unitsize, get_available_samples(element));
+
+   gst_buffer_ref(inbuf);  /* don't let calling code free buffer */
+   gst_audioadapter_push(element->adapter, inbuf);
+   GST_INFO_OBJECT(element, "transform() adapter size is now %d (samples)", get_available_samples(element));
+
+   expected_output_length = get_output_length(element, 0); // just check the output length based on what we have in the adapter already
+   expected_output_size = get_output_size(element, 0); // Ditto here
+   GST_INFO_OBJECT(element, "transform() expect an output buffer with size %d (%d samples): got one with size %d", expected_output_length, expected_output_size, GST_BUFFER_SIZE(outbuf));
+   g_assert_cmpuint(expected_output_size, ==, GST_BUFFER_SIZE(outbuf));
+
+
+   /*
+    * Handle the different possiblilities
+    */
+
+   if (GST_BUFFER_SIZE(outbuf) == 0)
+       set_metadata(element, outbuf, 0, FALSE);
+   else {
+       guint flushed = 0;
+       guint processed = 0;
+       float *last = NULL;
+       float *output = (float *) GST_BUFFER_DATA(outbuf);
+
+       // FIXME actually handle gaps properly
+       while ((get_available_samples(element) >= element->nrin) && (processed < expected_output_length)) {
+
+           /* First copy the data we need */
+           gst_audioadapter_copy_samples(element->adapter, element->data, element->nrin, NULL, NULL);
+
+           /* Resample each channel */
+           for (guint i = 0; i < element->channels; i++) {
+
+               last = element->last + i * element->tapersampsout; // get a pointer to the history for this channel
+
+               /* Adapt output for FFTW FIXME FIXME FIXME we
+                * need to avoid this copying around by making use of the advanced FFTW
+                * interface
+                */
+               for (guint j = 0; j < element->nrin; j++) {
+                   element->rin[j] = element->data[i+j*element->channels];
+               }
+               /* Clear the output */
+               memset(element->cout, 0, sizeof(fftwf_complex) * element->ncout);
+
+               /* taper */
+               applytaper(element->rin, element->tapersampsin, element->up);
+               applytaper(element->rin + element->nrin - element->tapersampsin, element->tapersampsin, element->down);
+
+               /* resample */
+               fftwf_execute(element->fwdplan_in);
+               memcpy(element->cout, element->cin, sizeof(fftwf_complex) * element->ncin);
+               fftwf_execute(element->revplan_out);
+
+               /* Blend the outputs */
+               blend(element->rout, last, 0, element->tapersampsout);
+               memcpy(last, element->rout + (element->nrout - element->tapersampsout), sizeof(*last) * element->tapersampsout);
+
+               /* Adapt output for FFTW FIXME FIXME FIXME we
+                * need to avoid this copying around by making use of the advanced FFTW
+                * interface
+                */
+               for (guint j = 0; j < element->nrout-element->tapersampsout; j++) {
+                   output[i+j*element->channels] = element->rout[j] / element->nrin;
+               }
+           }
+
+           /* Then flush the data we will never need again */
+           // FIXME add a check that we have processed the correct number of samples
+           gst_audioadapter_flush_samples(element->adapter, element->nrin - element->tapersampsin);
+           flushed += element->nrin - element->tapersampsin;
+           processed += element->nrout - element->tapersampsout;
+           // Dont forget to increase the output pointer
+           output = (float *) GST_BUFFER_DATA(outbuf) + processed * element->channels;
+       }
+       GST_INFO_OBJECT(element, "flushed %d processed %d expected output length %d", flushed, processed, expected_output_length);
+       set_metadata(element, outbuf, expected_output_length, FALSE);
+
+   }
+   return result;
+}
+
+
+static void finalize(GObject *object)
+{
+        GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(object);
+
+        /*
+         * free resources
+         */
+   
+   free(element->up);
+   free(element->down);
+   free(element->last);
+   fftwf_free(element->rin);
+   fftwf_free(element->cin);
+   fftwf_free(element->rout);
+   fftwf_free(element->cout);
+   fftwf_free(element->data);
+
+        /*
+         * chain to parent class' finalize() method
+         */
+
+        G_OBJECT_CLASS(parent_class)->finalize(object);
+}
+
diff --git a/gstlal/gst/lal/gstlal_interpolator.h b/gstlal/gst/lal/gstlal_interpolator.h
new file mode 100644
index 0000000..e8be6cc
--- /dev/null
+++ b/gstlal/gst/lal/gstlal_interpolator.h
@@ -0,0 +1,127 @@
+/*
+ * Copyright (C) 2011 Chad Hanna <chad.hanna@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 __GSTLAL_INTERPOLATOR_H__
+#define __GSTLAL_INTERPOLATOR_H__
+
+
+#include <glib.h>
+#include <gst/gst.h>
+#include <gst/base/gstadapter.h>
+#include <gst/base/gstbasetransform.h>
+#include <fftw3.h>
+
+#include <gstlal/gstlal.h>
+#include <gstlal/gstaudioadapter.h>
+
+G_BEGIN_DECLS
+
+
+#define GSTLAL_INTERPOLATOR_TYPE \
+   (gstlal_interpolator_get_type())
+
+#define GSTLAL_INTERPOLATOR(obj) \
+   (G_TYPE_CHECK_INSTANCE_CAST((obj), GSTLAL_INTERPOLATOR_TYPE, GSTLALInterpolator))
+
+#define GSTLAL_INTERPOLATOR_CLASS(klass) \
+   (G_TYPE_CHECK_CLASS_CAST((klass), GSTLAL_INTERPOLATOR_TYPE, GSTLALInterpolatorClass))
+
+#define GST_IS_GSTLAL_INTERPOLATOR(obj) \
+   (G_TYPE_CHECK_INSTANCE_TYPE((obj), GSTLAL_INTERPOLATOR_TYPE))
+
+#define GST_IS_GSTLAL_INTERPOLATOR_CLASS(klass) \
+   (G_TYPE_CHECK_CLASS_TYPE((klass), GSTLAL_INTERPOLATOR_TYPE))
+
+
+typedef struct _GSTLALInterpolator GSTLALInterpolator;
+typedef struct _GSTLALInterpolatorClass GSTLALInterpolatorClass;
+
+
+
+
+
+struct _GSTLALInterpolator {
+   GstBaseTransform element;
+
+   gint inrate;
+   gint outrate;
+   guint channels;
+   GstAudioAdapter *adapter;
+   
+   /* Timestamp and offset bookeeping */
+   guint64 t0;
+   guint64 offset0;
+   guint64 next_input_offset;
+   GstClockTime next_input_timestamp;
+   guint64 next_output_offset;
+   GstClockTime next_output_timestamp;
+   gboolean need_discont;
+
+   /* Variables to control the size of transforms */
+   guint nrin;
+   guint ncin;
+   guint nrout;
+   guint ncout;
+   guint tapersampsin;
+   guint tapersampsout;
+   guint blocksampsin;
+   guint blocksampsout;
+   guint unitsize;
+
+   float *data;
+   float *up;
+   float *down;
+   float *last;
+   float *rin;
+   fftwf_complex *cin;
+   float *rout;
+   fftwf_complex *cout;
+   fftwf_plan fwdplan_in;
+   fftwf_plan revplan_out;
+};
+
+
+
+
+
+struct _GSTLALInterpolatorClass {
+        GstBaseTransformClass parent_class;
+};
+
+
+GType gstlal_interpolator_get_type(void);
+
+
+G_END_DECLS
+
+
+#endif /* __GSTLAL_INTERPOLATOR_H__ */
-- 
1.8.3.2


From 044e38137c40e312f798f60294acfebc422f7f8d Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Thu, 26 Feb 2015 16:23:20 -0500
Subject: [PATCH 2/8] update interpolator to time domain

---
 gstlal/gst/lal/gstlal_interpolator.c | 335 +++++++++++++++--------------------
 1 file changed, 138 insertions(+), 197 deletions(-)

diff --git a/gstlal/gst/lal/gstlal_interpolator.c b/gstlal/gst/lal/gstlal_interpolator.c
index 731d161..b4d9980 100644
--- a/gstlal/gst/lal/gstlal_interpolator.c
+++ b/gstlal/gst/lal/gstlal_interpolator.c
@@ -1,7 +1,7 @@
 /*
  * An interpolator element
  *
- * Copyright (C) 2011  Chad Hanna, Kipp Cannon
+ * Copyright (C) 2015  Chad Hanna, Kipp Cannon
  *
  * 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
@@ -62,40 +62,50 @@
 
 #define PI 3.141592653589793
 
-static float* taperup(int samps) {
-   float x;
-   float *out = (float *) calloc(samps, sizeof(float));
-   for (int i = 0; i < samps; i++) {
-       x = cos(PI / 2. * (float) i / samps);
-       out[i] = 1. - x * x;
+static float* kernel(int half_length, int factor) {
+   int kernel_length = (2 * half_length + 1) * factor;
+   int domain = kernel_length / 2; // kernel length is gauranteed to be even
+   float *out = fftwf_alloc_real(kernel_length + 1);
+   float norm = 0.;
+
+   for (int j = 0; j < kernel_length + 1; j++) {
+       float x = j - domain;
+       if (j == domain)
+           out[j] = 1.;
+       else
+           out[j] = sin(PI * x / factor) / (PI * x / factor) * (1. - x*x / domain / domain);
    }
-   return out;
-}
 
-static float* taperdown(int samps) {
-   float x;
-   float *out = (float *) calloc(samps, sizeof(float));
-   for (int i = 0; i < samps; i++) {
-       x = cos(PI / 2. * (float) i / samps);
-       out[i] = x * x;
-   }
+   for (int i = 0; i < kernel_length + 1; i++)
+       norm += out[i] * out[i];
+
+   for (int i = 0; i < kernel_length + 1; i++)
+       out[i] /= sqrt(norm / factor);
+
    return out;
 }
 
-static int applytaper(float *in, int end, float *taper) {
-   for (int i = 0; i < end; i++) {
-       in[i] *= taper[i];
+void convolve(float *output, float *thiskernel, float *input, guint kernel_length, guint factor, guint channels) {
+   for (guint i = 1; i < kernel_length; i++) {
+       *output += (*thiskernel) * (*input);
+       input += channels;
+       thiskernel += factor;
    }
-   return 0;
+   return;
 }
 
-static int blend(float *in1, float *in2, int start, int end) {
-   for (int i = start; i < end; i++)
-       in1[i] += in2[i];
-   return 0;
+void resample(float *output, float *thiskernel, float *input, guint kernel_length, guint factor, guint channels, guint blockstrideout) {
+   guint kernel_offset, output_offset, input_offset;
+   for (gint samp = 0; samp < blockstrideout; samp++) {
+       kernel_offset = factor - samp % factor + factor / 2;
+       output_offset = samp * channels;
+       input_offset = (samp / factor) * channels; // first input sample
+       for (gint chan = 0; chan < channels; chan++)
+           convolve(output + output_offset + chan, thiskernel + kernel_offset, input + input_offset + chan + channels, kernel_length, factor, channels);
+   }
+   return;
 }
 
-
 /*
  * gstreamer boiler plate
  */
@@ -166,7 +176,7 @@ static void gstlal_interpolator_class_init(GSTLALInterpolatorClass *klass)
         GstElementClass *element_class = GST_ELEMENT_CLASS(klass);
         GstBaseTransformClass *transform_class = GST_BASE_TRANSFORM_CLASS(klass);
 
-        gst_element_class_set_details_simple(element_class, "Interpolator", "Filter/Audio", "Interpolates multichannel audio data using FFTs", "Chad Hanna <chad.hanna@ligo.org>");
+        gst_element_class_set_details_simple(element_class, "Interpolator", "Filter/Audio", "Interpolates multichannel audio data using FFTs", "Chad Hanna <chad.hanna@ligo.org>, Patrick Brockill <brockill@uwm.edu>");
 
    gobject_class->finalize = GST_DEBUG_FUNCPTR(finalize);
    transform_class->get_unit_size = GST_DEBUG_FUNCPTR(get_unit_size);
@@ -190,20 +200,17 @@ static void gstlal_interpolator_init(GSTLALInterpolator *element, GSTLALInterpol
    element->inrate = 0;
    element->outrate = 0;
 
-   element->nrin = 0; // size of real input to FFT
-   element->ncin = 0; // size of complex input to FFFT
-   element->nrout = 0; // size of real output to FFT
-   element->ncout = 0; // size of complex output to FFT
-   element->tapersampsin = 0;
-   element->tapersampsout = 0;
-   element->up = NULL;
-   element->down = NULL;
-   element->last = NULL;
-   element->rin = NULL;
-   element->cin = NULL;
-   element->rout = NULL;
-   element->cout = NULL;
-   element->data = NULL;
+   element->factor = 0; // size of complex output to FFT
+   element->kernel = NULL;
+   element->workspace = NULL;
+
+   // hardcoded kernel size
+   element->half_length = 8;
+   element->kernel_length = element->half_length * 2 + 1;
+
+   // Always initialize with a discont
+   element->need_discont = TRUE;
+   element->need_pretend = TRUE;
 
    element->adapter = g_object_new(GST_TYPE_AUDIOADAPTER, NULL);
 }
@@ -247,6 +254,7 @@ static gboolean set_caps (GstBaseTransform * base, GstCaps * incaps, GstCaps * o
    element->inrate = inrate;
    element->outrate = outrate;
    element->channels = inchannels;
+   element->factor = outrate / inrate;
 
    /* Timestamp and offset bookeeping */
 
@@ -255,59 +263,24 @@ static gboolean set_caps (GstBaseTransform * base, GstCaps * incaps, GstCaps * o
         element->next_input_offset = GST_BUFFER_OFFSET_NONE;
         element->next_output_offset = GST_BUFFER_OFFSET_NONE;
    element->need_discont = TRUE;
+   element->need_pretend = TRUE;
 
-   /* 
-    * NOTE: This element hardcodes transform sizes to be 1 second so we
-    * can set up a lot of stuff right now. Might as well.
-    */
+   if (element->kernel)
+       free(element->kernel);
+   element->kernel = kernel(element->half_length, element->factor);
 
-   element->nrin = inrate;
-   element->ncin = element->nrin / 2 + 1; // FFTW documentation for complex size
-   element->nrout = outrate;
-   element->ncout = element->nrout / 2 + 1; // FFTW documentation for complex size
+   // Assume that we process inrate worth of samples at a time (e.g. 1s stride)
+   element->blockstridein = 128;//element->inrate;
+   element->blocksampsin = element->blockstridein + element->kernel_length;
+   element->blockstrideout = element->blockstridein * element->factor;//element->outrate;
+   element->blocksampsout = element->blockstrideout + (element->kernel_length) * element->factor;
 
-   element->tapersampsin = element->nrin / 4;
-   element->tapersampsout = element->nrout / 4;
-   element->blocksampsin = element->nrin - element->tapersampsin;
-   element->blocksampsout = element->nrout - element->tapersampsout;
+   GST_INFO_OBJECT(element, "blocksampsin %d, blocksampsout %d, blockstridein %d, blockstrideout %d", element->blocksampsin, element->blocksampsout, element->blockstridein, element->blockstrideout);
 
-   if (element->up)
-       free(element->up);
-   element->up = taperup(element->tapersampsin);
-
-   if (element->down)
-       free(element->down);
-   element->down = taperdown(element->tapersampsin);
-
-   if (element->last)
-       free(element->last);
-   element->last = (float *) fftw_alloc_real(element->tapersampsout * element->channels);
-   memset(element->last, 0, sizeof(float) * element->tapersampsout * element->channels);
-       
-   if (!element->data)
-       free(element->data);
-   element->data = (float*) fftwf_alloc_real(element->nrin * element->channels);
-
-   if (element->rin)
-       fftwf_free(element->rin);
-   element->rin = (float*) fftwf_alloc_real(element->nrin);
-
-   if (element->cin)
-       fftwf_free(element->cin);
-   element->cin = (fftwf_complex*) fftwf_alloc_complex(element->ncin);
-
-   if (element->rout)
-       fftwf_free(element->rout);
-   element->rout = (float*) fftwf_alloc_real(element->nrout);
-
-   if (element->cout)
-       fftwf_free(element->cout);
-   element->cout = (fftwf_complex*) fftwf_alloc_complex(element->ncout);
-
-   gstlal_fftw_lock();
-   element->fwdplan_in = fftwf_plan_dft_r2c_1d(element->nrin, element->rin, element->cin, FFTW_PATIENT);
-   element->revplan_out = fftwf_plan_dft_c2r_1d(element->nrout, element->cout, element->rout, FFTW_PATIENT);
-   gstlal_fftw_unlock();
+   if (element->workspace)
+       free(element->workspace);
+   element->workspace = (float *) fftw_alloc_real(element->blocksampsin * element->channels);
+   memset(element->workspace, 0, sizeof(float) * element->blocksampsin * element->channels);
 
    return TRUE;
 }
@@ -343,7 +316,7 @@ static gboolean get_unit_size(GstBaseTransform *trans, GstCaps *caps, guint *siz
        *size = width / 8 * channels;
        element->unitsize = *size;
        g_object_set(element->adapter, "unit-size", *size, NULL);
-       GST_INFO_OBJECT(element, "get_unit_size(): channels %d, width %d", channels, width);
+       GST_INFO_OBJECT(element, "channels %d, width %d", channels, width);
    }
    else
        GST_WARNING_OBJECT(trans, "unable to parse channels from %" GST_PTR_FORMAT, caps);
@@ -351,14 +324,6 @@ static gboolean get_unit_size(GstBaseTransform *trans, GstCaps *caps, guint *siz
    return success;
 }
 
-static guint minimum_input_length(GSTLALInterpolator *element, guint samps) {
-   return ceil(samps / element->blocksampsin) * element->blocksampsin + element->tapersampsin;
-}
-
-static guint minimum_input_size(GSTLALInterpolator *element, guint size) {
-   return minimum_input_length(element, size / element->unitsize) * element->unitsize;
-}
-
 static guint64 get_available_samples(GSTLALInterpolator *element)
 {
    //FIXME is size here really samples, I guess so??
@@ -367,20 +332,25 @@ static guint64 get_available_samples(GSTLALInterpolator *element)
    return size;
 }
 
+
+static guint minimum_input_length(GSTLALInterpolator *element, guint samps) {
+   return samps / element->factor + element->kernel_length; // FIXME check this
+}
+
+static guint minimum_input_size(GSTLALInterpolator *element, guint size) {
+   return minimum_input_length(element, size / element->unitsize) * element->unitsize;
+}
+
 static guint get_output_length(GSTLALInterpolator *element, guint samps) {
-   guint remainder;
-   guint numinsamps = get_available_samples(element) + samps;
-   if (numinsamps == 0)
+   // Pretend that we have a half_length set of samples if we are at a discont
+   guint pretend_samps = element->need_pretend ? element->half_length : 0;
+   guint numinsamps = get_available_samples(element) + samps + pretend_samps;
+   if (numinsamps <= element->kernel_length)
        return 0;
-   guint numoutsamps = numinsamps * element->outrate / element->inrate;
-   guint numblocks = numoutsamps / element->blocksampsout; //truncation
-   if (numblocks != 0)
-       remainder = numoutsamps % numblocks; // what doesn't fit into a multiple of block sizes
-   else
-       remainder = numoutsamps;
-   if ((remainder < element->tapersampsout) && (numblocks > 0)) // we can't actually produce output for the last tapersampsout worth of samples since those will be blended
-       numblocks -= 1;
-   return numblocks * element->blocksampsout; // Could be zero
+   guint numoutsamps = (numinsamps - element->kernel_length) * element->factor; // FIXME check this
+   guint numblocks = numoutsamps / element->blockstrideout; //truncation
+
+   return numblocks * element->blockstrideout; // Could be zero
 }
 
 static guint get_output_size(GSTLALInterpolator *element, guint size) {
@@ -410,7 +380,7 @@ static gboolean transform_size(GstBaseTransform *trans, GstPadDirection directio
         * buffer of (at least) the requested size
         */
        *othersize = minimum_input_size(element, size);
-       GST_INFO_OBJECT(element, "transform_size() producing %d (bytes) buffer for request on src pad", *othersize);
+       GST_INFO_OBJECT(element, "producing %d (bytes) buffer for request on SRC pad", *othersize);
        break;
 
    case GST_PAD_SINK:
@@ -420,7 +390,7 @@ static gboolean transform_size(GstBaseTransform *trans, GstPadDirection directio
         */
 
        *othersize = get_output_size(element, size);
-       GST_INFO_OBJECT(element, "transform_size() SINK pad buffer of size %d (bytes) %d (samples) provided. Transforming to size %d (bytes) %d (samples).", size, size / element->unitsize,  *othersize, *othersize / element->unitsize);
+       GST_INFO_OBJECT(element, "SINK pad buffer of size %d (bytes) %d (samples) provided. Transforming to size %d (bytes) %d (samples).", size, size / element->unitsize,  *othersize, *othersize / element->unitsize);
        break;
 
    case GST_PAD_UNKNOWN:
@@ -440,15 +410,17 @@ static gboolean start(GstBaseTransform *trans)
    element->next_input_offset = GST_BUFFER_OFFSET_NONE;
    element->next_output_offset = GST_BUFFER_OFFSET_NONE;
    element->need_discont = TRUE;
-   //element->need_new_segment = TRUE;
+   element->need_pretend = TRUE;
+   // FIXME properly handle segments
+   // element->need_new_segment = TRUE;
    return TRUE;
 }
 
 static gboolean stop(GstBaseTransform *trans)
 {
-        GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
-        g_object_unref(element->adapter);
-        element->adapter = NULL;
+   GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
+   g_object_unref(element->adapter);
+   element->adapter = NULL;
         return TRUE;
 }
 
@@ -472,22 +444,21 @@ static void set_metadata(GSTLALInterpolator *element, GstBuffer *buf, guint64 ou
        GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_GAP);
    else
        GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP);
-   GST_INFO_OBJECT(element, "set_metadata() %s%s output buffer %p spans %" GST_BUFFER_BOUNDARIES_FORMAT, gap ? "gap" : "nongap", GST_BUFFER_FLAG_IS_SET(buf, GST_BUFFER_FLAG_DISCONT) ? "+discont" : "", buf, GST_BUFFER_BOUNDARIES_ARGS(buf));
+   GST_INFO_OBJECT(element, "%s%s output_buffer %p spans %" GST_BUFFER_BOUNDARIES_FORMAT, gap ? "gap" : "nongap", GST_BUFFER_FLAG_IS_SET(buf, GST_BUFFER_FLAG_DISCONT) ? "+discont" : "", buf, GST_BUFFER_BOUNDARIES_ARGS(buf));
 }
 
 
 static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuffer *outbuf)
 {
-// FIXME, finish this. It is just (partially) copied from firbank
    GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
-   guint expected_output_length, expected_output_size;
+   guint input_length, output_length, expected_output_size;
    GstFlowReturn result = GST_FLOW_OK;
 
    g_assert(GST_BUFFER_TIMESTAMP_IS_VALID(inbuf));
    g_assert(GST_BUFFER_DURATION_IS_VALID(inbuf));
    g_assert(GST_BUFFER_OFFSET_IS_VALID(inbuf));
    g_assert(GST_BUFFER_OFFSET_END_IS_VALID(inbuf));
-
+       
    if(G_UNLIKELY(GST_BUFFER_IS_DISCONT(inbuf) || GST_BUFFER_OFFSET(inbuf) != element->next_input_offset || !GST_CLOCK_TIME_IS_VALID(element->t0))) {
        /*
         * flush any previous history and clear the adapter
@@ -501,104 +472,80 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
 
        element->t0 = GST_BUFFER_TIMESTAMP(inbuf);
        element->offset0 = GST_BUFFER_OFFSET(inbuf);
-       element->next_output_offset = element->offset0 + get_output_length(element, GST_BUFFER_SIZE(inbuf) / element->unitsize);
+       element->next_output_offset = element->offset0;
 
        /*
         * be sure to flag the next output buffer as a discontinuity
         */
 
        element->need_discont = TRUE;
-       GST_INFO_OBJECT(element, "transform() A discontinuity was detected. The offset has been reset to %" G_GUINT64_FORMAT " and the timestamp has been reset to %" GST_TIME_SECONDS_FORMAT, element->offset0, element->t0);
+       element->need_pretend = TRUE;
+       GST_INFO_OBJECT(element, "A discontinuity was detected. The offset has been reset to %" G_GUINT64_FORMAT " and the timestamp has been reset to %" GST_TIME_SECONDS_FORMAT, element->offset0, element->t0);
 
-   } else
+   }
+   else {
        g_assert_cmpuint(GST_BUFFER_TIMESTAMP(inbuf), ==, gst_audioadapter_expected_timestamp(element->adapter));
+   }
+
 
    element->next_input_offset = GST_BUFFER_OFFSET_END(inbuf);
 
-   /*
-    * Put the input buffer into the adapter first
-    * Then check the output buffer size that was expected.
-    * Note that the transform size function tells you what you can produce
-    * *after* receiving the next input buffer so this order is important
-    */
+   GST_INFO_OBJECT(element, "%s input_buffer %p spans %" GST_BUFFER_BOUNDARIES_FORMAT, GST_BUFFER_FLAG_IS_SET(inbuf, GST_BUFFER_FLAG_DISCONT) ? "+discont" : "", inbuf, GST_BUFFER_BOUNDARIES_ARGS(inbuf));
 
-   GST_INFO_OBJECT(element, "transform() pushing %d (bytes) %d (samples) sample buffer into adapter with size %d (samples)", GST_BUFFER_SIZE(inbuf), GST_BUFFER_SIZE(inbuf) / element->unitsize, get_available_samples(element));
+   GST_INFO_OBJECT(element, "pushing %d (bytes) %d (samples) sample buffer into adapter with size %d (samples)", GST_BUFFER_SIZE(inbuf), GST_BUFFER_SIZE(inbuf) / element->unitsize, get_available_samples(element));
 
    gst_buffer_ref(inbuf);  /* don't let calling code free buffer */
    gst_audioadapter_push(element->adapter, inbuf);
-   GST_INFO_OBJECT(element, "transform() adapter size is now %d (samples)", get_available_samples(element));
-
-   expected_output_length = get_output_length(element, 0); // just check the output length based on what we have in the adapter already
-   expected_output_size = get_output_size(element, 0); // Ditto here
-   GST_INFO_OBJECT(element, "transform() expect an output buffer with size %d (%d samples): got one with size %d", expected_output_length, expected_output_size, GST_BUFFER_SIZE(outbuf));
-   g_assert_cmpuint(expected_output_size, ==, GST_BUFFER_SIZE(outbuf));
-
+   GST_INFO_OBJECT(element, "adapter_size %d (samples)", get_available_samples(element));
+   
+   // FIXME check the sanity of the output buffer
 
    /*
     * Handle the different possiblilities
     */
 
+   output_length = GST_BUFFER_SIZE(outbuf) / element->unitsize;
+
    if (GST_BUFFER_SIZE(outbuf) == 0)
        set_metadata(element, outbuf, 0, FALSE);
    else {
+
+
        guint flushed = 0;
        guint processed = 0;
-       float *last = NULL;
+       gint input_offset;
+       gint kernel_offset;
+       gint output_offset;
+       gint f;
+       gint i;
        float *output = (float *) GST_BUFFER_DATA(outbuf);
+       memset(GST_BUFFER_DATA(outbuf), 0, GST_BUFFER_SIZE(outbuf));  // FIXME necesary?
 
-       // FIXME actually handle gaps properly
-       while ((get_available_samples(element) >= element->nrin) && (processed < expected_output_length)) {
-
-           /* First copy the data we need */
-           gst_audioadapter_copy_samples(element->adapter, element->data, element->nrin, NULL, NULL);
-
-           /* Resample each channel */
-           for (guint i = 0; i < element->channels; i++) {
-
-               last = element->last + i * element->tapersampsout; // get a pointer to the history for this channel
-
-               /* Adapt output for FFTW FIXME FIXME FIXME we
-                * need to avoid this copying around by making use of the advanced FFTW
-                * interface
-                */
-               for (guint j = 0; j < element->nrin; j++) {
-                   element->rin[j] = element->data[i+j*element->channels];
-               }
-               /* Clear the output */
-               memset(element->cout, 0, sizeof(fftwf_complex) * element->ncout);
-
-               /* taper */
-               applytaper(element->rin, element->tapersampsin, element->up);
-               applytaper(element->rin + element->nrin - element->tapersampsin, element->tapersampsin, element->down);
-
-               /* resample */
-               fftwf_execute(element->fwdplan_in);
-               memcpy(element->cout, element->cin, sizeof(fftwf_complex) * element->ncin);
-               fftwf_execute(element->revplan_out);
-
-               /* Blend the outputs */
-               blend(element->rout, last, 0, element->tapersampsout);
-               memcpy(last, element->rout + (element->nrout - element->tapersampsout), sizeof(*last) * element->tapersampsout);
-
-               /* Adapt output for FFTW FIXME FIXME FIXME we
-                * need to avoid this copying around by making use of the advanced FFTW
-                * interface
-                */
-               for (guint j = 0; j < element->nrout-element->tapersampsout; j++) {
-                   output[i+j*element->channels] = element->rout[j] / element->nrin;
-               }
-           }
+       GST_INFO_OBJECT(element, "Processing a %d sample output buffer from %d input", output_length);
+
+       while (processed < output_length) {
+
+           memset(element->workspace, 0, sizeof(float) * element->blocksampsin * element->channels); // FIXME necessary?
 
-           /* Then flush the data we will never need again */
-           // FIXME add a check that we have processed the correct number of samples
-           gst_audioadapter_flush_samples(element->adapter, element->nrin - element->tapersampsin);
-           flushed += element->nrin - element->tapersampsin;
-           processed += element->nrout - element->tapersampsout;
-           // Dont forget to increase the output pointer
-           output = (float *) GST_BUFFER_DATA(outbuf) + processed * element->channels;
+           // Special case to handle discontinuities: effectively zero pad. FIXME make this more elegant
+           if (element->need_pretend)
+               gst_audioadapter_copy_samples(element->adapter, element->workspace + (element->half_length) * element->channels, element->blocksampsin - element->half_length, NULL, NULL);
+           else
+               gst_audioadapter_copy_samples(element->adapter, element->workspace, element->blocksampsin, NULL, NULL);
+
+           resample(output, element->kernel, element->workspace, element->kernel_length, element->factor, element->channels, element->blockstrideout);
+
+           if (element->need_pretend) {
+               element->need_pretend = FALSE;
+               gst_audioadapter_flush_samples(element->adapter, element->blockstridein - element->half_length);
+           }
+           else
+               gst_audioadapter_flush_samples(element->adapter, element->blockstridein);
+           output += element->blockstrideout * element->channels;
+           processed += element->blockstrideout;
        }
-       GST_INFO_OBJECT(element, "flushed %d processed %d expected output length %d", flushed, processed, expected_output_length);
-       set_metadata(element, outbuf, expected_output_length, FALSE);
+       GST_INFO_OBJECT(element, "Processed a %d samples", processed);
+       set_metadata(element, outbuf, output_length, FALSE);
 
    }
    return result;
@@ -613,14 +560,8 @@ static void finalize(GObject *object)
          * free resources
          */
    
-   free(element->up);
-   free(element->down);
-   free(element->last);
-   fftwf_free(element->rin);
-   fftwf_free(element->cin);
-   fftwf_free(element->rout);
-   fftwf_free(element->cout);
-   fftwf_free(element->data);
+   free(element->kernel);
+   free(element->workspace);
 
         /*
          * chain to parent class' finalize() method
-- 
1.8.3.2


From 74a62f546d2e895a908d7acba4ffc560cdd9134b Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Thu, 26 Feb 2015 16:23:22 -0500
Subject: [PATCH 3/8] update interpolator to time domain

---
 gstlal/gst/lal/gstlal_interpolator.h | 30 +++++++++++-------------------
 1 file changed, 11 insertions(+), 19 deletions(-)

diff --git a/gstlal/gst/lal/gstlal_interpolator.h b/gstlal/gst/lal/gstlal_interpolator.h
index e8be6cc..4860c79 100644
--- a/gstlal/gst/lal/gstlal_interpolator.h
+++ b/gstlal/gst/lal/gstlal_interpolator.h
@@ -69,10 +69,11 @@ typedef struct _GSTLALInterpolatorClass GSTLALInterpolatorClass;
 struct _GSTLALInterpolator {
    GstBaseTransform element;
 
+        GstAudioAdapter *adapter;
+
    gint inrate;
    gint outrate;
    guint channels;
-   GstAudioAdapter *adapter;
    
    /* Timestamp and offset bookeeping */
    guint64 t0;
@@ -82,28 +83,19 @@ struct _GSTLALInterpolator {
    guint64 next_output_offset;
    GstClockTime next_output_timestamp;
    gboolean need_discont;
+   gboolean need_pretend;
 
    /* Variables to control the size of transforms */
-   guint nrin;
-   guint ncin;
-   guint nrout;
-   guint ncout;
-   guint tapersampsin;
-   guint tapersampsout;
+   guint unitsize;
    guint blocksampsin;
    guint blocksampsout;
-   guint unitsize;
-
-   float *data;
-   float *up;
-   float *down;
-   float *last;
-   float *rin;
-   fftwf_complex *cin;
-   float *rout;
-   fftwf_complex *cout;
-   fftwf_plan fwdplan_in;
-   fftwf_plan revplan_out;
+   guint blockstridein;
+   guint blockstrideout;
+   guint factor;
+   guint half_length;
+   guint kernel_length;
+   float *kernel;
+   float *workspace;
 };
 
 
-- 
1.8.3.2


From 2a5224390fa5ebef5e680891eb979bbad882452e Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Fri, 27 Feb 2015 14:09:35 -0500
Subject: [PATCH 4/8] gstlal_interpolator: additional vectorization

---
 gstlal/gst/lal/gstlal_interpolator.c | 97 +++++++++++++++++++-----------------
 1 file changed, 52 insertions(+), 45 deletions(-)

diff --git a/gstlal/gst/lal/gstlal_interpolator.c b/gstlal/gst/lal/gstlal_interpolator.c
index b4d9980..97741ab 100644
--- a/gstlal/gst/lal/gstlal_interpolator.c
+++ b/gstlal/gst/lal/gstlal_interpolator.c
@@ -46,6 +46,8 @@
 #include <math.h>
 #include <string.h>
 #include <fftw3.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
 
 /*
  * our own stuff
@@ -62,46 +64,56 @@
 
 #define PI 3.141592653589793
 
-static float* kernel(int half_length, int factor) {
+gsl_vector_float** kernel(int half_length, int factor) {
    int kernel_length = (2 * half_length + 1) * factor;
    int domain = kernel_length / 2; // kernel length is gauranteed to be even
-   float *out = fftwf_alloc_real(kernel_length + 1);
+
+   gsl_vector_float **vecs = malloc(sizeof(gsl_vector_float *) * factor);
+   for (int i = 0; i < factor; i++)
+       vecs[i] = gsl_vector_float_calloc(2 * half_length + 1);
+
+   float *out = fftwf_malloc(sizeof(float) * (kernel_length + factor / 2));
+   memset(out, 0, (kernel_length + factor / 2) * sizeof(float));
+
    float norm = 0.;
 
-   for (int j = 0; j < kernel_length + 1; j++) {
-       float x = j - domain;
-       if (j == domain)
-           out[j] = 1.;
-       else
-           out[j] = sin(PI * x / factor) / (PI * x / factor) * (1. - x*x / domain / domain);
+   for (int j = 0; j < 2 * half_length + 1; j++) {
+       for (int i = 0; i < factor; i++) {
+           int x = j * factor + i - domain;
+           if (x == 0)
+               out[x + domain] = 1.;
+           else
+               out[x + domain] = sin(PI * x / factor) / (PI * x / factor) * (1. - (float) x*x / domain / domain);
+           norm += out[x + domain] * out[x + domain];
+       }
    }
 
-   for (int i = 0; i < kernel_length + 1; i++)
-       norm += out[i] * out[i];
-
-   for (int i = 0; i < kernel_length + 1; i++)
-       out[i] /= sqrt(norm / factor);
+   for (int i = 0; i < 2 * half_length+1; i++) {
+       for (int j = 0; j < factor; j++) {
+           int index = i * factor + j + factor / 2; //FIXME FIXME FIXME this offset belongs above
+           gsl_vector_float_set(vecs[factor - j - 1], i, out[index] / sqrt(norm / factor));
+       }
+   }
 
-   return out;
+   free(out);
+   return vecs;
 }
 
-void convolve(float *output, float *thiskernel, float *input, guint kernel_length, guint factor, guint channels) {
-   for (guint i = 1; i < kernel_length; i++) {
-       *output += (*thiskernel) * (*input);
-       input += channels;
-       thiskernel += factor;
-   }
+void convolve(float *output, gsl_vector_float *thiskernel, float *input, guint kernel_length, guint factor, guint channels) {
+   gsl_vector_float_view output_vector = gsl_vector_float_view_array(output, channels);
+   gsl_matrix_float_view input_matrix = gsl_matrix_float_view_array(input, kernel_length, channels);
+   gsl_blas_sgemv (CblasTrans, 1.0, &(input_matrix.matrix), thiskernel, 0, &(output_vector.vector));
    return;
 }
 
-void resample(float *output, float *thiskernel, float *input, guint kernel_length, guint factor, guint channels, guint blockstrideout) {
+void resample(float *output, gsl_vector_float **thiskernel, float *input, guint kernel_length, guint factor, guint channels, guint blockstrideout) {
    guint kernel_offset, output_offset, input_offset;
-   for (gint samp = 0; samp < blockstrideout; samp++) {
-       kernel_offset = factor - samp % factor + factor / 2;
+   for (guint samp = 0; samp < blockstrideout; samp++) {
+       kernel_offset = samp % factor;
        output_offset = samp * channels;
        input_offset = (samp / factor) * channels; // first input sample
-       for (gint chan = 0; chan < channels; chan++)
-           convolve(output + output_offset + chan, thiskernel + kernel_offset, input + input_offset + chan + channels, kernel_length, factor, channels);
+       // FIXME FIXME FIXME should this + channels really be there??
+       convolve(output + output_offset, thiskernel[kernel_offset], input + input_offset + channels, kernel_length, factor, channels);
    }
    return;
 }
@@ -270,7 +282,7 @@ static gboolean set_caps (GstBaseTransform * base, GstCaps * incaps, GstCaps * o
    element->kernel = kernel(element->half_length, element->factor);
 
    // Assume that we process inrate worth of samples at a time (e.g. 1s stride)
-   element->blockstridein = 128;//element->inrate;
+   element->blockstridein = element->inrate;
    element->blocksampsin = element->blockstridein + element->kernel_length;
    element->blockstrideout = element->blockstridein * element->factor;//element->outrate;
    element->blocksampsout = element->blockstrideout + (element->kernel_length) * element->factor;
@@ -278,9 +290,8 @@ static gboolean set_caps (GstBaseTransform * base, GstCaps * incaps, GstCaps * o
    GST_INFO_OBJECT(element, "blocksampsin %d, blocksampsout %d, blockstridein %d, blockstrideout %d", element->blocksampsin, element->blocksampsout, element->blockstridein, element->blockstrideout);
 
    if (element->workspace)
-       free(element->workspace);
-   element->workspace = (float *) fftw_alloc_real(element->blocksampsin * element->channels);
-   memset(element->workspace, 0, sizeof(float) * element->blocksampsin * element->channels);
+       gsl_matrix_float_free(element->workspace);
+   element->workspace = gsl_matrix_float_calloc (element->blocksampsin, element->channels);
 
    return TRUE;
 }
@@ -451,7 +462,7 @@ static void set_metadata(GSTLALInterpolator *element, GstBuffer *buf, guint64 ou
 static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuffer *outbuf)
 {
    GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
-   guint input_length, output_length, expected_output_size;
+   guint output_length;
    GstFlowReturn result = GST_FLOW_OK;
 
    g_assert(GST_BUFFER_TIMESTAMP_IS_VALID(inbuf));
@@ -511,29 +522,24 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
    else {
 
 
-       guint flushed = 0;
        guint processed = 0;
-       gint input_offset;
-       gint kernel_offset;
-       gint output_offset;
-       gint f;
-       gint i;
        float *output = (float *) GST_BUFFER_DATA(outbuf);
-       memset(GST_BUFFER_DATA(outbuf), 0, GST_BUFFER_SIZE(outbuf));  // FIXME necesary?
+       //memset(GST_BUFFER_DATA(outbuf), 0, GST_BUFFER_SIZE(outbuf));  // FIXME necesary?
 
        GST_INFO_OBJECT(element, "Processing a %d sample output buffer from %d input", output_length);
 
        while (processed < output_length) {
 
-           memset(element->workspace, 0, sizeof(float) * element->blocksampsin * element->channels); // FIXME necessary?
 
            // Special case to handle discontinuities: effectively zero pad. FIXME make this more elegant
-           if (element->need_pretend)
-               gst_audioadapter_copy_samples(element->adapter, element->workspace + (element->half_length) * element->channels, element->blocksampsin - element->half_length, NULL, NULL);
+           if (element->need_pretend) {
+               memset(element->workspace->data, 0, sizeof(float) * element->workspace->size1 * element->workspace->size2); // FIXME necessary?
+               gst_audioadapter_copy_samples(element->adapter, element->workspace->data + (element->half_length) * element->channels, element->blocksampsin - element->half_length, NULL, NULL);
+           }
            else
-               gst_audioadapter_copy_samples(element->adapter, element->workspace, element->blocksampsin, NULL, NULL);
+               gst_audioadapter_copy_samples(element->adapter, element->workspace->data, element->blocksampsin, NULL, NULL);
 
-           resample(output, element->kernel, element->workspace, element->kernel_length, element->factor, element->channels, element->blockstrideout);
+           resample(output, element->kernel, element->workspace->data, element->kernel_length, element->factor, element->channels, element->blockstrideout);
 
            if (element->need_pretend) {
                element->need_pretend = FALSE;
@@ -559,9 +565,10 @@ static void finalize(GObject *object)
         /*
          * free resources
          */
-   
-   free(element->kernel);
-   free(element->workspace);
+
+   for (guint i = 0; i < element->factor; i++) 
+       gsl_vector_float_free(element->kernel[i]);
+   gsl_matrix_float_free(element->workspace);
 
         /*
          * chain to parent class' finalize() method
-- 
1.8.3.2


From dec37f9ffa00c2d4fc86b2ffeaea76ef228f94b4 Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Fri, 27 Feb 2015 14:09:36 -0500
Subject: [PATCH 5/8] gstlal_interpolator: additional vectorization

---
 gstlal/gst/lal/gstlal_interpolator.h | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/gstlal/gst/lal/gstlal_interpolator.h b/gstlal/gst/lal/gstlal_interpolator.h
index 4860c79..2d28cfa 100644
--- a/gstlal/gst/lal/gstlal_interpolator.h
+++ b/gstlal/gst/lal/gstlal_interpolator.h
@@ -34,6 +34,7 @@
 #include <gst/base/gstadapter.h>
 #include <gst/base/gstbasetransform.h>
 #include <fftw3.h>
+#include <gsl/gsl_matrix.h>
 
 #include <gstlal/gstlal.h>
 #include <gstlal/gstaudioadapter.h>
@@ -94,8 +95,9 @@ struct _GSTLALInterpolator {
    guint factor;
    guint half_length;
    guint kernel_length;
-   float *kernel;
-   float *workspace;
+   //float *kernel;
+   gsl_vector_float **kernel;
+   gsl_matrix_float *workspace;
 };
 
 
-- 
1.8.3.2


From edd27dd3429ab1f7382227bdf2853b2b0fa6c488 Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Fri, 27 Feb 2015 15:09:35 -0500
Subject: [PATCH 6/8] lal_interpolator: fix off by one

---
 gstlal/gst/lal/gstlal_interpolator.c | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/gstlal/gst/lal/gstlal_interpolator.c b/gstlal/gst/lal/gstlal_interpolator.c
index 97741ab..1c4dee8 100644
--- a/gstlal/gst/lal/gstlal_interpolator.c
+++ b/gstlal/gst/lal/gstlal_interpolator.c
@@ -88,9 +88,11 @@ gsl_vector_float** kernel(int half_length, int factor) {
        }
    }
 
-   for (int i = 0; i < 2 * half_length+1; i++) {
+   for (int i = 1; i < 2 * half_length+1; i++) {
        for (int j = 0; j < factor; j++) {
-           int index = i * factor + j + factor / 2; //FIXME FIXME FIXME this offset belongs above
+           int index = i * factor + j - factor / 2; //FIXME FIXME FIXME this offset belongs above
+           if (factor > 1)
+               index += 1;
            gsl_vector_float_set(vecs[factor - j - 1], i, out[index] / sqrt(norm / factor));
        }
    }
@@ -113,7 +115,8 @@ void resample(float *output, gsl_vector_float **thiskernel, float *input, guint
        output_offset = samp * channels;
        input_offset = (samp / factor) * channels; // first input sample
        // FIXME FIXME FIXME should this + channels really be there??
-       convolve(output + output_offset, thiskernel[kernel_offset], input + input_offset + channels, kernel_length, factor, channels);
+       //convolve(output + output_offset, thiskernel[kernel_offset], input + input_offset + channels, kernel_length, factor, channels);
+       convolve(output + output_offset, thiskernel[kernel_offset], input + input_offset, kernel_length, factor, channels);
    }
    return;
 }
-- 
1.8.3.2


From 84bee4ed8256f7054bf40c687c5ea4ea7a8bac1a Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Sat, 28 Feb 2015 09:12:37 -0500
Subject: [PATCH 7/8] gstlal_interpolator: add minimal gap support

---
 gstlal/gst/lal/gstlal_interpolator.c | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/gstlal/gst/lal/gstlal_interpolator.c b/gstlal/gst/lal/gstlal_interpolator.c
index 1c4dee8..c099256 100644
--- a/gstlal/gst/lal/gstlal_interpolator.c
+++ b/gstlal/gst/lal/gstlal_interpolator.c
@@ -108,14 +108,16 @@ void convolve(float *output, gsl_vector_float *thiskernel, float *input, guint k
    return;
 }
 
-void resample(float *output, gsl_vector_float **thiskernel, float *input, guint kernel_length, guint factor, guint channels, guint blockstrideout) {
+void resample(float *output, gsl_vector_float **thiskernel, float *input, guint kernel_length, guint factor, guint channels, guint blockstrideout, gboolean nongap) {
+   if (!nongap) {
+       memset(output, 0, sizeof(float) * blockstrideout * channels);
+       return;
+   }
    guint kernel_offset, output_offset, input_offset;
    for (guint samp = 0; samp < blockstrideout; samp++) {
        kernel_offset = samp % factor;
        output_offset = samp * channels;
        input_offset = (samp / factor) * channels; // first input sample
-       // FIXME FIXME FIXME should this + channels really be there??
-       //convolve(output + output_offset, thiskernel[kernel_offset], input + input_offset + channels, kernel_length, factor, channels);
        convolve(output + output_offset, thiskernel[kernel_offset], input + input_offset, kernel_length, factor, channels);
    }
    return;
@@ -467,6 +469,7 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
    GSTLALInterpolator *element = GSTLAL_INTERPOLATOR(trans);
    guint output_length;
    GstFlowReturn result = GST_FLOW_OK;
+   gboolean copied_nongap;
 
    g_assert(GST_BUFFER_TIMESTAMP_IS_VALID(inbuf));
    g_assert(GST_BUFFER_DURATION_IS_VALID(inbuf));
@@ -537,12 +540,12 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
            // Special case to handle discontinuities: effectively zero pad. FIXME make this more elegant
            if (element->need_pretend) {
                memset(element->workspace->data, 0, sizeof(float) * element->workspace->size1 * element->workspace->size2); // FIXME necessary?
-               gst_audioadapter_copy_samples(element->adapter, element->workspace->data + (element->half_length) * element->channels, element->blocksampsin - element->half_length, NULL, NULL);
+               gst_audioadapter_copy_samples(element->adapter, element->workspace->data + (element->half_length) * element->channels, element->blocksampsin - element->half_length, NULL, &copied_nongap);
            }
            else
                gst_audioadapter_copy_samples(element->adapter, element->workspace->data, element->blocksampsin, NULL, NULL);
 
-           resample(output, element->kernel, element->workspace->data, element->kernel_length, element->factor, element->channels, element->blockstrideout);
+           resample(output, element->kernel, element->workspace->data, element->kernel_length, element->factor, element->channels, element->blockstrideout, copied_nongap);
 
            if (element->need_pretend) {
                element->need_pretend = FALSE;
-- 
1.8.3.2


From 816cab915aeaff7ceca017c299eda894fde92432 Mon Sep 17 00:00:00 2001
From: CHAD RICHARD HANNA <crh184@psu.edu>
Date: Sat, 28 Feb 2015 10:01:15 -0500
Subject: [PATCH 8/8] gstlal_interpolator: actually allow output buffers to be
 marked as gap

---
 gstlal/gst/lal/gstlal_interpolator.c | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/gstlal/gst/lal/gstlal_interpolator.c b/gstlal/gst/lal/gstlal_interpolator.c
index c099256..c091d3e 100644
--- a/gstlal/gst/lal/gstlal_interpolator.c
+++ b/gstlal/gst/lal/gstlal_interpolator.c
@@ -286,8 +286,9 @@ static gboolean set_caps (GstBaseTransform * base, GstCaps * incaps, GstCaps * o
        free(element->kernel);
    element->kernel = kernel(element->half_length, element->factor);
 
-   // Assume that we process inrate worth of samples at a time (e.g. 1s stride)
-   element->blockstridein = element->inrate;
+   // Keep blockstride small to prevent GAPS from growing to be large
+   // FIXME probably this should be decoupled 
+   element->blockstridein = 32;//element->inrate;
    element->blocksampsin = element->blockstridein + element->kernel_length;
    element->blockstrideout = element->blockstridein * element->factor;//element->outrate;
    element->blocksampsout = element->blockstrideout + (element->kernel_length) * element->factor;
@@ -557,7 +558,7 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
            processed += element->blockstrideout;
        }
        GST_INFO_OBJECT(element, "Processed a %d samples", processed);
-       set_metadata(element, outbuf, output_length, FALSE);
+       set_metadata(element, outbuf, output_length, !copied_nongap);
 
    }
    return result;
-- 
1.8.3.2