diff options
author | Tom Rondeau <trondeau@vt.edu> | 2013-06-13 13:32:22 -0400 |
---|---|---|
committer | Tom Rondeau <trondeau@vt.edu> | 2013-06-13 15:27:26 -0400 |
commit | 47ab3e90cd9d2a5a1932e56f10c794724ad0ce1d (patch) | |
tree | 29cd469fef63c0f16fd2430468f5ad8617b85767 | |
parent | 4c7192af4fb76de36bcce6a428d33ec8fbe1bd0e (diff) |
filter: Updated PFB arbitrary resamplers (both ccf and fff) to use new kernel, created more accessors, and have them now calculating their own group delay and expected phase offset.
QA code for resampler uses resamplers' delay and phase offset to set up expected offsets.
-rw-r--r-- | gr-filter/include/gnuradio/filter/pfb_arb_resampler.h | 209 | ||||
-rw-r--r-- | gr-filter/include/gnuradio/filter/pfb_arb_resampler_ccf.h | 27 | ||||
-rw-r--r-- | gr-filter/include/gnuradio/filter/pfb_arb_resampler_fff.h | 36 | ||||
-rw-r--r-- | gr-filter/lib/pfb_arb_resampler.cc | 268 | ||||
-rw-r--r-- | gr-filter/lib/pfb_arb_resampler_ccf_impl.cc | 43 | ||||
-rw-r--r-- | gr-filter/lib/pfb_arb_resampler_ccf_impl.h | 12 | ||||
-rw-r--r-- | gr-filter/lib/pfb_arb_resampler_fff_impl.cc | 205 | ||||
-rw-r--r-- | gr-filter/lib/pfb_arb_resampler_fff_impl.h | 43 | ||||
-rwxr-xr-x | gr-filter/python/filter/qa_pfb_arb_resampler.py | 82 |
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") |