[Commit-gnuradio] r11577 - gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter

[trondeau]

Author: trondeau
Date: 2009-08-11 17:22:20 -0600 (Tue, 11 Aug 2009)
New Revision: 11577

Modified:
   
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
   
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
   
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
Log:
Bug fixes and documentation for the PFB arbitrary resampler. This was the most 
problematic but now has nice supression in the images.

Modified: 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
===================================================================
--- 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
  2009-08-11 22:56:56 UTC (rev 11576)
+++ 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
  2009-08-11 23:22:20 UTC (rev 11577)
@@ -30,23 +30,38 @@
 #include <gr_io_signature.h>
 
 gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate, 
-                                                            const 
std::vector<float> &taps)
+                                                            const 
std::vector<float> &taps,
+                                                            unsigned int 
filter_size)
 {
-  return gr_pfb_arb_resampler_ccf_sptr (new gr_pfb_arb_resampler_ccf (rate, 
taps));
+  return gr_pfb_arb_resampler_ccf_sptr (new gr_pfb_arb_resampler_ccf (rate, 
taps,
+                                                                     
filter_size));
 }
 
 
 gr_pfb_arb_resampler_ccf::gr_pfb_arb_resampler_ccf (float rate, 
-                                                   const std::vector<float> 
&taps)
+                                                   const std::vector<float> 
&taps,
+                                                   unsigned int filter_size)
   : gr_block ("pfb_arb_resampler_ccf",
              gr_make_io_signature (1, 1, sizeof(gr_complex)),
              gr_make_io_signature (1, 1, sizeof(gr_complex))),
     d_updated (false)
 {
-  d_int_rate = 32;
+  /* 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;
   d_dec_rate = (unsigned int)floor(d_int_rate/rate);
   d_flt_rate = (d_int_rate/rate) - d_dec_rate;
+
+  // The accumulator keeps track of overflow to increment the stride correctly.
   d_acc = 0;
+
+  // Store the last filter between calls to work
   d_last_filter = 0;
   
   d_filters = std::vector<gr_fir_ccf*>(d_int_rate);
@@ -78,13 +93,21 @@
 
   // Create d_numchan vectors to store each channel's taps
   d_taps.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 = taps;
+  while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
+    tmp_taps.push_back(0.0);
+  }
   
   // Partition the filter
   for(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
     d_taps[i] = std::vector<float>(d_taps_per_filter, 0);
     for(j = 0; j < d_taps_per_filter; j++) {
-      d_taps[i][j] = taps[i + j*d_int_rate];  // add taps to channels in 
reverse order
+      d_taps[i][j] = tmp_taps[i + j*d_int_rate];  // add taps to channels in 
reverse order
     }
     
     // Build a filter for each channel and add it's taps to it
@@ -124,34 +147,48 @@
     return 0;               // history requirements may have changed.
   }
 
-  int i = 0, j = 0, count = 0;
+  int i = 0, j, count = 0;
   gr_complex o0, o1;
 
+  // Restore the last filter position
   j = d_last_filter;
-  while(i < noutput_items) {
+
+  // produce output as long as we can and there are enough input samples
+  while((i < noutput_items) && (count < ninput_items[0]-1)) {
+
     // start j by wrapping around mod the number of channels
-    
-    while(j < d_int_rate) {
+    while((j < d_int_rate) && (i < noutput_items)) {
       // Take the current filter output
       o0 = d_filters[j]->filter(&in[count]);
 
       // Take the next filter output; wrap around to 0 if necessary
       if(j+1 == d_int_rate)
+       // Use the sample of the next input item through the first filter
        o1 = d_filters[0]->filter(&in[count+1]);
-      else
-       o1 = d_filters[j]->filter(&in[count]);
+      else {
+       // Use the sample from the current input item through the nex filter
+       o1 = d_filters[j+1]->filter(&in[count]);
+      }
 
-      out[i] = o0 + (o1 - o0)*d_flt_rate;
+      //out[i] = o0;                     // nearest-neighbor approach
+      out[i] = o0 + (o1 - o0)*d_acc;     // linearly interpolate between 
samples
       i++;
-
-      d_acc += d_dec_rate + d_flt_rate;
-      j = (int)floor(d_acc);
+      
+      // Accumulate the position in the stream for the interpolated point.
+      // If it goes above 1, roll around to zero and increment the stride
+      // length this time by the decimation rate plus 1 for the increment
+      // due to the acculated position.
+      d_acc += d_flt_rate;
+      j += d_dec_rate + (int)floor(d_acc);
+      d_acc = fmodf(d_acc, 1.0);
     }
-    count++;
-    j = j % d_int_rate;
-    d_acc = d_acc - floor(d_acc);    
+    if(i < noutput_items) {              // keep state for next entry
+      count++;                           // we have fully consumed another 
input
+      j = j % d_int_rate;                // roll filter around
+    }
   }
 
+  // Store the current filter position
   d_last_filter = j;
 
   consume_each(count);

Modified: 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
===================================================================
--- 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
   2009-08-11 22:56:56 UTC (rev 11576)
+++ 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
   2009-08-11 23:22:20 UTC (rev 11577)
@@ -29,41 +29,124 @@
 class gr_pfb_arb_resampler_ccf;
 typedef boost::shared_ptr<gr_pfb_arb_resampler_ccf> 
gr_pfb_arb_resampler_ccf_sptr;
 gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate,
-                                                            const 
std::vector<float> &taps);
+                                                            const 
std::vector<float> &taps,
+                                                            unsigned int 
filter_size=32);
 
 class gr_fir_ccf;
 
 /*!
- * \brief FIR filter with gr_complex input, gr_complex output and float taps
+ * \class gr_pfb_arb_resampler_ccf
+ *
+ * \brief Polyphase filterbank arbitrary resampler with 
+ *        gr_complex input, gr_complex output and float taps
+ *
  * \ingroup filter_blk
+ * 
+ * This block 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 = gr.firdes.low_pass_2(32, 32*fs, BW, TB, 
+ *           attenuation_dB=ATT, window=gr.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 gr_pfb_arb_resampler_ccf : public gr_block
 {
  private:
+  /*!
+   * Build the polyphase filterbank arbitray resampler.
+   * \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.
+   */
   friend gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float 
rate,
-                                                                     const 
std::vector<float> &taps);
+                                                                     const 
std::vector<float> &taps,
+                                                                     unsigned 
int filter_size);
 
   std::vector<gr_fir_ccf*> d_filters;
   std::vector< std::vector<float> > d_taps;
-  unsigned int             d_int_rate;
-  unsigned int             d_dec_rate;
-  float                    d_flt_rate;
+  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;
   unsigned int             d_taps_per_filter;
   bool                    d_updated;
 
   /*!
-   * Construct a Polyphase filterbank for channelization with the given 
-   * number of channels and taps
+   * Build the polyphase filterbank arbitray resampler.
+   * \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.
    */
   gr_pfb_arb_resampler_ccf (float rate, 
-                          const std::vector<float> &taps);
+                           const std::vector<float> &taps,
+                           unsigned int filter_size);
   
 public:
   ~gr_pfb_arb_resampler_ccf ();
   
+  /*!
+   * Resets the filterbank's filter taps with the new prototype filter
+   * \param taps    (vector/list of floats) The prototype filter to populate 
the filterbank. The taps
+   *                                        should be generated at the 
interpolated sampling rate.
+   */
   void set_taps (const std::vector<float> &taps);
+
+  /*!
+   * Print all of the filterbank taps to screen.
+   */
   void print_taps();
   
   int general_work (int noutput_items,

Modified: 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
===================================================================
--- 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
   2009-08-11 22:56:56 UTC (rev 11576)
+++ 
gnuradio/branches/developers/trondeau/pfb/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
   2009-08-11 23:22:20 UTC (rev 11577)
@@ -23,13 +23,15 @@
 GR_SWIG_BLOCK_MAGIC(gr,pfb_arb_resampler_ccf);
 
 gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate,
-                                                            const 
std::vector<float> &taps);
+                                                            const 
std::vector<float> &taps,
+                                                            unsigned int 
filter_size=32);
 
 class gr_pfb_arb_resampler_ccf : public gr_block
 {
  private:
   gr_pfb_arb_resampler_ccf (float rate,
-                           const std::vector<float> &taps);
+                           const std::vector<float> &taps,
+                           unsigned int filter_size);
 
  public:
   ~gr_pfb_arb_resampler_ccf ();