git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2009 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 #ifdef HAVE_CONFIG_H
  24 #include "config.h"
  25 #endif
  26
  27 #include <gr_pfb_arb_resampler_ccf.h>
  28 #include <gr_fir_ccf.h>
  29 #include <gr_fir_util.h>
  30 #include <gr_io_signature.h>
  31
  32 gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate,
  33                                                              const std::vector<float> &taps,
  34                                                              unsigned int filter_size)
  35 {
  36   return gr_pfb_arb_resampler_ccf_sptr (new gr_pfb_arb_resampler_ccf (rate, taps,
  37                                                                       filter_size));
  38 }
  39
  40
  41 gr_pfb_arb_resampler_ccf::gr_pfb_arb_resampler_ccf (float rate,
  42                                                     const std::vector<float> &taps,
  43                                                     unsigned int filter_size)
  44   : gr_block ("pfb_arb_resampler_ccf",
  45               gr_make_io_signature (1, 1, sizeof(gr_complex)),
  46               gr_make_io_signature (1, 1, sizeof(gr_complex))),
  47     d_updated (false)
  48 {
  49   /* The number of filters is specified by the user as the filter size;
  50      this is also the interpolation rate of the filter. We use it and the
  51      rate provided to determine the decimation rate. This acts as a
  52      rational resampler. The flt_rate is calculated as the residual
  53      between the integer decimation rate and the real decimation rate and
  54      will be used to determine to interpolation point of the resampling
  55      process.
  56   */
  57   d_int_rate = filter_size;
  58   d_dec_rate = (unsigned int)floor(d_int_rate/rate);
  59   d_flt_rate = (d_int_rate/rate) - d_dec_rate;
  60
  61   // The accumulator keeps track of overflow to increment the stride correctly.
  62   d_acc = 0;
  63
  64   // Store the last filter between calls to work
  65   d_last_filter = 0;
  66
  67   d_filters = std::vector<gr_fir_ccf*>(d_int_rate);
  68
  69   // Create an FIR filter for each channel and zero out the taps
  70   std::vector<float> vtaps(0, d_int_rate);
  71   for(unsigned int i = 0; i < d_int_rate; i++) {
  72     d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps);
  73   }
  74
  75   // Now, actually set the filters' taps
  76   set_taps(taps);
  77 }
  78
  79 gr_pfb_arb_resampler_ccf::~gr_pfb_arb_resampler_ccf ()
  80 {
  81   for(unsigned int i = 0; i < d_int_rate; i++) {
  82     delete d_filters[i];
  83   }
  84 }
  85
  86 void
  87 gr_pfb_arb_resampler_ccf::set_taps (const std::vector<float> &taps)
  88 {
  89   unsigned int i,j;
  90
  91   unsigned int ntaps = taps.size();
  92   d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate);
  93
  94   // Create d_numchan vectors to store each channel's taps
  95   d_taps.resize(d_int_rate);
  96
  97   // Make a vector of the taps plus fill it out with 0's to fill
  98   // each polyphase filter with exactly d_taps_per_filter
  99   std::vector<float> tmp_taps;
 100   tmp_taps = taps;
 101   while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
 102     tmp_taps.push_back(0.0);
 103   }
 104
 105   // Partition the filter
 106   for(i = 0; i < d_int_rate; i++) {
 107     // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
 108     d_taps[i] = std::vector<float>(d_taps_per_filter, 0);
 109     for(j = 0; j < d_taps_per_filter; j++) {
 110       d_taps[i][j] = tmp_taps[i + j*d_int_rate];  // add taps to channels in reverse order
 111     }
 112
 113     // Build a filter for each channel and add it's taps to it
 114     d_filters[i]->set_taps(d_taps[i]);
 115   }
 116
 117   // Set the history to ensure enough input items for each filter
 118   set_history (d_taps_per_filter);
 119
 120   d_updated = true;
 121 }
 122
 123 void
 124 gr_pfb_arb_resampler_ccf::print_taps()
 125 {
 126   unsigned int i, j;
 127   for(i = 0; i < d_int_rate; i++) {
 128     printf("filter[%d]: [", i);
 129     for(j = 0; j < d_taps_per_filter; j++) {
 130       printf(" %.4e", d_taps[i][j]);
 131     }
 132     printf("]\n");
 133   }
 134 }
 135
 136 int
 137 gr_pfb_arb_resampler_ccf::general_work (int noutput_items,
 138                                         gr_vector_int &ninput_items,
 139                                         gr_vector_const_void_star &input_items,
 140                                         gr_vector_void_star &output_items)
 141 {
 142   gr_complex *in = (gr_complex *) input_items[0];
 143   gr_complex *out = (gr_complex *) output_items[0];
 144
 145   if (d_updated) {
 146     d_updated = false;
 147     return 0;                // history requirements may have changed.
 148   }
 149
 150   int i = 0, j, count = 0;
 151   gr_complex o0, o1;
 152
 153   // Restore the last filter position
 154   j = d_last_filter;
 155
 156   // produce output as long as we can and there are enough input samples
 157   while((i < noutput_items) && (count < ninput_items[0]-1)) {
 158
 159     // start j by wrapping around mod the number of channels
 160     while((j < d_int_rate) && (i < noutput_items)) {
 161       // Take the current filter output
 162       o0 = d_filters[j]->filter(&in[count]);
 163
 164       // Take the next filter output; wrap around to 0 if necessary
 165       if(j+1 == d_int_rate)
 166         // Use the sample of the next input item through the first filter
 167         o1 = d_filters[0]->filter(&in[count+1]);
 168       else {
 169         // Use the sample from the current input item through the nex filter
 170         o1 = d_filters[j+1]->filter(&in[count]);
 171       }
 172
 173       //out[i] = o0;                     // nearest-neighbor approach
 174       out[i] = o0 + (o1 - o0)*d_acc;     // linearly interpolate between samples
 175       i++;
 176
 177       // Accumulate the position in the stream for the interpolated point.
 178       // If it goes above 1, roll around to zero and increment the stride
 179       // length this time by the decimation rate plus 1 for the increment
 180       // due to the acculated position.
 181       d_acc += d_flt_rate;
 182       j += d_dec_rate + (int)floor(d_acc);
 183       d_acc = fmodf(d_acc, 1.0);
 184     }
 185     if(i < noutput_items) {              // keep state for next entry
 186       count++;                           // we have fully consumed another input
 187       j = j % d_int_rate;                // roll filter around
 188     }
 189   }
 190
 191   // Store the current filter position
 192   d_last_filter = j;
 193
 194   consume_each(count);
 195   return i;
 196 }