9 files changed, 681 insertions, 244 deletions
diff --git a/gr-filter/include/gnuradio/filter/pfb_arb_resampler.h b/gr-filter/include/gnuradio/filter/pfb_arb_resampler.h
index 69bc5408bc..47f55ad7e5 100644
--- a/gr-filter/include/gnuradio/filter/pfb_arb_resampler.h
+++ b/gr-filter/include/gnuradio/filter/pfb_arb_resampler.h
@@ -22,7 +22,7 @@
 
 
 #ifndef INCLUDED_PFB_ARB_RESAMPLER_H
-#define	INCLUDED_PFB_ARB_RESAMPLER__H
+#define	INCLUDED_PFB_ARB_RESAMPLER_H
 
 #include <gnuradio/filter/fir_filter.h>
 
@@ -104,6 +104,8 @@ namespace gr {
         float        d_acc;               // accumulator; holds fractional part of sample
         unsigned int d_last_filter;       // stores filter for re-entry
         unsigned int d_taps_per_filter;   // num taps for each arm of the filterbank
+        int d_delay;                      // filter's group delay
+        float d_est_phase_change;         // est. of phase change of a sine wave through filt.
 
         /*!
          * Takes in the taps and convolves them with [-1,0,1], which
@@ -177,6 +179,22 @@ namespace gr {
          */
         unsigned int taps_per_filter() const;
 
+        unsigned int interpolation_rate() const { return d_int_rate; }
+        unsigned int decimation_rate() const { return d_dec_rate; }
+        float fractional_rate() const { return d_flt_rate; }
+
+        /*!
+         * Get the group delay of the filter.
+         */
+        int group_delay() const { return d_delay; }
+
+        /*!
+         * Calculates the phase offset expected by a sine wave of
+         * frequency \p freq and sampling rate \p fs (assuming input
+         * sine wave has 0 degree phase).
+         */
+        float phase_offset(float freq, float fs);
+
         /*!
          * Performs the filter operation that resamples the signal.
          *
@@ -195,8 +213,195 @@ namespace gr {
                    int n_to_read, int &n_read);
       };
 
+
+      /**************************************************************/
+
+
+      /*!
+       * \brief Polyphase filterbank arbitrary resampler with
+       *        float input, float output and float taps
+       * \ingroup resamplers_blk
+       *
+       * \details
+       * This  takes in a signal stream and performs arbitrary
+       * resampling. The resampling rate can be any real number
+       * <EM>r</EM>. The resampling is done by constructing <EM>N</EM>
+       * filters where <EM>N</EM> is the interpolation rate.  We then
+       * calculate <EM>D</EM> where <EM>D = floor(N/r)</EM>.
+       *
+       * Using <EM>N</EM> and <EM>D</EM>, we can perform rational
+       * resampling where <EM>N/D</EM> is a rational number close to
+       * the input rate <EM>r</EM> where we have <EM>N</EM> filters
+       * and we cycle through them as a polyphase filterbank with a
+       * stride of <EM>D</EM> so that <EM>i+1 = (i + D) % N</EM>.
+       *
+       * To get the arbitrary rate, we want to interpolate between two
+       * points. For each value out, we take an output from the
+       * current filter, <EM>i</EM>, and the next filter <EM>i+1</EM>
+       * and then linearly interpolate between the two based on the
+       * real resampling rate we want.
+       *
+       * The linear interpolation only provides us with an
+       * approximation to the real sampling rate specified. The error
+       * is a quantization error between the two filters we used as
+       * our interpolation points.  To this end, the number of
+       * filters, <EM>N</EM>, used determines the quantization error;
+       * the larger <EM>N</EM>, the smaller the noise. You can design
+       * for a specified noise floor by setting the filter size
+       * (parameters <EM>filter_size</EM>). The size defaults to 32
+       * filters, which is about as good as most implementations need.
+       *
+       * The trick with designing this filter is in how to specify the
+       * taps of the prototype filter. Like the PFB interpolator, the
+       * taps are specified using the interpolated filter rate. In
+       * this case, that rate is the input sample rate multiplied by
+       * the number of filters in the filterbank, which is also the
+       * interpolation rate. All other values should be relative to
+       * this rate.
+       *
+       * For example, for a 32-filter arbitrary resampler and using
+       * the GNU Radio's firdes utility to build the filter, we build
+       * a low-pass filter with a sampling rate of <EM>fs</EM>, a 3-dB
+       * bandwidth of <EM>BW</EM> and a transition bandwidth of
+       * <EM>TB</EM>. We can also specify the out-of-band attenuation
+       * to use, <EM>ATT</EM>, and the filter window function (a
+       * Blackman-harris window in this case). The first input is the
+       * gain of the filter, which we specify here as the
+       * interpolation rate (<EM>32</EM>).
+       *
+       *   <B><EM>self._taps = filter.firdes.low_pass_2(32, 32*fs, BW, TB,
+       *      attenuation_dB=ATT, window=filter.firdes.WIN_BLACKMAN_hARRIS)</EM></B>
+       *
+       * The theory behind this block can be found in Chapter 7.5 of
+       * the following book.
+       *
+       *   <B><EM>f. harris, "Multirate Signal Processing for Communication
+       *      Systems", Upper Saddle River, NJ: Prentice Hall, Inc. 2004.</EM></B>
+       */
+      class FILTER_API pfb_arb_resampler_fff
+      {
+      private:
+        std::vector<fir_filter_fff*> d_filters;
+        std::vector<fir_filter_fff*> d_diff_filters;
+        std::vector< std::vector<float> > d_taps;
+        std::vector< std::vector<float> > d_dtaps;
+        unsigned int d_int_rate;          // the number of filters (interpolation rate)
+        unsigned int d_dec_rate;          // the stride through the filters (decimation rate)
+        float        d_flt_rate;          // residual rate for the linear interpolation
+        float        d_acc;               // accumulator; holds fractional part of sample
+        unsigned int d_last_filter;       // stores filter for re-entry
+        unsigned int d_taps_per_filter;   // num taps for each arm of the filterbank
+        int d_delay;                      // filter's group delay
+        float d_est_phase_change;         // est. of phase change of a sine wave through filt.
+
+        /*!
+         * Takes in the taps and convolves them with [-1,0,1], which
+         * creates a differential set of taps that are used in the
+         * difference filterbank.
+         * \param newtaps (vector of floats) The prototype filter.
+         * \param difftaps (vector of floats) (out) The differential filter taps.
+         */
+        void create_diff_taps(const std::vector<float> &newtaps,
+                              std::vector<float> &difftaps);
+
+        /*!
+         * Resets the filterbank's filter taps with the new prototype filter
+         * \param newtaps    (vector of floats) The prototype filter to populate the filterbank.
+         *                   The taps should be generated at the interpolated sampling rate.
+         * \param ourtaps    (vector of floats) Reference to our internal member of holding the taps.
+         * \param ourfilter  (vector of filters) Reference to our internal filter to set the taps for.
+         */
+        void create_taps(const std::vector<float> &newtaps,
+                         std::vector< std::vector<float> > &ourtaps,
+                         std::vector<kernel::fir_filter_fff*> &ourfilter);
+
+      public:
+        /*!
+         * Creates a kernel to perform arbitrary resampling on a set of samples.
+         * \param rate  (float) Specifies the resampling rate to use
+         * \param taps  (vector/list of floats) The prototype filter to populate the filterbank. The taps       *              should be generated at the filter_size sampling rate.
+         * \param filter_size (unsigned int) The number of filters in the filter bank. This is directly
+         *                    related to quantization noise introduced during the resampling.
+         *                    Defaults to 32 filters.
+         */
+        pfb_arb_resampler_fff(float rate,
+                              const std::vector<float> &taps,
+                              unsigned int filter_size);
+
+        ~pfb_arb_resampler_fff();
+
+        /*!
+         * Resets the filterbank's filter taps with the new prototype filter
+         * \param taps (vector/list of floats) The prototype filter to populate the filterbank.
+         */
+        void set_taps(const std::vector<float> &taps);
+
+        /*!
+         * Return a vector<vector<>> of the filterbank taps
+         */
+        std::vector<std::vector<float> > taps() const;
+
+        /*!
+         * Print all of the filterbank taps to screen.
+         */
+        void print_taps();
+
+        /*!
+         * Sets the resampling rate of the block.
+         */
+        void set_rate(float rate);
+
+        /*!
+         * Sets the current phase offset in radians (0 to 2pi).
+         */
+        void set_phase(float ph);
+
+        /*!
+         * Gets the current phase of the resampler in radians (2 to 2pi).
+         */
+        float phase() const;
+
+        /*!
+         * Gets the number of taps per filter.
+         */
+        unsigned int taps_per_filter() const;
+
+        unsigned int interpolation_rate() const { return d_int_rate; }
+        unsigned int decimation_rate() const { return d_dec_rate; }
+        float fractional_rate() const { return d_flt_rate; }
+
+        /*!
+         * Get the group delay of the filter.
+         */
+        int group_delay() const { return d_delay; }
+
+        /*!
+         * Calculates the phase offset expected by a sine wave of
+         * frequency \p freq and sampling rate \p fs (assuming input
+         * sine wave has 0 degree phase).
+         */
+        float phase_offset(float freq, float fs);
+
+        /*!
+         * Performs the filter operation that resamples the signal.
+         *
+         * This block takes in a stream of samples and outputs a
+         * resampled and filtered stream. This block should be called
+         * such that the output has \p rate * \p n_to_read amount of
+         * space available in the \p output buffer.
+         *
+         * \param input An input vector of samples to be resampled
+         * \param output The output samples at the new sample rate.
+         * \param n_to_read Number of samples to read from \p input.
+         * \param n_read (out) Number of samples actually read from \p input.
+         * \return Number of samples put into \p output.
+         */
+        int filter(float *input, float *output,
+                   int n_to_read, int &n_read);
+      };
+
     } /* namespace kernel */
   } /* namespace filter */
 } /* namespace gr */
 
-#endif /* INCLUDED_PFB_ARB_RESAMPLER_CCF_IMPL_H */
+#endif /* INCLUDED_PFB_ARB_RESAMPLER_H */
diff --git a/gr-filter/include/gnuradio/filter/pfb_arb_resampler_ccf.h b/gr-filter/include/gnuradio/filter/pfb_arb_resampler_ccf.h
index a519b2d096..1c5dd49628 100644
--- a/gr-filter/include/gnuradio/filter/pfb_arb_resampler_ccf.h
+++ b/gr-filter/include/gnuradio/filter/pfb_arb_resampler_ccf.h
@@ -97,6 +97,33 @@ namespace gr {
        * Gets the number of taps per filter.
        */
       virtual unsigned int taps_per_filter() const = 0;
+
+      /*!
+       * Gets the interpolation rate of the filter.
+       */
+      virtual unsigned int interpolation_rate() const = 0;
+
+      /*!
+       * Gets the decimation rate of the filter.
+       */    
+      virtual unsigned int decimation_rate() const =0;
+      
+      /*!
+       * Gets the fractional rate of the filter.
+       */    
+      virtual float fractional_rate() const = 0;
+
+      /*!
+       * Get the group delay of the filter.
+       */
+      virtual int group_delay() const = 0;
+
+      /*!
+       * Calculates the phase offset expected by a sine wave of
+       * frequency \p freq and sampling rate \p fs (assuming input
+       * sine wave has 0 degree phase).
+       */
+      virtual float phase_offset(float freq, float fs) = 0;
     };
 
   } /* namespace filter */
diff --git a/gr-filter/include/gnuradio/filter/pfb_arb_resampler_fff.h b/gr-filter/include/gnuradio/filter/pfb_arb_resampler_fff.h
index c28e1dc2a6..70016bc656 100644
--- a/gr-filter/include/gnuradio/filter/pfb_arb_resampler_fff.h
+++ b/gr-filter/include/gnuradio/filter/pfb_arb_resampler_fff.h
@@ -108,8 +108,8 @@ namespace gr {
        *			           Defaults to 32 filters.
        */
       static sptr make(float rate,
-				  const std::vector<float> &taps,
-				  unsigned int filter_size=32);
+                       const std::vector<float> &taps,
+                       unsigned int filter_size=32);
 
       /*!
        * Resets the filterbank's filter taps with the new prototype filter
@@ -141,6 +141,38 @@ namespace gr {
        * Gets the current phase of the resampler in radians (2 to 2pi).
        */
       virtual float phase() const = 0;
+
+      /*!
+       * Gets the number of taps per filter.
+       */
+      virtual unsigned int taps_per_filter() const = 0;
+
+      /*!
+       * Gets the interpolation rate of the filter.
+       */
+      virtual unsigned int interpolation_rate() const = 0;
+
+      /*!
+       * Gets the decimation rate of the filter.
+       */    
+      virtual unsigned int decimation_rate() const =0;
+      
+      /*!
+       * Gets the fractional rate of the filter.
+       */    
+      virtual float fractional_rate() const = 0;
+
+      /*!
+       * Get the group delay of the filter.
+       */
+      virtual int group_delay() const = 0;
+
+      /*!
+       * Calculates the phase offset expected by a sine wave of
+       * frequency \p freq and sampling rate \p fs (assuming input
+       * sine wave has 0 degree phase).
+       */
+      virtual float phase_offset(float freq, float fs) = 0;
     };
 
   } /* namespace filter */
diff --git a/gr-filter/lib/pfb_arb_resampler.cc b/gr-filter/lib/pfb_arb_resampler.cc
index 7e92839296..9f66276d70 100644
--- a/gr-filter/lib/pfb_arb_resampler.cc
+++ b/gr-filter/lib/pfb_arb_resampler.cc
@@ -28,6 +28,7 @@
 #include <gnuradio/logger.h>
 #include <cstdio>
 #include <stdexcept>
+#include <boost/math/special_functions/round.hpp>
 
 namespace gr {
   namespace filter {
@@ -51,7 +52,7 @@ namespace gr {
         d_int_rate = filter_size;
         set_rate(rate);
 
-        d_last_filter = 0;
+        d_last_filter = (taps.size()/2) % filter_size;
         
         d_filters = std::vector<fir_filter_ccf*>(d_int_rate);
         d_diff_filters = std::vector<fir_filter_ccf*>(d_int_rate);
@@ -65,6 +66,25 @@ namespace gr {
 
         // Now, actually set the filters' taps
         set_taps(taps);
+
+        // Delay is based on number of taps per filter arm. Round to
+        // the nearest integer.
+        float delay = -rate * (taps_per_filter() - 1.0) / 2.0;
+        d_delay = static_cast<int>(boost::math::iround(delay));
+
+        // This calculation finds the phase offset induced by the
+        // arbitrary resampling. It's based on which filter arm we are
+        // at the filter's group delay plus the fractional offset
+        // between the samples. Calculated here based on the rotation
+        // around nfilts starting at start_filter.
+        float accum = -d_delay * d_flt_rate;
+        int   accum_int = static_cast<int>(accum);
+        float accum_frac = accum - accum_int;
+        int end_filter = static_cast<int>
+          (boost::math::iround(fmodf(d_last_filter - d_delay * d_dec_rate + accum_int, \
+                                     static_cast<float>(d_int_rate))));
+
+        d_est_phase_change = d_last_filter - (end_filter + accum_frac);
       }
 
       pfb_arb_resampler_ccf::~pfb_arb_resampler_ccf()
@@ -94,33 +114,36 @@ namespace gr {
           tmp_taps.push_back(0.0);
         }
 
-        // Partition the filter
         for(unsigned int i = 0; i < d_int_rate; i++) {
           // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
-          ourtaps[d_int_rate-1-i] = std::vector<float>(d_taps_per_filter, 0);
+          ourtaps[i] = std::vector<float>(d_taps_per_filter, 0);
           for(unsigned int j = 0; j < d_taps_per_filter; j++) {
-            ourtaps[d_int_rate - 1 - i][j] = tmp_taps[i + j*d_int_rate];
+            ourtaps[i][j] = tmp_taps[i + j*d_int_rate];
           }
 
           // Build a filter for each channel and add it's taps to it
-          ourfilter[i]->set_taps(ourtaps[d_int_rate-1-i]);
+          ourfilter[i]->set_taps(ourtaps[i]);
         }
       }
 
       void
       pfb_arb_resampler_ccf::create_diff_taps(const std::vector<float> &newtaps,
-                                                   std::vector<float> &difftaps)
+                                              std::vector<float> &difftaps)
       {
-        // Calculate the differential taps (derivative filter) by
-        // taking the difference between two taps. Duplicate the last
-        // one to make both filters the same length.
-        float tap;
-        difftaps.clear();
+        // Calculate the differential taps using a derivative filter
+        std::vector<float> diff_filter(2);
+        diff_filter[0] = -1;
+        diff_filter[1] = 1;
+
+        difftaps.push_back(0);
         for(unsigned int i = 0; i < newtaps.size()-1; i++) {
-          tap = newtaps[i+1] - newtaps[i];
+          float tap = 0;
+          for(unsigned int j = 0; j < diff_filter.size(); j++) {
+            tap += diff_filter[j]*newtaps[i+j];
+          }
           difftaps.push_back(tap);
         }
-        difftaps.push_back(tap);
+        difftaps.push_back(0);
       }
 
       void
@@ -182,6 +205,13 @@ namespace gr {
         return d_taps_per_filter;
       }
 
+      float
+      pfb_arb_resampler_ccf::phase_offset(float freq, float fs)
+      {
+        float adj = (2.0*M_PI)*(freq/fs)/static_cast<float>(d_int_rate);
+        return -adj * d_est_phase_change;
+      }
+
       int
       pfb_arb_resampler_ccf::filter(gr_complex *output, gr_complex *input,
                                     int n_to_read, int &n_read)
@@ -213,6 +243,218 @@ namespace gr {
         n_read = i_in;   // return how much we've actually read
         return i_out;    // return how much we've produced
       }
+    
+      /****************************************************************/
+
+      pfb_arb_resampler_fff::pfb_arb_resampler_fff(float rate,
+                                                   const std::vector<float> &taps,
+                                                   unsigned int filter_size)
+      {
+        d_acc = 0; // start accumulator at 0
+
+        /* The number of filters is specified by the user as the
+           filter size; this is also the interpolation rate of the
+           filter. We use it and the rate provided to determine the
+           decimation rate. This acts as a rational resampler. The
+           flt_rate is calculated as the residual between the integer
+           decimation rate and the real decimation rate and will be
+           used to determine to interpolation point of the resampling
+           process.
+        */
+        d_int_rate = filter_size;
+        set_rate(rate);
+
+        d_last_filter = (taps.size()/2) % filter_size;
+        
+        d_filters = std::vector<fir_filter_fff*>(d_int_rate);
+        d_diff_filters = std::vector<fir_filter_fff*>(d_int_rate);
+
+        // Create an FIR filter for each channel and zero out the taps
+        std::vector<float> vtaps(0, d_int_rate);
+        for(unsigned int i = 0; i < d_int_rate; i++) {
+          d_filters[i] = new fir_filter_fff(1, vtaps);
+          d_diff_filters[i] = new fir_filter_fff(1, vtaps);
+        }
+
+        // Now, actually set the filters' taps
+        set_taps(taps);
+
+        // Delay is based on number of taps per filter arm. Round to
+        // the nearest integer.
+        float delay = -rate * (taps_per_filter() - 1.0) / 2.0;
+        d_delay = static_cast<int>(boost::math::iround(delay));
+
+        // This calculation finds the phase offset induced by the
+        // arbitrary resampling. It's based on which filter arm we are
+        // at the filter's group delay plus the fractional offset
+        // between the samples. Calculated here based on the rotation
+        // around nfilts starting at start_filter.
+        float accum = -d_delay * d_flt_rate;
+        int   accum_int = static_cast<int>(accum);
+        float accum_frac = accum - accum_int;
+        int end_filter = static_cast<int>
+          (boost::math::iround(fmodf(d_last_filter - d_delay * d_dec_rate + accum_int, \
+                                     static_cast<float>(d_int_rate))));
+
+        d_est_phase_change = d_last_filter - (end_filter + accum_frac);
+      }
+
+      pfb_arb_resampler_fff::~pfb_arb_resampler_fff()
+      {
+        for(unsigned int i = 0; i < d_int_rate; i++) {
+          delete d_filters[i];
+          delete d_diff_filters[i];
+        }
+      }
+
+      void
+      pfb_arb_resampler_fff::create_taps(const std::vector<float> &newtaps,
+                                         std::vector< std::vector<float> > &ourtaps,
+                                         std::vector<fir_filter_fff*> &ourfilter)
+      {
+        unsigned int ntaps = newtaps.size();
+        d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate);
+
+        // Create d_numchan vectors to store each channel's taps
+        ourtaps.resize(d_int_rate);
+
+        // Make a vector of the taps plus fill it out with 0's to fill
+        // each polyphase filter with exactly d_taps_per_filter
+        std::vector<float> tmp_taps;
+        tmp_taps = newtaps;
+        while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
+          tmp_taps.push_back(0.0);
+        }
+
+        for(unsigned int i = 0; i < d_int_rate; i++) {
+          // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
+          ourtaps[i] = std::vector<float>(d_taps_per_filter, 0);
+          for(unsigned int j = 0; j < d_taps_per_filter; j++) {
+            ourtaps[i][j] = tmp_taps[i + j*d_int_rate];
+          }
+
+          // Build a filter for each channel and add it's taps to it
+          ourfilter[i]->set_taps(ourtaps[i]);
+        }
+      }
+
+      void
+      pfb_arb_resampler_fff::create_diff_taps(const std::vector<float> &newtaps,
+                                              std::vector<float> &difftaps)
+      {
+        // Calculate the differential taps using a derivative filter
+        std::vector<float> diff_filter(2);
+        diff_filter[0] = -1;
+        diff_filter[1] = 1;
+
+        difftaps.push_back(0);
+        for(unsigned int i = 0; i < newtaps.size()-1; i++) {
+          float tap = 0;
+          for(unsigned int j = 0; j < diff_filter.size(); j++) {
+            tap += diff_filter[j]*newtaps[i+j];
+          }
+          difftaps.push_back(tap);
+        }
+        difftaps.push_back(0);
+      }
+
+      void
+      pfb_arb_resampler_fff::set_taps(const std::vector<float> &taps)
+      {
+        std::vector<float> dtaps;
+        create_diff_taps(taps, dtaps);
+        create_taps(taps, d_taps, d_filters);
+        create_taps(dtaps, d_dtaps, d_diff_filters);
+      }
+ 
+      std::vector<std::vector<float> >
+      pfb_arb_resampler_fff::taps() const
+      {
+        return d_taps;
+      }
+
+      void
+      pfb_arb_resampler_fff::print_taps()
+      {
+        unsigned int i, j;
+        for(i = 0; i < d_int_rate; i++) {
+          printf("filter[%d]: [", i);
+          for(j = 0; j < d_taps_per_filter; j++) {
+            printf(" %.4e", d_taps[i][j]);
+          }
+          printf("]\n");
+        }
+      }
+
+      void
+      pfb_arb_resampler_fff::set_rate(float rate)
+      {
+        d_dec_rate = (unsigned int)floor(d_int_rate/rate);
+        d_flt_rate = (d_int_rate/rate) - d_dec_rate;
+      }
+
+      void
+      pfb_arb_resampler_fff::set_phase(float ph)
+      {
+        if((ph < 0) || (ph >= 2.0*M_PI)) {
+          throw std::runtime_error("pfb_arb_resampler_fff: set_phase value out of bounds [0, 2pi).\n");
+        }
+        
+        float ph_diff = 2.0*M_PI / (float)d_filters.size();
+        d_last_filter = static_cast<int>(ph / ph_diff);
+      }
+
+      float
+      pfb_arb_resampler_fff::phase() const
+      {
+        float ph_diff = 2.0*M_PI / static_cast<float>(d_filters.size());
+        return d_last_filter * ph_diff;
+      }
+
+      unsigned int
+      pfb_arb_resampler_fff::taps_per_filter() const
+      {
+        return d_taps_per_filter;
+      }
+
+      float
+      pfb_arb_resampler_fff::phase_offset(float freq, float fs)
+      {
+        float adj = (2.0*M_PI)*(freq/fs)/static_cast<float>(d_int_rate);
+        return -adj * d_est_phase_change;
+      }
+
+      int
+      pfb_arb_resampler_fff::filter(float *output, float *input,
+                                    int n_to_read, int &n_read)
+      {
+        int i_out = 0, i_in = 0;
+        unsigned int j = d_last_filter;;
+        float o0, o1;
+        
+        while(i_in < n_to_read) {
+          // start j by wrapping around mod the number of channels
+          while(j < d_int_rate) {
+            // Take the current filter and derivative filter output
+            o0 = d_filters[j]->filter(&input[i_in]);
+            o1 = d_diff_filters[j]->filter(&input[i_in]);
+
+            output[i_out] = o0 + o1*d_acc;     // linearly interpolate between samples
+            i_out++;
+
+            // Adjust accumulator and index into filterbank
+            d_acc += d_flt_rate;
+            j += d_dec_rate + (int)floor(d_acc);
+            d_acc = fmodf(d_acc, 1.0);
+          }
+          i_in += (int)(j / d_int_rate);
+          j = j % d_int_rate;
+        }
+        d_last_filter = j; // save last filter state for re-entry
+
+        n_read = i_in;   // return how much we've actually read
+        return i_out;    // return how much we've produced
+      }
 
     } /* namespace kernel */
   } /* namespace filter */
diff --git a/gr-filter/lib/pfb_arb_resampler_ccf_impl.cc b/gr-filter/lib/pfb_arb_resampler_ccf_impl.cc
index a34df89e87..23950f36eb 100644
--- a/gr-filter/lib/pfb_arb_resampler_ccf_impl.cc
+++ b/gr-filter/lib/pfb_arb_resampler_ccf_impl.cc
@@ -48,7 +48,6 @@ namespace gr {
               io_signature::make(1, 1, sizeof(gr_complex)),
               io_signature::make(1, 1, sizeof(gr_complex)))
     {
-      d_start_index = 0;
       d_updated = false;
 
       d_resamp = new kernel::pfb_arb_resampler_ccf(rate, taps, filter_size);
@@ -66,8 +65,14 @@ namespace gr {
                                          gr_vector_int &ninput_items_required)
     {
       unsigned ninputs = ninput_items_required.size();
-      for(unsigned i = 0; i < ninputs; i++)
-        ninput_items_required[i] = noutput_items/relative_rate() + history() - 1;
+      if(noutput_items / relative_rate() < 1) {
+        for(unsigned i = 0; i < ninputs; i++)
+          ninput_items_required[i] = max_output_buffer(i)-1;
+      }
+      else {
+        for(unsigned i = 0; i < ninputs; i++)
+          ninput_items_required[i] = noutput_items/relative_rate() + history() - 1;
+      }
     }
 
     void
@@ -76,7 +81,7 @@ namespace gr {
       gr::thread::scoped_lock guard(d_mutex);
 
       d_resamp->set_taps(taps);
-      set_history(d_resamp->taps_per_filter() + 1);
+      set_history(d_resamp->taps_per_filter());
       d_updated = true;
     }
  
@@ -115,12 +120,42 @@ namespace gr {
     }
 
     unsigned int
+    pfb_arb_resampler_ccf_impl::interpolation_rate() const
+    { 
+      return d_resamp->interpolation_rate();
+    }
+
+    unsigned int
+    pfb_arb_resampler_ccf_impl::decimation_rate() const
+    {
+      return d_resamp->decimation_rate();
+    }
+
+    float
+    pfb_arb_resampler_ccf_impl::fractional_rate() const
+    {
+      return d_resamp->fractional_rate();
+    }
+
+    unsigned int
     pfb_arb_resampler_ccf_impl::taps_per_filter() const
     {
       return d_resamp->taps_per_filter();
     }
 
     int
+    pfb_arb_resampler_ccf_impl::group_delay() const
+    {
+      return d_resamp->group_delay();
+    }
+
+    float
+    pfb_arb_resampler_ccf_impl::phase_offset(float freq, float fs)
+    {
+      return d_resamp->phase_offset(freq, fs);
+    }
+
+    int
     pfb_arb_resampler_ccf_impl::general_work(int noutput_items,
 					     gr_vector_int &ninput_items,
 					     gr_vector_const_void_star &input_items,
diff --git a/gr-filter/lib/pfb_arb_resampler_ccf_impl.h b/gr-filter/lib/pfb_arb_resampler_ccf_impl.h
index 6d7a5a9619..aaafaf9e75 100644
--- a/gr-filter/lib/pfb_arb_resampler_ccf_impl.h
+++ b/gr-filter/lib/pfb_arb_resampler_ccf_impl.h
@@ -1,6 +1,6 @@
 /* -*- c++ -*- */
 /*
- * Copyright 2009,2010,2012 Free Software Foundation, Inc.
+ * Copyright 2009,2010,2013 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -35,8 +35,7 @@ namespace gr {
     {
     private:
       kernel::pfb_arb_resampler_ccf *d_resamp;
-      int          d_start_index;
-      bool         d_updated;
+      bool d_updated;
       gr::thread::mutex d_mutex; // mutex to protect set/work access
 
     public:
@@ -55,8 +54,15 @@ namespace gr {
       void set_rate(float rate);
       void set_phase(float ph);
       float phase() const;
+
+      unsigned int interpolation_rate() const;
+      unsigned int decimation_rate() const;
+      float fractional_rate() const;
       unsigned int taps_per_filter() const;
 
+      int group_delay() const;
+      float phase_offset(float freq, float fs);
+
       int general_work(int noutput_items,
 		       gr_vector_int &ninput_items,
 		       gr_vector_const_void_star &input_items,
diff --git a/gr-filter/lib/pfb_arb_resampler_fff_impl.cc b/gr-filter/lib/pfb_arb_resampler_fff_impl.cc
index c17f637776..55985d0d1f 100644
--- a/gr-filter/lib/pfb_arb_resampler_fff_impl.cc
+++ b/gr-filter/lib/pfb_arb_resampler_fff_impl.cc
@@ -46,127 +46,55 @@ namespace gr {
 							   unsigned int filter_size)
       : block("pfb_arb_resampler_fff",
 		 io_signature::make(1, 1, sizeof(float)),
-		 io_signature::make(1, 1, sizeof(float))),
-	d_updated(false)
+		 io_signature::make(1, 1, sizeof(float)))
     {
-      d_acc = 0; // start accumulator at 0
+      d_updated = false;
 
-      /* The number of filters is specified by the user as the filter
-	 size; this is also the interpolation rate of the filter. We
-	 use it and the rate provided to determine the decimation
-	 rate. This acts as a rational resampler. The flt_rate is
-	 calculated as the residual between the integer decimation
-	 rate and the real decimation rate and will be used to
-	 determine to interpolation point of the resampling process.
-      */
-      d_int_rate = filter_size;
-      set_rate(rate);
-
-      // Store the last filter between calls to work
-      d_last_filter = 0;
-
-      d_start_index = 0;
-
-      d_filters = std::vector<kernel::fir_filter_fff*>(d_int_rate);
-      d_diff_filters = std::vector<kernel::fir_filter_fff*>(d_int_rate);
-
-      // Create an FIR filter for each channel and zero out the taps
-      std::vector<float> vtaps(0, d_int_rate);
-      for(unsigned int i = 0; i < d_int_rate; i++) {
-	d_filters[i] = new kernel::fir_filter_fff(1, vtaps);
-	d_diff_filters[i] = new kernel::fir_filter_fff(1, vtaps);
-      }
-
-      // Now, actually set the filters' taps
-      std::vector<float> dtaps;
-      create_diff_taps(taps, dtaps);
-      create_taps(taps, d_taps, d_filters);
-      create_taps(dtaps, d_dtaps, d_diff_filters);
+      d_resamp = new kernel::pfb_arb_resampler_fff(rate, taps, filter_size);
+      set_history(d_resamp->taps_per_filter());
+      set_relative_rate(rate);
     }
 
     pfb_arb_resampler_fff_impl::~pfb_arb_resampler_fff_impl()
     {
-      for(unsigned int i = 0; i < d_int_rate; i++) {
-	delete d_filters[i];
-	delete d_diff_filters[i];
-      }
+      delete d_resamp;
     }
 
     void
-    pfb_arb_resampler_fff_impl::create_taps(const std::vector<float> &newtaps,
-					    std::vector< std::vector<float> > &ourtaps,
-					    std::vector<kernel::fir_filter_fff*> &ourfilter)
+    pfb_arb_resampler_fff_impl::forecast(int noutput_items,
+                                         gr_vector_int &ninput_items_required)
     {
-      unsigned int ntaps = newtaps.size();
-      d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate);
-
-      // Create d_numchan vectors to store each channel's taps
-      ourtaps.resize(d_int_rate);
-
-      // Make a vector of the taps plus fill it out with 0's to fill
-      // each polyphase filter with exactly d_taps_per_filter
-      std::vector<float> tmp_taps;
-      tmp_taps = newtaps;
-      while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
-	tmp_taps.push_back(0.0);
+      unsigned ninputs = ninput_items_required.size();
+      if(noutput_items / relative_rate() < 1) {
+        for(unsigned i = 0; i < ninputs; i++)
+          ninput_items_required[i] = max_output_buffer(i)-1;
       }
-
-      // Partition the filter
-      for(unsigned int i = 0; i < d_int_rate; i++) {
-	// Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
-	ourtaps[d_int_rate-1-i] = std::vector<float>(d_taps_per_filter, 0);
-	for(unsigned int j = 0; j < d_taps_per_filter; j++) {
-	  ourtaps[d_int_rate - 1 - i][j] = tmp_taps[i + j*d_int_rate];
-	}
-
-	// Build a filter for each channel and add it's taps to it
-	ourfilter[i]->set_taps(ourtaps[d_int_rate-1-i]);
+      else {
+        for(unsigned i = 0; i < ninputs; i++)
+          ninput_items_required[i] = noutput_items/relative_rate() + history() - 1;
       }
-
-      // Set the history to ensure enough input items for each filter
-      set_history(d_taps_per_filter + 1);
-
-      d_updated = true;
-    }
-
-    void
-    pfb_arb_resampler_fff_impl::create_diff_taps(const std::vector<float> &newtaps,
-						 std::vector<float> &difftaps)
-    {
-      // Calculate the differential taps (derivative filter) by taking the difference
-      // between two taps. Duplicate the last one to make both filters the same length.
-      float tap;
-      difftaps.clear();
-      for(unsigned int i = 0; i < newtaps.size()-1; i++) {
-	tap = newtaps[i+1] - newtaps[i];
-	difftaps.push_back(tap);
-      }
-      difftaps.push_back(tap);
     }
 
     void
     pfb_arb_resampler_fff_impl::set_taps(const std::vector<float> &taps)
     {
       gr::thread::scoped_lock guard(d_mutex);
+
+      d_resamp->set_taps(taps);
+      set_history(d_resamp->taps_per_filter());
+      d_updated = true;
     }
  
     std::vector<std::vector<float> >
     pfb_arb_resampler_fff_impl::taps() const
     {
-      return d_taps;
+      return d_resamp->taps();
     }
 
     void
     pfb_arb_resampler_fff_impl::print_taps()
     {
-      unsigned int i, j;
-      for(i = 0; i < d_int_rate; i++) {
-	printf("filter[%d]: [", i);
-	for(j = 0; j < d_taps_per_filter; j++) {
-	  printf(" %.4e", d_taps[i][j]);
-	}
-	printf("]\n");
-      }
+      d_resamp->print_taps();
     }
 
     void
@@ -174,8 +102,7 @@ namespace gr {
     {
       gr::thread::scoped_lock guard(d_mutex);
 
-      d_dec_rate = (unsigned int)floor(d_int_rate/rate);
-      d_flt_rate = (d_int_rate/rate) - d_dec_rate;
+      d_resamp->set_rate(rate);
       set_relative_rate(rate);
     }
 
@@ -183,19 +110,49 @@ namespace gr {
     pfb_arb_resampler_fff_impl::set_phase(float ph)
     {
       gr::thread::scoped_lock guard(d_mutex);
-      if((ph < 0) || (ph >= 2.0*M_PI)) {
-	throw std::runtime_error("pfb_arb_resampler_ccf: set_phase value out of bounds [0, 2pi).\n");
-      }
-      
-      float ph_diff = 2.0*M_PI / (float)d_filters.size();
-      d_last_filter = static_cast<int>(ph / ph_diff);
+      d_resamp->set_phase(ph);
     }
 
     float
     pfb_arb_resampler_fff_impl::phase() const
     {
-      float ph_diff = 2.0*M_PI / static_cast<float>(d_filters.size());
-      return d_last_filter * ph_diff;
+      return d_resamp->phase();
+    }
+
+    unsigned int
+    pfb_arb_resampler_fff_impl::interpolation_rate() const
+    { 
+      return d_resamp->interpolation_rate();
+    }
+
+    unsigned int
+    pfb_arb_resampler_fff_impl::decimation_rate() const
+    {
+      return d_resamp->decimation_rate();
+    }
+
+    float
+    pfb_arb_resampler_fff_impl::fractional_rate() const
+    {
+      return d_resamp->fractional_rate();
+    }
+
+    unsigned int
+    pfb_arb_resampler_fff_impl::taps_per_filter() const
+    {
+      return d_resamp->taps_per_filter();
+    }
+
+    int
+    pfb_arb_resampler_fff_impl::group_delay() const
+    {
+      return d_resamp->group_delay();
+    }
+
+    float
+    pfb_arb_resampler_fff_impl::phase_offset(float freq, float fs)
+    {
+      return d_resamp->phase_offset(freq, fs);
     }
 
     int
@@ -214,44 +171,12 @@ namespace gr {
 	return 0;		     // history requirements may have changed.
       }
 
-      int i = 0, count = d_start_index;
-      unsigned int j;
-      float o0, o1;
-
-      // Restore the last filter position
-      j = d_last_filter;
-
-      // produce output as long as we can and there are enough input samples
-      int max_input = ninput_items[0] - (int)d_taps_per_filter;
-      while((i < noutput_items) && (count < max_input)) {
-	// start j by wrapping around mod the number of channels
-	while((j < d_int_rate) && (i < noutput_items)) {
-	  // Take the current filter and derivative filter output
-	  o0 = d_filters[j]->filter(&in[count]);
-	  o1 = d_diff_filters[j]->filter(&in[count]);
-
-	  out[i] = o0 + o1*d_acc;     // linearly interpolate between samples
-	  i++;
-
-	  // Adjust accumulator and index into filterbank
-	  d_acc += d_flt_rate;
-	  j += d_dec_rate + (int)floor(d_acc);
-	  d_acc = fmodf(d_acc, 1.0);
-	}
-	if(i < noutput_items) {              // keep state for next entry
-	  float ss = (int)(j / d_int_rate);  // number of items to skip ahead by
-	  count += ss;                       // we have fully consumed another input
-	  j = j % d_int_rate;                // roll filter around
-	}
-      }
-
-      // Store the current filter position and start of next sample
-      d_last_filter = j;
-      d_start_index = std::max(0, count - ninput_items[0]);
+      int nitems_read;
+      int nitems = floorf((float)noutput_items / relative_rate());
+      int processed = d_resamp->filter(out, in, nitems, nitems_read);
 
-      // consume all we've processed but no more than we can
-      consume_each(std::min(count, ninput_items[0]));
-      return i;
+      consume_each(nitems_read);
+      return processed;
     }
 
   } /* namespace filter */
diff --git a/gr-filter/lib/pfb_arb_resampler_fff_impl.h b/gr-filter/lib/pfb_arb_resampler_fff_impl.h
index e6022f1c18..69cc439caa 100644
--- a/gr-filter/lib/pfb_arb_resampler_fff_impl.h
+++ b/gr-filter/lib/pfb_arb_resampler_fff_impl.h
@@ -1,6 +1,6 @@
 /* -*- c++ -*- */
 /*
- * Copyright 2009-2012 Free Software Foundation, Inc.
+ * Copyright 2009-2013 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -25,7 +25,7 @@
 #define	INCLUDED_PFB_ARB_RESAMPLER_FFF_IMPL_H
 
 #include <gnuradio/filter/pfb_arb_resampler_fff.h>
-#include <gnuradio/filter/fir_filter.h>
+#include <gnuradio/filter/pfb_arb_resampler.h>
 #include <gnuradio/thread/thread.h>
 
 namespace gr {
@@ -34,33 +34,10 @@ namespace gr {
     class FILTER_API pfb_arb_resampler_fff_impl : public pfb_arb_resampler_fff
     {
     private:
-      std::vector<kernel::fir_filter_fff*> d_filters;
-      std::vector<kernel::fir_filter_fff*> d_diff_filters;
-      std::vector< std::vector<float> > d_taps;
-      std::vector< std::vector<float> > d_dtaps;
-      unsigned int d_int_rate;          // the number of filters (interpolation rate)
-      unsigned int d_dec_rate;          // the stride through the filters (decimation rate)
-      float        d_flt_rate;          // residual rate for the linear interpolation
-      float        d_acc;
-      unsigned int d_last_filter;
-      int          d_start_index;
-      unsigned int d_taps_per_filter;
-      bool         d_updated;
+      kernel::pfb_arb_resampler_fff *d_resamp;
+      bool d_updated;
       gr::thread::mutex d_mutex; // mutex to protect set/work access
 
-      void create_diff_taps(const std::vector<float> &newtaps,
-			    std::vector<float> &difftaps);
-
-      /*!
-       * Resets the filterbank's filter taps with the new prototype filter
-       * \param newtaps    (vector of floats) The prototype filter to populate the filterbank.
-       *                   The taps should be generated at the interpolated sampling rate.
-       * \param ourtaps    (vector of floats) Reference to our internal member of holding the taps.
-       * \param ourfilter  (vector of filters) Reference to our internal filter to set the taps for.
-       */
-      void create_taps(const std::vector<float> &newtaps,
-		       std::vector< std::vector<float> > &ourtaps,
-		       std::vector<kernel::fir_filter_fff*> &ourfilter);
     public:
       pfb_arb_resampler_fff_impl(float rate,
 				 const std::vector<float> &taps,
@@ -68,14 +45,24 @@ namespace gr {
 
       ~pfb_arb_resampler_fff_impl();
 
+      void forecast(int noutput_items, gr_vector_int &ninput_items_required);
+
       void set_taps(const std::vector<float> &taps);
       std::vector<std::vector<float> > taps() const;
       void print_taps();
-      void set_rate(float rate);
 
+      void set_rate(float rate);
       void set_phase(float ph);
       float phase() const;
 
+      unsigned int interpolation_rate() const;
+      unsigned int decimation_rate() const;
+      float fractional_rate() const;
+      unsigned int taps_per_filter() const;
+
+      int group_delay() const;
+      float phase_offset(float freq, float fs);
+
       int general_work(int noutput_items,
 		       gr_vector_int &ninput_items,
 		       gr_vector_const_void_star &input_items,
diff --git a/gr-filter/python/filter/qa_pfb_arb_resampler.py b/gr-filter/python/filter/qa_pfb_arb_resampler.py
index 15827be357..674d508e45 100755
--- a/gr-filter/python/filter/qa_pfb_arb_resampler.py
+++ b/gr-filter/python/filter/qa_pfb_arb_resampler.py
@@ -21,9 +21,7 @@
 #
 
 from gnuradio import gr, gr_unittest, filter, blocks
-
 import math
-from scipy import signal
 
 def sig_source_c(samp_rate, freq, amp, N):
     t = map(lambda x: float(x)/samp_rate, xrange(N))
@@ -45,19 +43,19 @@ class test_pfb_arb_resampler(gr_unittest.TestCase):
         self.tb = None
 
     def test_fff_000(self):
-        N = 1000         # number of samples to use
-        fs = 1000        # baseband sampling rate
-        rrate = 1.123    # resampling rate
+        N = 500         # number of samples to use
+        fs = 5000.0     # baseband sampling rate
+        rrate = 2.3421  # resampling rate
 
         nfilts = 32
         taps = filter.firdes.low_pass_2(nfilts, nfilts*fs, fs/2, fs/10,
                                         attenuation_dB=80,
                                         window=filter.firdes.WIN_BLACKMAN_hARRIS)
 
-        freq = 100
+        freq = 121.213
         data = sig_source_f(fs, freq, 1, N)
         signal = blocks.vector_source_f(data)
-        pfb = filter.pfb_arb_resampler_fff(rrate, taps)
+        pfb = filter.pfb_arb_resampler_fff(rrate, taps, nfilts)
         snk = blocks.vector_sink_f()
 
         self.tb.connect(signal, pfb, snk)
@@ -65,77 +63,57 @@ class test_pfb_arb_resampler(gr_unittest.TestCase):
 
         Ntest = 50
         L = len(snk.data())
-        t = map(lambda x: float(x)/(fs*rrate), xrange(L))
 
-        phase = 0.53013
+        # Get group delay and estimate of phase offset from the filter itself.
+        delay = pfb.group_delay()
+        phase = pfb.phase_offset(freq, fs)
+        
+        # Create a timeline offset by the filter's group delay
+        t = map(lambda x: float(x)/(fs*rrate), xrange(delay, L+delay))
+
+        # Data of the sinusoid at frequency freq with the delay and phase offset.
         expected_data = map(lambda x: math.sin(2.*math.pi*freq*x+phase), t)
 
         dst_data = snk.data()
-        self.assertFloatTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 3)
+
+        self.assertFloatTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 2)
 
     def test_ccf_000(self):
-        N = 1000         # number of samples to use
-        fs = 1000        # baseband sampling rate
-        rrate = 1.123    # resampling rate
+        N = 5000         # number of samples to use
+        fs = 5000.0      # baseband sampling rate
+        rrate = 2.4321   # resampling rate
 
         nfilts = 32
         taps = filter.firdes.low_pass_2(nfilts, nfilts*fs, fs/2, fs/10,
                                         attenuation_dB=80,
                                         window=filter.firdes.WIN_BLACKMAN_hARRIS)
 
-        freq = 100
+        freq = 211.123
         data = sig_source_c(fs, freq, 1, N)
         signal = blocks.vector_source_c(data)
-        pfb = filter.pfb_arb_resampler_ccf(rrate, taps)
+        pfb = filter.pfb_arb_resampler_ccf(rrate, taps, nfilts)
         snk = blocks.vector_sink_c()
-
+        
         self.tb.connect(signal, pfb, snk)
         self.tb.run() 
 
         Ntest = 50
         L = len(snk.data())
-        t = map(lambda x: float(x)/(fs*rrate), xrange(L))
 
-        phase = 0.53013
-        delay = ((len(taps)/nfilts) - 1.0) / 2.0
-        print "DELAY: ", delay
-        print delay/(fs*rrate)
+        # Get group delay and estimate of phase offset from the filter itself.
+        delay = pfb.group_delay()
+        phase = pfb.phase_offset(freq, fs)
+        
+        # Create a timeline offset by the filter's group delay
+        t = map(lambda x: float(x)/(fs*rrate), xrange(delay, L+delay))
+
+        # Data of the sinusoid at frequency freq with the delay and phase offset.
         expected_data = map(lambda x: math.cos(2.*math.pi*freq*x+phase) + \
                                 1j*math.sin(2.*math.pi*freq*x+phase), t)
 
         dst_data = snk.data()
-        self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 3)
 
-#    def test_ccf_001(self):
-#        N = 1000         # number of samples to use
-#        fs = 1000        # baseband sampling rate
-#        rrate = 1.123    # resampling rate
-#
-#        nfilts = 32
-#        taps = filter.firdes.low_pass_2(nfilts, nfilts*fs, fs/2, fs/10,
-#                                        attenuation_dB=80,
-#                                        window=filter.firdes.WIN_BLACKMAN_hARRIS)
-#
-#        freq = 100
-#        data = sig_source_c(fs, freq, 1, N)
-#        src = blocks.vector_source_c(data)
-#        pfb = filter.pfb_arb_resampler_ccf(rrate, taps)
-#        snk = blocks.vector_sink_c()
-#
-#        self.tb.connect(src, pfb, snk)
-#        self.tb.run() 
-#
-#        Ntest = 50
-#        L = len(snk.data())
-#        t = map(lambda x: float(x)/(fs*rrate), xrange(L))
-#
-#        resamp = len(data)*rrate
-#        expected_data = signal.resample(data, resamp)
-#
-#        dst_data = snk.data()
-#        print expected_data[-Ntest:]
-#        print dst_data[-Ntest:]
-#        self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 3)
+        self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 2)
 
 if __name__ == '__main__':
     gr_unittest.run(test_pfb_arb_resampler, "test_pfb_arb_resampler.xml")