1 files changed, 794 insertions, 0 deletions
diff --git a/gr-digital/lib/interpolating_resampler.cc b/gr-digital/lib/interpolating_resampler.cc
new file mode 100644
index 0000000000..3ec27b05de
--- /dev/null
+++ b/gr-digital/lib/interpolating_resampler.cc
@@ -0,0 +1,794 @@
+/* -*- c++ -*- */
+/*
+ * Copyright (C) 2017 Free Software Foundation, Inc.
+ *
+ * This file is part of GNU Radio
+ *
+ * GNU Radio 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 3, or (at your option)
+ * any later version.
+ *
+ * GNU Radio 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 file; see the file COPYING.  If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "interpolating_resampler.h"
+#include <gnuradio/math.h>
+#include <stdexcept>
+#include <deque>
+
+namespace gr {
+  namespace digital {
+
+    interpolating_resampler::interpolating_resampler(enum ir_type type,
+                                                     bool derivative)
+      : d_type(type),
+        d_derivative(derivative),
+        d_phase(0.0f),
+        d_phase_wrapped(0.0f),
+        d_phase_n(0),
+        d_prev_phase(0.0f),
+        d_prev_phase_wrapped(0.0f),
+        d_prev_phase_n(0)
+    {
+        switch (d_type) {
+        case IR_MMSE_8TAP:
+            break;
+        case IR_PFB_NO_MF:
+            break;
+        case IR_PFB_MF:
+            break;
+        case IR_NONE:
+        default: 
+            throw std::invalid_argument(
+              "interpolating_resampler: invalid interpolating resampler type.");
+            break;
+        }
+
+        sync_reset(0.0f);
+    }
+
+    void 
+    interpolating_resampler::next_phase(float increment,
+                                        float &phase,
+                                        int   &phase_n,
+                                        float &phase_wrapped)
+    {
+        float n;
+
+        phase = d_phase_wrapped + increment;
+        n = floorf(phase);
+        phase_wrapped = phase - n;
+        phase_n = static_cast<int>(n);
+    }
+
+    void 
+    interpolating_resampler::advance_phase(float increment)
+    {
+        d_prev_phase = d_phase;
+        d_prev_phase_wrapped = d_phase_wrapped;
+        d_prev_phase_n = d_phase_n;
+
+        next_phase(increment, d_phase, d_phase_n, d_phase_wrapped);
+    }
+
+    void
+    interpolating_resampler::revert_phase()
+    {
+        d_phase = d_prev_phase;
+        d_phase_wrapped = d_prev_phase_wrapped;
+        d_phase_n = d_prev_phase_n;
+    }
+
+    void
+    interpolating_resampler::sync_reset(float phase)
+    {
+        float n;
+
+        d_phase = phase;
+        n = floorf(d_phase);
+        d_phase_wrapped = d_phase - n;
+        d_phase_n = static_cast<int>(n);
+
+        d_prev_phase = d_phase;
+        d_prev_phase_wrapped = d_phase_wrapped;
+        d_prev_phase_n = d_phase_n;
+    }
+
+    /*************************************************************************/
+
+    interpolating_resampler_ccf *
+    interpolating_resampler_ccf::make(enum ir_type type,
+                                      bool derivative,
+                                      int nfilts,
+                                      const std::vector<float> &taps)
+    {
+        interpolating_resampler_ccf *ret = NULL;
+        switch (type) {
+        case IR_MMSE_8TAP:
+            ret = new interp_resampler_mmse_8tap_cc(derivative);
+            break;
+        case IR_PFB_NO_MF:
+            ret = new interp_resampler_pfb_no_mf_cc(derivative, nfilts);
+            break;
+        case IR_PFB_MF:
+            ret = new interp_resampler_pfb_mf_ccf(taps, nfilts, derivative);
+            break;
+        case IR_NONE:
+        default: 
+            throw std::invalid_argument("interpolating_resampler_ccf: invalid "
+                                        "interpolating resampler type.");
+            break;
+        }
+        return ret;
+    }
+
+    /*************************************************************************/
+
+    interpolating_resampler_fff *
+    interpolating_resampler_fff::make(enum ir_type type,
+                                      bool derivative,
+                                      int nfilts,
+                                      const std::vector<float> &taps)
+    {
+        interpolating_resampler_fff *ret = NULL;
+        switch (type) {
+        case IR_MMSE_8TAP:
+            ret = new interp_resampler_mmse_8tap_ff(derivative);
+            break;
+        case IR_PFB_NO_MF:
+            ret = new interp_resampler_pfb_no_mf_ff(derivative, nfilts);
+            break;
+        case IR_PFB_MF:
+            ret = new interp_resampler_pfb_mf_fff(taps, nfilts, derivative);
+            break;
+        case IR_NONE:
+        default: 
+            throw std::invalid_argument("interpolating_resampler_fff: invalid "
+                                        "interpolating resampler type.");
+            break;
+        }
+        return ret;
+    }
+
+    /*************************************************************************/
+
+    interp_resampler_mmse_8tap_cc::interp_resampler_mmse_8tap_cc(
+                                                                bool derivative)
+    : interpolating_resampler_ccf(IR_MMSE_8TAP, derivative),
+      d_interp(NULL),
+      d_interp_diff(NULL)
+    {
+        d_interp = new filter::mmse_fir_interpolator_cc();
+        if (d_interp == NULL)
+          throw std::runtime_error("unable to create mmse_fir_interpolator_cc");
+
+        if (d_derivative) {
+            d_interp_diff = new filter::mmse_interp_differentiator_cc();
+            if (d_interp_diff == NULL)
+                throw std::runtime_error("unable to create "
+                                         "mmse_interp_differentiator_cc");
+        }
+    }
+
+    interp_resampler_mmse_8tap_cc::~interp_resampler_mmse_8tap_cc()
+    {
+        delete d_interp;
+        if (d_derivative)
+            delete d_interp_diff;
+    }
+
+    gr_complex
+    interp_resampler_mmse_8tap_cc::interpolate(const gr_complex input[],
+                                               float mu) const
+    {
+        return d_interp->interpolate(input, mu);
+    }
+
+    gr_complex
+    interp_resampler_mmse_8tap_cc::differentiate(const gr_complex input[],
+                                                 float mu) const
+    {
+        return d_interp_diff->differentiate(input, mu);
+    }
+
+    unsigned int
+    interp_resampler_mmse_8tap_cc::ntaps() const
+    {
+        return d_interp->ntaps();
+    }
+
+    /*************************************************************************/
+
+    interp_resampler_mmse_8tap_ff::interp_resampler_mmse_8tap_ff(
+                                                                bool derivative)
+    : interpolating_resampler_fff(IR_MMSE_8TAP, derivative),
+      d_interp(NULL),
+      d_interp_diff(NULL)
+    {
+        d_interp = new filter::mmse_fir_interpolator_ff();
+        if (d_interp == NULL)
+          throw std::runtime_error("unable to create mmse_fir_interpolator_ff");
+
+        if (d_derivative) {
+            d_interp_diff = new filter::mmse_interp_differentiator_ff();
+            if (d_interp_diff == NULL)
+                throw std::runtime_error("unable to create "
+                                         "mmse_interp_differentiator_ff");
+        }
+    }
+
+    interp_resampler_mmse_8tap_ff::~interp_resampler_mmse_8tap_ff()
+    {
+        delete d_interp;
+        if (d_derivative)
+            delete d_interp_diff;
+    }
+
+    float
+    interp_resampler_mmse_8tap_ff::interpolate(const float input[],
+                                               float mu) const
+    {
+        return d_interp->interpolate(input, mu);
+    }
+
+    float
+    interp_resampler_mmse_8tap_ff::differentiate(const float input[],
+                                                 float mu) const
+    {
+        return d_interp_diff->differentiate(input, mu);
+    }
+
+    unsigned int
+    interp_resampler_mmse_8tap_ff::ntaps() const
+    {
+        return d_interp->ntaps();
+    }
+
+    /*************************************************************************/
+
+#include "gnuradio/filter/interpolator_taps.h"
+#include "gnuradio/filter/interp_differentiator_taps.h"
+
+    interp_resampler_pfb_no_mf_cc::interp_resampler_pfb_no_mf_cc(
+                                                                bool derivative,
+                                                                int nfilts)
+    : interpolating_resampler_ccf(IR_PFB_NO_MF, derivative),
+      d_nfilters(0),
+      d_filters(),
+      d_diff_filters()
+    {
+        if (nfilts <= 1)
+            throw std::invalid_argument("interpolating_resampler_pfb_no_mf_cc: "
+                                        "number of polyphase filter arms "
+                                        "must be greater than 1");
+
+        // Round up the number of filter arms to the current or next power of 2
+        d_nfilters = 
+               1 << (static_cast<int>(log2f(static_cast<float>(nfilts-1))) + 1);
+
+        // N.B. We assume in this class: NSTEPS == DNSTEPS and NTAPS == DNTAPS
+
+        // Limit to the maximum number of precomputed MMSE tap sets
+        if (d_nfilters <= 0 or d_nfilters > NSTEPS)
+            d_nfilters = NSTEPS;
+
+        // Create our polyphase filter arms for the steps from 0.0 to 1.0 from
+        // the MMSE interpolating filter and MMSE interpolating differentiator
+        // taps rows.
+        // N.B. We create an extra final row for an offset of 1.0, because it's
+        // easier than dealing with wrap around from 0.99... to 0.0 shifted
+        // by 1 tap.
+        d_filters = std::vector<filter::kernel::fir_filter_ccf*>(d_nfilters + 1,
+                                                                 NULL);
+        d_diff_filters = std::vector<filter::kernel::fir_filter_ccf*>(
+                                                                 d_nfilters + 1,
+                                                                 NULL);
+
+        std::vector<float> t(NTAPS, 0);
+        int incr = NSTEPS/d_nfilters;
+        int src, dst;
+        for (src = 0, dst = 0; src <= NSTEPS; src += incr, dst++) {
+
+            t.assign(&taps[src][0], &taps[src][NTAPS]);
+            d_filters[dst] = new filter::kernel::fir_filter_ccf(1, t);
+            if (d_filters[dst] == NULL)
+                throw std::runtime_error("unable to create fir_filter_ccf");
+
+            if (d_derivative) {
+                t.assign(&Dtaps[src][0], &Dtaps[src][DNTAPS]);
+                d_diff_filters[dst] = new filter::kernel::fir_filter_ccf(1, t);
+                if (d_diff_filters[dst] == NULL)
+                    throw std::runtime_error("unable to create fir_filter_ccf");
+            }
+        }
+    }
+
+    interp_resampler_pfb_no_mf_cc::~interp_resampler_pfb_no_mf_cc()
+    {
+        for (int i = 0; i <= d_nfilters; i++) {
+            delete d_filters[i];
+            if (d_derivative)
+                delete d_diff_filters[i];
+        }
+    }
+
+    gr_complex
+    interp_resampler_pfb_no_mf_cc::interpolate(const gr_complex input[],
+                                               float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_no_mf_cc: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_filters[arm]->filter(input);
+    }
+
+    gr_complex
+    interp_resampler_pfb_no_mf_cc::differentiate(const gr_complex input[],
+                                                 float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_no_mf_cc: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_diff_filters[arm]->filter(input);
+    }
+
+    unsigned int
+    interp_resampler_pfb_no_mf_cc::ntaps() const
+    {
+        return NTAPS;
+    }
+
+    /*************************************************************************/
+
+    interp_resampler_pfb_no_mf_ff::interp_resampler_pfb_no_mf_ff(
+                                                                bool derivative,
+                                                                int nfilts)
+    : interpolating_resampler_fff(IR_PFB_NO_MF, derivative),
+      d_nfilters(0),
+      d_filters(),
+      d_diff_filters()
+    {
+        if (nfilts <= 1)
+            throw std::invalid_argument("interpolating_resampler_pfb_no_mf_ff: "
+                                        "number of polyphase filter arms "
+                                        "must be greater than 1");
+
+        // Round up the number of filter arms to the current or next power of 2
+        d_nfilters = 
+               1 << (static_cast<int>(log2f(static_cast<float>(nfilts-1))) + 1);
+
+        // N.B. We assume in this class: NSTEPS == DNSTEPS and NTAPS == DNTAPS
+
+        // Limit to the maximum number of precomputed MMSE tap sets
+        if (d_nfilters <= 0 or d_nfilters > NSTEPS)
+            d_nfilters = NSTEPS;
+
+        // Create our polyphase filter arms for the steps from 0.0 to 1.0 from
+        // the MMSE interpolating filter and MMSE interpolating differentiator
+        // taps rows.
+        // N.B. We create an extra final row for an offset of 1.0, because it's
+        // easier than dealing with wrap around from 0.99... to 0.0 shifted
+        // by 1 tap.
+        d_filters = std::vector<filter::kernel::fir_filter_fff*>(d_nfilters + 1,
+                                                                 NULL);
+        d_diff_filters = std::vector<filter::kernel::fir_filter_fff*>(
+                                                                 d_nfilters + 1,
+                                                                 NULL);
+
+        std::vector<float> t(NTAPS, 0);
+        int incr = NSTEPS/d_nfilters;
+        int src, dst;
+        for (src = 0, dst = 0; src <= NSTEPS; src += incr, dst++) {
+
+            t.assign(&taps[src][0], &taps[src][NTAPS]);
+            d_filters[dst] = new filter::kernel::fir_filter_fff(1, t);
+            if (d_filters[dst] == NULL)
+                throw std::runtime_error("unable to create fir_filter_fff");
+
+            if (d_derivative) {
+                t.assign(&Dtaps[src][0], &Dtaps[src][DNTAPS]);
+                d_diff_filters[dst] = new filter::kernel::fir_filter_fff(1, t);
+                if (d_diff_filters[dst] == NULL)
+                    throw std::runtime_error("unable to create fir_filter_fff");
+            }
+        }
+    }
+
+    interp_resampler_pfb_no_mf_ff::~interp_resampler_pfb_no_mf_ff()
+    {
+        for (int i = 0; i <= d_nfilters; i++) {
+            delete d_filters[i];
+            if (d_derivative)
+                delete d_diff_filters[i];
+        }
+    }
+
+    float
+    interp_resampler_pfb_no_mf_ff::interpolate(const float input[],
+                                               float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_no_mf_ff: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_filters[arm]->filter(input);
+    }
+
+    float
+    interp_resampler_pfb_no_mf_ff::differentiate(const float input[],
+                                                 float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_no_mf_ff: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_diff_filters[arm]->filter(input);
+    }
+
+    unsigned int
+    interp_resampler_pfb_no_mf_ff::ntaps() const
+    {
+        return NTAPS;
+    }
+
+    /*************************************************************************/
+
+    interp_resampler_pfb_mf_ccf::interp_resampler_pfb_mf_ccf(
+                                                 const std::vector<float> &taps,
+                                                 int nfilts,
+                                                 bool derivative)
+    : interpolating_resampler_ccf(IR_PFB_MF, derivative),
+      d_nfilters(nfilts),
+      d_taps_per_filter(static_cast<unsigned int>(
+                                     ceil(  static_cast<double>(taps.size())
+                                          / static_cast<double>(nfilts     )))),
+      d_filters(),
+      d_diff_filters(),
+      d_taps(),
+      d_diff_taps()
+    {
+        if (d_nfilters <= 1)
+            throw std::invalid_argument("interpolating_resampler_pfb_mf_ccf: "
+                                        "number of polyphase filter arms "
+                                        "must be greater than 1");
+        if (taps.size() < static_cast<unsigned int>(d_nfilters))
+            throw std::invalid_argument("interpolating_resampler_pfb_mf_ccf: "
+                                        "length of the prototype filter taps"
+                                        " must be greater than or equal to "
+                                        "the number of polyphase filter arms.");
+
+        // Create a derivative filter from the provided taps
+
+        // First create a truncated ideal differentiator filter
+        int ideal_diff_filt_len = 3; // Must be odd; rest of init assumes odd.
+        std::vector<float> ideal_diff_taps(ideal_diff_filt_len, 0.0f);
+        int i, n;
+        n = ideal_diff_taps.size()/2;
+        for (i = -n; i < 0; i++) {
+            ideal_diff_taps[i+n] = (-i & 1) == 1 ? -1.0f/i : 1.0f/i;
+            ideal_diff_taps[n-i] = -ideal_diff_taps[i+n];
+        }
+        ideal_diff_taps[n] = 0.0f;
+
+        // Perform linear convolution of prototype filter taps and the truncated
+        // ideal differentiator taps to generate a derivative matched filter.
+        // N.B. the truncated ideal differentiator taps must have an odd length
+        int j, k, l, m;
+        m = ideal_diff_taps.size();
+        n = taps.size();
+        l = m + n - 1; // length of convolution
+        std::deque<float> diff_taps(l, 0.0f);
+        for (i = 0; i < l; i++) {
+            for (j = 0; j < m;  j++) {
+                k = i + j - (m - 1);
+                if (k < 0 or k >= n)
+                    continue;
+                diff_taps[i] += ideal_diff_taps[(m - 1) - j] * taps[k];
+            }
+        }
+
+        // Trim the convolution so the prototype derivative filter is the same
+        // length as the passed in prototype filter taps.
+        n = ideal_diff_taps.size()/2;
+        for (i = 0; i < n; i++) {
+            diff_taps.pop_back();
+            diff_taps.pop_front();
+        }
+
+        // Squash the differentiation noise spikes at the filter ends.
+        diff_taps[0] = 0.0f;
+        diff_taps[diff_taps.size()-1] = 0.0f;
+
+        // Normalize the prototype derviative filter gain to the number of
+        // filter arms
+        n = diff_taps.size();
+        float mag = 0.0f;
+        for (i = 0; i < n; i++)
+            mag += fabsf(diff_taps[i]);
+        for (i = 0; i < n; i++) {
+            diff_taps[i] *= d_nfilters/mag;
+            if (d_derivative and std::isnan(diff_taps[i]))
+                throw std::runtime_error("interpolating_resampler_pfb_mf_ccf: "
+                                         "NaN error creating derivative filter."
+                                        ); 
+        }
+
+        // Create our polyphase filter arms for the steps from 0.0 to 1.0 from
+        // the prototype matched filter.
+        // N.B. We create an extra final row for an offset of 1.0, because it's
+        // easier than dealing with wrap around from 0.99... to 0.0 shifted
+        // by 1 tap.
+        d_filters = std::vector<filter::kernel::fir_filter_ccf*>(d_nfilters + 1,
+                                                                 NULL);
+        d_diff_filters = std::vector<filter::kernel::fir_filter_ccf*>(
+                                                                 d_nfilters + 1,
+                                                                 NULL);
+
+        m = taps.size();
+        n = diff_taps.size();
+        d_taps.resize(d_nfilters+1);
+        d_diff_taps.resize(d_nfilters+1);
+        signed int taps_per_filter = static_cast<signed int>(d_taps_per_filter);
+
+        for (i = 0; i <= d_nfilters; i++) {
+            d_taps[i] = std::vector<float>(d_taps_per_filter, 0.0f);
+            for (j = 0; j < taps_per_filter; j++) {
+                k = i+j*d_nfilters;
+                if (k < m)
+                    d_taps[i][j] = taps[k];
+            }
+            d_filters[i] = new filter::kernel::fir_filter_ccf(1, d_taps[i]);
+            if (d_filters[i] == NULL)
+                throw std::runtime_error("unable to create fir_filter_ccf");
+
+            if (not d_derivative)
+                continue;
+
+            d_diff_taps[i] = std::vector<float>(d_taps_per_filter, 0.0f);
+            for (j = 0; j < taps_per_filter; j++) {
+                k = i+j*d_nfilters;
+                if (k < n)
+                    d_diff_taps[i][j] = diff_taps[k];
+            }
+            d_diff_filters[i] = new filter::kernel::fir_filter_ccf(
+                                                                1,
+                                                                d_diff_taps[i]);
+            if (d_diff_filters[i] == NULL)
+                throw std::runtime_error("unable to create fir_filter_ccf");
+        }
+    }
+
+    interp_resampler_pfb_mf_ccf::~interp_resampler_pfb_mf_ccf()
+    {
+        for (int i = 0; i <= d_nfilters; i++) {
+            delete d_filters[i];
+            if (d_derivative)
+                delete d_diff_filters[i];
+        }
+    }
+
+    gr_complex
+    interp_resampler_pfb_mf_ccf::interpolate(const gr_complex input[],
+                                             float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_mf_ccf: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_filters[arm]->filter(input);
+    }
+
+    gr_complex
+    interp_resampler_pfb_mf_ccf::differentiate(const gr_complex input[],
+                                               float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_mf_ccf: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_diff_filters[arm]->filter(input);
+    }
+
+    unsigned int
+    interp_resampler_pfb_mf_ccf::ntaps() const
+    {
+        return d_taps_per_filter;
+    }
+
+    /*************************************************************************/
+
+    interp_resampler_pfb_mf_fff::interp_resampler_pfb_mf_fff(
+                                                 const std::vector<float> &taps,
+                                                 int nfilts,
+                                                 bool derivative)
+    : interpolating_resampler_fff(IR_PFB_MF, derivative),
+      d_nfilters(nfilts),
+      d_taps_per_filter(static_cast<unsigned int>(
+                                     ceil(  static_cast<double>(taps.size())
+                                          / static_cast<double>(nfilts     )))),
+      d_filters(),
+      d_diff_filters(),
+      d_taps(),
+      d_diff_taps()
+    {
+        if (d_nfilters <= 1)
+            throw std::invalid_argument("interpolating_resampler_pfb_mf_fff: "
+                                        "number of polyphase filter arms "
+                                        "must be greater than 1");
+        if (taps.size() < static_cast<unsigned int>(d_nfilters))
+            throw std::invalid_argument("interpolating_resampler_pfb_mf_fff: "
+                                        "length of the prototype filter taps"
+                                        " must be greater than or equal to "
+                                        "the number of polyphase filter arms.");
+
+        // Create a derivative filter from the provided taps
+
+        // First create a truncated ideal differentiator filter
+        int ideal_diff_filt_len = 3; // Must be odd; rest of init assumes odd.
+        std::vector<float> ideal_diff_taps(ideal_diff_filt_len, 0.0f);
+        int i, n;
+        n = ideal_diff_taps.size()/2;
+        for (i = -n; i < 0; i++) {
+            ideal_diff_taps[i+n] = (-i & 1) == 1 ? -1.0f/i : 1.0f/i;
+            ideal_diff_taps[n-i] = -ideal_diff_taps[i+n];
+        }
+        ideal_diff_taps[n] = 0.0f;
+
+        // Perform linear convolution of prototype filter taps and the truncated
+        // ideal differentiator taps to generate a derivative matched filter.
+        // N.B. the truncated ideal differentiator taps must have an odd length
+        int j, k, l, m;
+        m = ideal_diff_taps.size();
+        n = taps.size();
+        l = m + n - 1; // length of convolution
+        std::deque<float> diff_taps(l, 0.0f);
+        for (i = 0; i < l; i++) {
+            for (j = 0; j < m;  j++) {
+                k = i + j - (m - 1);
+                if (k < 0 or k >= n)
+                    continue;
+                diff_taps[i] += ideal_diff_taps[(m - 1) - j] * taps[k];
+            }
+        }
+
+        // Trim the convolution so the prototype derivative filter is the same
+        // length as the passed in prototype filter taps.
+        n = ideal_diff_taps.size()/2;
+        for (i = 0; i < n; i++) {
+            diff_taps.pop_back();
+            diff_taps.pop_front();
+        }
+
+        // Squash the differentiation noise spikes at the filter ends.
+        diff_taps[0] = 0.0f;
+        diff_taps[diff_taps.size()-1] = 0.0f;
+
+        // Normalize the prototype derviative filter gain to the number of
+        // filter arms
+        n = diff_taps.size();
+        float mag = 0.0f;
+        for (i = 0; i < n; i++)
+            mag += fabsf(diff_taps[i]);
+        for (i = 0; i < n; i++) {
+            diff_taps[i] *= d_nfilters/mag;
+            if (d_derivative and std::isnan(diff_taps[i]))
+                throw std::runtime_error("interpolating_resampler_pfb_mf_fff: "
+                                         "NaN error creating derivative filter."
+                                        ); 
+        }
+
+        // Create our polyphase filter arms for the steps from 0.0 to 1.0 from
+        // the prototype matched filter.
+        // N.B. We create an extra final row for an offset of 1.0, because it's
+        // easier than dealing with wrap around from 0.99... to 0.0 shifted
+        // by 1 tap.
+        d_filters = std::vector<filter::kernel::fir_filter_fff*>(d_nfilters + 1,
+                                                                 NULL);
+        d_diff_filters = std::vector<filter::kernel::fir_filter_fff*>(
+                                                                 d_nfilters + 1,
+                                                                 NULL);
+
+        m = taps.size();
+        n = diff_taps.size();
+        d_taps.resize(d_nfilters+1);
+        d_diff_taps.resize(d_nfilters+1);
+        signed int taps_per_filter = static_cast<signed int>(d_taps_per_filter);
+
+        for (i = 0; i <= d_nfilters; i++) {
+            d_taps[i] = std::vector<float>(d_taps_per_filter, 0.0f);
+            for (j = 0; j < taps_per_filter; j++) {
+                k = i+j*d_nfilters;
+                if (k < m)
+                    d_taps[i][j] = taps[k];
+            }
+            d_filters[i] = new filter::kernel::fir_filter_fff(1, d_taps[i]);
+            if (d_filters[i] == NULL)
+                throw std::runtime_error("unable to create fir_filter_fff");
+
+            if (not d_derivative)
+                continue;
+
+            d_diff_taps[i] = std::vector<float>(d_taps_per_filter, 0.0f);
+            for (j = 0; j < taps_per_filter; j++) {
+                k = i+j*d_nfilters;
+                if (k < n)
+                    d_diff_taps[i][j] = diff_taps[k];
+            }
+            d_diff_filters[i] = new filter::kernel::fir_filter_fff(
+                                                                1,
+                                                                d_diff_taps[i]);
+            if (d_diff_filters[i] == NULL)
+                throw std::runtime_error("unable to create fir_filter_fff");
+        }
+    }
+
+    interp_resampler_pfb_mf_fff::~interp_resampler_pfb_mf_fff()
+    {
+        for (int i = 0; i <= d_nfilters; i++) {
+            delete d_filters[i];
+            if (d_derivative)
+                delete d_diff_filters[i];
+        }
+    }
+
+    float
+    interp_resampler_pfb_mf_fff::interpolate(const float input[],
+                                             float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_mf_fff: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_filters[arm]->filter(input);
+    }
+
+    float
+    interp_resampler_pfb_mf_fff::differentiate(const float input[],
+                                               float mu) const
+    {
+        int arm = static_cast<int>(rint(mu * d_nfilters));
+
+        if (arm < 0 or arm > d_nfilters)
+            throw std::runtime_error("interp_resampler_pfb_mf_fff: mu is not "
+                                     "in the range [0.0, 1.0]");
+
+        return d_diff_filters[arm]->filter(input);
+    }
+
+    unsigned int
+    interp_resampler_pfb_mf_fff::ntaps() const
+    {
+        return d_taps_per_filter;
+    }
+
+  } /* namespace digital */
+} /* namespace gr */