From a0d13b42bfb3fd081d77e9d73cf4db9695a6d88b Mon Sep 17 00:00:00 2001 From: trondeau Date: Wed, 12 Aug 2009 03:39:03 +0000 Subject: [PATCH] Merging trondeau/pfb r11249:11581 into trunk. This adds a few polyphase filterbank implementations that do (integer) decimation, (integer) interpolation, arbitrary resampling, and channelizing. gnuradio-example/python/pfb includes a number of different examples of how to use these blocks. git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@11583 221aa14e-8319-0410-a670-987f0aec2ac5 --- config/grc_gnuradio_examples.m4 | 1 + gnuradio-core/src/lib/filter/Makefile.am | 16 +- gnuradio-core/src/lib/filter/filter.i | 13 +- .../lib/filter/gr_pfb_arb_resampler_ccf.cc | 196 +++++++++++++++ .../src/lib/filter/gr_pfb_arb_resampler_ccf.h | 158 ++++++++++++ .../src/lib/filter/gr_pfb_arb_resampler_ccf.i | 41 ++++ .../src/lib/filter/gr_pfb_channelizer_ccf.cc | 150 ++++++++++++ .../src/lib/filter/gr_pfb_channelizer_ccf.h | 144 +++++++++++ .../src/lib/filter/gr_pfb_channelizer_ccf.i | 38 +++ .../src/lib/filter/gr_pfb_decimator_ccf.cc | 174 ++++++++++++++ .../src/lib/filter/gr_pfb_decimator_ccf.h | 148 ++++++++++++ .../src/lib/filter/gr_pfb_decimator_ccf.i | 41 ++++ .../src/lib/filter/gr_pfb_interpolator_ccf.cc | 142 +++++++++++ .../src/lib/filter/gr_pfb_interpolator_ccf.h | 129 ++++++++++ .../src/lib/filter/gr_pfb_interpolator_ccf.i | 39 +++ .../src/python/gnuradio/blks2impl/Makefile.am | 4 + .../gnuradio/blks2impl/pfb_arb_resampler.py | 50 ++++ .../gnuradio/blks2impl/pfb_channelizer.py | 57 +++++ .../gnuradio/blks2impl/pfb_decimator.py | 49 ++++ .../gnuradio/blks2impl/pfb_interpolator.py | 49 ++++ gnuradio-examples/python/Makefile.am | 1 + gnuradio-examples/python/pfb/Makefile.am | 31 +++ gnuradio-examples/python/pfb/channelize.py | 177 ++++++++++++++ .../python/pfb/chirp_channelize.py | 192 +++++++++++++++ gnuradio-examples/python/pfb/decimate.py | 171 +++++++++++++ gnuradio-examples/python/pfb/fmtest.py | 197 +++++++++++++++ gnuradio-examples/python/pfb/interpolate.py | 226 ++++++++++++++++++ gr-utils/src/python/gr_plot_const.py | 8 +- gr-utils/src/python/gr_plot_fft.py | 13 +- gr-utils/src/python/gr_plot_iq.py | 7 +- gr-utils/src/python/gr_plot_psd.py | 13 +- gr-utils/src/python/plot_data.py | 3 - 32 files changed, 2642 insertions(+), 36 deletions(-) create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.i create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.i create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h create mode 100644 gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.i create mode 100644 gnuradio-core/src/python/gnuradio/blks2impl/pfb_arb_resampler.py create mode 100644 gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py create mode 100644 gnuradio-core/src/python/gnuradio/blks2impl/pfb_decimator.py create mode 100644 gnuradio-core/src/python/gnuradio/blks2impl/pfb_interpolator.py create mode 100644 gnuradio-examples/python/pfb/Makefile.am create mode 100755 gnuradio-examples/python/pfb/channelize.py create mode 100755 gnuradio-examples/python/pfb/chirp_channelize.py create mode 100755 gnuradio-examples/python/pfb/decimate.py create mode 100755 gnuradio-examples/python/pfb/fmtest.py create mode 100755 gnuradio-examples/python/pfb/interpolate.py diff --git a/config/grc_gnuradio_examples.m4 b/config/grc_gnuradio_examples.m4 index 7de55666..3225f3ab 100644 --- a/config/grc_gnuradio_examples.m4 +++ b/config/grc_gnuradio_examples.m4 @@ -39,6 +39,7 @@ AC_DEFUN([GRC_GNURADIO_EXAMPLES],[ gnuradio-examples/python/multi_usrp/Makefile \ gnuradio-examples/python/network/Makefile \ gnuradio-examples/python/ofdm/Makefile \ + gnuradio-examples/python/pfb/Makefile \ gnuradio-examples/python/usrp/Makefile \ gnuradio-examples/python/usrp2/Makefile \ ]) diff --git a/gnuradio-core/src/lib/filter/Makefile.am b/gnuradio-core/src/lib/filter/Makefile.am index e230e88b..838f69b9 100644 --- a/gnuradio-core/src/lib/filter/Makefile.am +++ b/gnuradio-core/src/lib/filter/Makefile.am @@ -201,7 +201,11 @@ libfilter_la_common_SOURCES = \ complex_dotprod_generic.cc \ ccomplex_dotprod_generic.cc \ float_dotprod_generic.c \ - short_dotprod_generic.c + short_dotprod_generic.c \ + gr_pfb_channelizer_ccf.cc \ + gr_pfb_decimator_ccf.cc \ + gr_pfb_interpolator_ccf.cc \ + gr_pfb_arb_resampler_ccf.cc libfilter_qa_la_common_SOURCES = \ qa_filter.cc \ @@ -276,7 +280,11 @@ grinclude_HEADERS = \ qa_filter.h \ short_dotprod_generic.h \ short_dotprod_x86.h \ - sse_debug.h + sse_debug.h \ + gr_pfb_channelizer_ccf.h \ + gr_pfb_decimator_ccf.h \ + gr_pfb_interpolator_ccf.h \ + gr_pfb_arb_resampler_ccf.h noinst_HEADERS = \ assembly.h \ @@ -327,6 +335,10 @@ swiginclude_HEADERS = \ gr_iir_filter_ffd.i \ gr_single_pole_iir_filter_ff.i \ gr_single_pole_iir_filter_cc.i \ + gr_pfb_channelizer_ccf.i \ + gr_pfb_decimator_ccf.i \ + gr_pfb_interpolator_ccf.i \ + gr_pfb_arb_resampler_ccf.i \ $(GENERATED_I) endif diff --git a/gnuradio-core/src/lib/filter/filter.i b/gnuradio-core/src/lib/filter/filter.i index 7931a4d5..16b8005b 100644 --- a/gnuradio-core/src/lib/filter/filter.i +++ b/gnuradio-core/src/lib/filter/filter.i @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2004,2005,2006,2007 Free Software Foundation, Inc. + * Copyright 2004,2005,2006,2007,2009 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -32,6 +32,10 @@ #include #include #include +#include +#include +#include +#include %} %include "gr_iir_filter_ffd.i" @@ -45,5 +49,12 @@ %include "gr_fractional_interpolator_cc.i" %include "gr_goertzel_fc.i" %include "gr_cma_equalizer_cc.i" +%include "gr_pfb_channelizer_ccf.i" +%include "gr_pfb_decimator_ccf.i" +%include "gr_pfb_interpolator_ccf.i" +%include "gr_pfb_arb_resampler_ccf.i" +%include "gr_pfb_decimator_ccf.i" +%include "gr_pfb_interpolator_ccf.i" +%include "gr_pfb_arb_resampler_ccf.i" %include "filter_generated.i" diff --git a/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc new file mode 100644 index 00000000..bfc4c046 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc @@ -0,0 +1,196 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include + +gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate, + const std::vector &taps, + unsigned int filter_size) +{ + 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 &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) +{ + /* 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(d_int_rate); + + // Create an FIR filter for each channel and zero out the taps + std::vector vtaps(0, d_int_rate); + for(unsigned int i = 0; i < d_int_rate; i++) { + d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps); + } + + // Now, actually set the filters' taps + set_taps(taps); +} + +gr_pfb_arb_resampler_ccf::~gr_pfb_arb_resampler_ccf () +{ + for(unsigned int i = 0; i < d_int_rate; i++) { + delete d_filters[i]; + } +} + +void +gr_pfb_arb_resampler_ccf::set_taps (const std::vector &taps) +{ + unsigned int i,j; + + unsigned int ntaps = taps.size(); + d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate); + + // 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 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(d_taps_per_filter, 0); + for(j = 0; j < d_taps_per_filter; j++) { + 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 + d_filters[i]->set_taps(d_taps[i]); + } + + // Set the history to ensure enough input items for each filter + set_history (d_taps_per_filter); + + d_updated = true; +} + +void +gr_pfb_arb_resampler_ccf::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"); + } +} + +int +gr_pfb_arb_resampler_ccf::general_work (int noutput_items, + gr_vector_int &ninput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + gr_complex *in = (gr_complex *) input_items[0]; + gr_complex *out = (gr_complex *) output_items[0]; + + if (d_updated) { + d_updated = false; + return 0; // history requirements may have changed. + } + + int i = 0, j, count = 0; + gr_complex o0, o1; + + // Restore the last filter position + j = d_last_filter; + + // 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) && (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 { + // Use the sample from the current input item through the nex filter + o1 = d_filters[j+1]->filter(&in[count]); + } + + //out[i] = o0; // nearest-neighbor approach + out[i] = o0 + (o1 - o0)*d_acc; // linearly interpolate between samples + i++; + + // 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); + } + 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); + return i; +} diff --git a/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h new file mode 100644 index 00000000..b79a89fe --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h @@ -0,0 +1,158 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + + +#ifndef INCLUDED_GR_PFB_ARB_RESAMPLER_CCF_H +#define INCLUDED_GR_PFB_ARB_RESAMPLER_CCF_H + +#include + +class gr_pfb_arb_resampler_ccf; +typedef boost::shared_ptr gr_pfb_arb_resampler_ccf_sptr; +gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate, + const std::vector &taps, + unsigned int filter_size=32); + +class gr_fir_ccf; + +/*! + * \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 r. The resampling is done by constructing + * N filters where N is the interpolation rate. We + * then calculate D where D = floor(N/r). + * + * Using N and D, we can perform rational resampling + * where N/D is a rational number close to the input rate + * r where we have N filters and we cycle through + * them as a polyphase filterbank with a stride of D so that + * i+1 = (i + D) % N. + * + * To get the arbitrary rate, we want to interpolate between two + * points. For each value out, we take an output from the current + * filter, i, and the next filter i+1 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, N, used determines the + * quantization error; the larger N, the smaller the + * noise. You can design for a specified noise floor by setting the + * filter size (parameters filter_size). 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 fs, a 3-dB bandwidth of + * BW and a transition bandwidth of TB. We can also + * specify the out-of-band attenuation to use, ATT, 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 (32). + * + * self._taps = gr.firdes.low_pass_2(32, 32*fs, BW, TB, + * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS) + * + * The theory behind this block can be found in Chapter 7.5 of + * the following book. + * + * f. harris, Multirate Signal Processing for Communication + * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004. + */ + +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 &taps, + unsigned int filter_size); + + std::vector d_filters; + std::vector< std::vector > d_taps; + 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; + + /*! + * 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 &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 &taps); + + /*! + * Print all of the filterbank taps to screen. + */ + void print_taps(); + + int general_work (int noutput_items, + gr_vector_int &ninput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); +}; + +#endif diff --git a/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i new file mode 100644 index 00000000..e365e031 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i @@ -0,0 +1,41 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +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 &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 &taps, + unsigned int filter_size); + + public: + ~gr_pfb_arb_resampler_ccf (); + + void set_taps (const std::vector &taps); + void print_taps(); +}; diff --git a/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc new file mode 100644 index 00000000..7be611e2 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc @@ -0,0 +1,150 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include +#include + +gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, + const std::vector &taps) +{ + return gr_pfb_channelizer_ccf_sptr (new gr_pfb_channelizer_ccf (numchans, taps)); +} + + +gr_pfb_channelizer_ccf::gr_pfb_channelizer_ccf (unsigned int numchans, + const std::vector &taps) + : gr_sync_block ("pfb_channelizer_ccf", + gr_make_io_signature (numchans, numchans, sizeof(gr_complex)), + gr_make_io_signature (1, 1, numchans*sizeof(gr_complex))), + d_updated (false) +{ + d_numchans = numchans; + d_filters = std::vector(d_numchans); + + // Create an FIR filter for each channel and zero out the taps + std::vector vtaps(0, d_numchans); + for(unsigned int i = 0; i < d_numchans; i++) { + d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps); + } + + // Now, actually set the filters' taps + set_taps(taps); + + // Create the FFT to handle the output de-spinning of the channels + d_fft = new gri_fft_complex (d_numchans, false); +} + +gr_pfb_channelizer_ccf::~gr_pfb_channelizer_ccf () +{ + for(unsigned int i = 0; i < d_numchans; i++) { + delete d_filters[i]; + } +} + +void +gr_pfb_channelizer_ccf::set_taps (const std::vector &taps) +{ + unsigned int i,j; + + unsigned int ntaps = taps.size(); + d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_numchans); + + // Create d_numchan vectors to store each channel's taps + d_taps.resize(d_numchans); + + // 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 tmp_taps; + tmp_taps = taps; + while((float)(tmp_taps.size()) < d_numchans*d_taps_per_filter) { + tmp_taps.push_back(0.0); + } + + // Partition the filter + for(i = 0; i < d_numchans; i++) { + // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out + d_taps[i] = std::vector(d_taps_per_filter, 0); + for(j = 0; j < d_taps_per_filter; j++) { + d_taps[i][j] = tmp_taps[i + j*d_numchans]; // add taps to channels in reverse order + } + + // Build a filter for each channel and add it's taps to it + d_filters[i]->set_taps(d_taps[i]); + } + + // Set the history to ensure enough input items for each filter + set_history (d_taps_per_filter); + + d_updated = true; +} + +void +gr_pfb_channelizer_ccf::print_taps() +{ + unsigned int i, j; + for(i = 0; i < d_numchans; i++) { + printf("filter[%d]: [", i); + for(j = 0; j < d_taps_per_filter; j++) { + printf(" %.4e", d_taps[i][j]); + } + printf("]\n\n"); + } +} + + +int +gr_pfb_channelizer_ccf::work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + gr_complex *in = (gr_complex *) input_items[0]; + gr_complex *out = (gr_complex *) output_items[0]; + + if (d_updated) { + d_updated = false; + return 0; // history requirements may have changed. + } + + for(int i = 0; i < noutput_items; i++) { + // Move through filters from bottom to top + for(int j = d_numchans-1; j >= 0; j--) { + // Take in the items from the first input stream to d_numchans + in = (gr_complex*)input_items[d_numchans - 1 - j]; + + // Filter current input stream from bottom filter to top + d_fft->get_inbuf()[j] = d_filters[j]->filter(&in[i]); + } + + // despin through FFT + d_fft->execute(); + memcpy(&out[d_numchans*i], d_fft->get_outbuf(), d_numchans*sizeof(gr_complex)); + } + + return noutput_items; +} diff --git a/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h new file mode 100644 index 00000000..7d0a31c5 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h @@ -0,0 +1,144 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + + +#ifndef INCLUDED_GR_PFB_CHANNELIZER_CCF_H +#define INCLUDED_GR_PFB_CHANNELIZER_CCF_H + +#include + +class gr_pfb_channelizer_ccf; +typedef boost::shared_ptr gr_pfb_channelizer_ccf_sptr; +gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, + const std::vector &taps); + +class gr_fir_ccf; +class gri_fft_complex; + + +/*! + * \class gr_pfb_channelizer_ccf + * + * \brief Polyphase filterbank channelizer with + * gr_complex input, gr_complex output and float taps + * + * \ingroup filter_blk + * + * This block takes in complex inputs and channelizes it to M + * channels of equal bandwidth. Each of the resulting channels is + * decimated to the new rate that is the input sampling rate + * fs divided by the number of channels, M. + * + * The PFB channelizer code takes the taps generated above and builds + * a set of filters. The set contains M number of filters + * and each filter contains ceil(taps.size()/decim) number of taps. + * Each tap from the filter prototype is sequentially inserted into + * the next filter. When all of the input taps are used, the remaining + * filters in the filterbank are filled out with 0's to make sure each + * filter has the same number of taps. + * + * Each filter operates using the gr_fir filter classs of GNU Radio, + * which takes the input stream at i and performs the inner + * product calculation to i+(n-1) where n is the + * number of filter taps. To efficiently handle this in the GNU Radio + * structure, each filter input must come from its own input + * stream. So the channelizer must be provided with M streams + * where the input stream has been deinterleaved. This is most easily + * done using the gr_stream_to_streams block. + * + * The output is then produced as a vector, where index i in + * the vector is the next sample from the ith channel. This + * is most easily handled by sending the output to a + * gr_vector_to_streams block to handle the conversion and passing + * M streams out. + * + * The input and output formatting is done using a hier_block2 called + * pfb_channelizer_ccf. This can take in a single stream and outputs + * M streams based on the behavior described above. + * + * The filter's taps should be based on the input sampling rate. + * + * For example, using the GNU Radio's firdes utility to building + * filters, we build a low-pass filter with a sampling rate of + * fs, a 3-dB bandwidth of BW and a transition + * bandwidth of TB. We can also specify the out-of-band + * attenuation to use, ATT, 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 unity. + * + * self._taps = gr.firdes.low_pass_2(1, fs, BW, TB, + * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS) + * + * The theory behind this block can be found in Chapter 6 of + * the following book. + * + * f. harris, Multirate Signal Processing for Communication + * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004. + * + */ + +class gr_pfb_channelizer_ccf : public gr_sync_block +{ + private: + /*! + * Build the polyphase filterbank decimator. + * \param numchans (unsigned integer) Specifies the number of channels M + * \param taps (vector/list of floats) The prototype filter to populate the filterbank. + */ + friend gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, + const std::vector &taps); + + std::vector d_filters; + std::vector< std::vector > d_taps; + gri_fft_complex *d_fft; + unsigned int d_numchans; + unsigned int d_taps_per_filter; + bool d_updated; + + /*! + * Build the polyphase filterbank decimator. + * \param numchans (unsigned integer) Specifies the number of channels M + * \param taps (vector/list of floats) The prototype filter to populate the filterbank. + */ + gr_pfb_channelizer_ccf (unsigned int numchans, + const std::vector &taps); + +public: + ~gr_pfb_channelizer_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. + */ + void set_taps (const std::vector &taps); + + /*! + * Print all of the filterbank taps to screen. + */ + void print_taps(); + + int work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); +}; + +#endif diff --git a/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.i b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.i new file mode 100644 index 00000000..4bef90e2 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.i @@ -0,0 +1,38 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +GR_SWIG_BLOCK_MAGIC(gr,pfb_channelizer_ccf); + +gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, + const std::vector &taps); + +class gr_pfb_channelizer_ccf : public gr_sync_block +{ + private: + gr_pfb_channelizer_ccf (unsigned int numchans, + const std::vector &taps); + + public: + ~gr_pfb_channelizer_ccf (); + + void set_taps (const std::vector &taps); +}; diff --git a/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc b/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc new file mode 100644 index 00000000..b334f587 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc @@ -0,0 +1,174 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include +#include +#include + +gr_pfb_decimator_ccf_sptr gr_make_pfb_decimator_ccf (unsigned int decim, + const std::vector &taps, + unsigned int channel) +{ + return gr_pfb_decimator_ccf_sptr (new gr_pfb_decimator_ccf (decim, taps, channel)); +} + + +gr_pfb_decimator_ccf::gr_pfb_decimator_ccf (unsigned int decim, + const std::vector &taps, + unsigned int channel) + : gr_sync_block ("pfb_decimator_ccf", + gr_make_io_signature (decim, decim, sizeof(gr_complex)), + gr_make_io_signature (1, 1, sizeof(gr_complex))), + d_updated (false) +{ + d_rate = decim; + d_filters = std::vector(d_rate); + d_chan = channel; + d_rotator = new gr_complex[d_rate]; + + // Create an FIR filter for each channel and zero out the taps + std::vector vtaps(0, d_rate); + for(unsigned int i = 0; i < d_rate; i++) { + d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps); + d_rotator[i] = gr_expj(i*2*M_PI*d_chan/d_rate); + } + + // Now, actually set the filters' taps + set_taps(taps); + + // Create the FFT to handle the output de-spinning of the channels + d_fft = new gri_fft_complex (d_rate, false); +} + +gr_pfb_decimator_ccf::~gr_pfb_decimator_ccf () +{ + for(unsigned int i = 0; i < d_rate; i++) { + delete d_filters[i]; + } +} + +void +gr_pfb_decimator_ccf::set_taps (const std::vector &taps) +{ + unsigned int i,j; + + unsigned int ntaps = taps.size(); + d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_rate); + + // Create d_numchan vectors to store each channel's taps + d_taps.resize(d_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 tmp_taps; + tmp_taps = taps; + while((float)(tmp_taps.size()) < d_rate*d_taps_per_filter) { + tmp_taps.push_back(0.0); + } + + // Partition the filter + for(i = 0; i < d_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(d_taps_per_filter, 0); + for(j = 0; j < d_taps_per_filter; j++) { + d_taps[i][j] = tmp_taps[i + j*d_rate]; // add taps to channels in reverse order + } + + // Build a filter for each channel and add it's taps to it + d_filters[i]->set_taps(d_taps[i]); + } + + // Set the history to ensure enough input items for each filter + set_history (d_taps_per_filter); + + d_updated = true; +} + +void +gr_pfb_decimator_ccf::print_taps() +{ + unsigned int i, j; + for(i = 0; i < d_rate; i++) { + printf("filter[%d]: [", i); + for(j = 0; j < d_taps_per_filter; j++) { + printf(" %.4e", d_taps[i][j]); + } + printf("]\n\n"); + } +} + +#define ROTATEFFT + +int +gr_pfb_decimator_ccf::work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + gr_complex *in; + gr_complex *out = (gr_complex *) output_items[0]; + + if (d_updated) { + d_updated = false; + return 0; // history requirements may have changed. + } + + int i; + for(i = 0; i < noutput_items; i++) { + // Move through filters from bottom to top + out[i] = 0; + for(int j = d_rate-1; j >= 0; j--) { + // Take in the items from the first input stream to d_rate + in = (gr_complex*)input_items[d_rate - 1 - j]; + + // Filter current input stream from bottom filter to top + // The rotate them by expj(j*k*2pi/M) where M is the number of filters + // (the decimation rate) and k is the channel number to extract + + // This is the real math that goes on; we abuse the FFT to do this quickly + // for decimation rates > N where N is a small number (~5): + // out[i] += d_filters[j]->filter(&in[i])*gr_expj(j*d_chan*2*M_PI/d_rate); +#ifdef ROTATEFFT + d_fft->get_inbuf()[j] = d_filters[j]->filter(&in[i]); +#else + out[i] += d_filters[j]->filter(&in[i])*d_rotator[i]; +#endif + } + +#ifdef ROTATEFFT + // Perform the FFT to do the complex multiply despinning for all channels + d_fft->execute(); + + // Select only the desired channel out + out[i] = d_fft->get_outbuf()[d_chan]; +#endif + + } + + return noutput_items; +} diff --git a/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h b/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h new file mode 100644 index 00000000..83997c0c --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h @@ -0,0 +1,148 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + + +#ifndef INCLUDED_GR_PFB_DECIMATOR_CCF_H +#define INCLUDED_GR_PFB_DECIMATOR_CCF_H + +#include + +class gr_pfb_decimator_ccf; +typedef boost::shared_ptr gr_pfb_decimator_ccf_sptr; +gr_pfb_decimator_ccf_sptr gr_make_pfb_decimator_ccf (unsigned int decim, + const std::vector &taps, + unsigned int channel=0); + +class gr_fir_ccf; +class gri_fft_complex; + +/*! + * \class gr_pfb_decimator_ccf + * \brief Polyphase filterbank bandpass decimator with gr_complex + * input, gr_complex output and float taps + * + * \ingroup filter_blk + * + * This block takes in a signal stream and performs interger down- + * sampling (decimation) with a polyphase filterbank. The first input + * is the integer specifying how much to decimate by. The second + * input is a vector (Python list) of floating-point taps of the + * prototype filter. The third input specifies the channel to extract. + * By default, the zeroth channel is used, which is the baseband + * channel (first Nyquist zone). + * + * The channel parameter specifies which channel to use since + * this class is capable of bandpass decimation. Given a complex input + * stream at a sampling rate of fs and a decimation rate of + * decim, the input frequency domain is split into + * decim channels that represent the Nyquist zones. Using the + * polyphase filterbank, we can select any one of these channels to + * decimate. + * + * The output signal will be the basebanded and decimated signal from + * that channel. This concept is very similar to the PFB channelizer + * (see #gr_pfb_channelizer_ccf) where only a single channel is + * extracted at a time. + * + * The filter's taps should be based on the sampling rate before + * decimation. + * + * For example, using the GNU Radio's firdes utility to building + * filters, we build a low-pass filter with a sampling rate of + * fs, a 3-dB bandwidth of BW and a transition + * bandwidth of TB. We can also specify the out-of-band + * attenuation to use, ATT, 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 unity. + * + * self._taps = gr.firdes.low_pass_2(1, fs, BW, TB, + * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS) + * + * The PFB decimator code takes the taps generated above and builds a + * set of filters. The set contains decim number of filters + * and each filter contains ceil(taps.size()/decim) number of taps. + * Each tap from the filter prototype is sequentially inserted into + * the next filter. When all of the input taps are used, the remaining + * filters in the filterbank are filled out with 0's to make sure each + * filter has the same number of taps. + * + * The theory behind this block can be found in Chapter 6 of + * the following book. + * + * f. harris, Multirate Signal Processing for Communication + * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004. + */ + +class gr_pfb_decimator_ccf : public gr_sync_block +{ + private: + /*! + * Build the polyphase filterbank decimator. + * \param decim (unsigned integer) Specifies the decimation rate to use + * \param taps (vector/list of floats) The prototype filter to populate the filterbank. + * \param channel (unsigned integer) Selects the channel to return [default=0]. + */ + friend gr_pfb_decimator_ccf_sptr gr_make_pfb_decimator_ccf (unsigned int decim, + const std::vector &taps, + unsigned int channel); + + std::vector d_filters; + std::vector< std::vector > d_taps; + gri_fft_complex *d_fft; + unsigned int d_rate; + unsigned int d_chan; + unsigned int d_taps_per_filter; + bool d_updated; + gr_complex *d_rotator; + + /*! + * Build the polyphase filterbank decimator. + * \param decim (unsigned integer) Specifies the decimation rate to use + * \param taps (vector/list of floats) The prototype filter to populate the filterbank. + * \param channel (unsigned integer) Selects the channel to return [default=0]. + */ + gr_pfb_decimator_ccf (unsigned int decim, + const std::vector &taps, + unsigned int channel); + +public: + ~gr_pfb_decimator_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. + */ + void set_taps (const std::vector &taps); + + /*! + * Print all of the filterbank taps to screen. + */ + void print_taps(); + + //void set_channel (unsigned int channel); + + int work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); +}; + +#endif diff --git a/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.i b/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.i new file mode 100644 index 00000000..c4215fce --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.i @@ -0,0 +1,41 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +GR_SWIG_BLOCK_MAGIC(gr,pfb_decimator_ccf); + +gr_pfb_decimator_ccf_sptr gr_make_pfb_decimator_ccf (unsigned int decim, + const std::vector &taps, + unsigned int channel); + +class gr_pfb_decimator_ccf : public gr_sync_block +{ + private: + gr_pfb_decimator_ccf (unsigned int decim, + const std::vector &taps, + unsigned int channel); + + public: + ~gr_pfb_decimator_ccf (); + + void set_taps (const std::vector &taps); + //void set_channel (unsigned int channel); +}; diff --git a/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc b/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc new file mode 100644 index 00000000..d5eba885 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc @@ -0,0 +1,142 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include + +gr_pfb_interpolator_ccf_sptr gr_make_pfb_interpolator_ccf (unsigned int interp, + const std::vector &taps) +{ + return gr_pfb_interpolator_ccf_sptr (new gr_pfb_interpolator_ccf (interp, taps)); +} + + +gr_pfb_interpolator_ccf::gr_pfb_interpolator_ccf (unsigned int interp, + const std::vector &taps) + : gr_sync_interpolator ("pfb_interpolator_ccf", + gr_make_io_signature (1, 1, sizeof(gr_complex)), + gr_make_io_signature (1, 1, sizeof(gr_complex)), + interp), + d_updated (false) +{ + d_rate = interp; + d_filters = std::vector(d_rate); + + // Create an FIR filter for each channel and zero out the taps + std::vector vtaps(0, d_rate); + for(unsigned int i = 0; i < d_rate; i++) { + d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps); + } + + // Now, actually set the filters' taps + set_taps(taps); +} + +gr_pfb_interpolator_ccf::~gr_pfb_interpolator_ccf () +{ + for(unsigned int i = 0; i < d_rate; i++) { + delete d_filters[i]; + } +} + +void +gr_pfb_interpolator_ccf::set_taps (const std::vector &taps) +{ + unsigned int i,j; + + unsigned int ntaps = taps.size(); + d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_rate); + + // Create d_numchan vectors to store each channel's taps + //std::vector< std::vector > vtaps(d_rate); + d_taps.resize(d_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 tmp_taps; + tmp_taps = taps; + while((float)(tmp_taps.size()) < d_rate*d_taps_per_filter) { + tmp_taps.push_back(0.0); + } + + // Partition the filter + for(i = 0; i < d_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(d_taps_per_filter, 0); + for(j = 0; j < d_taps_per_filter; j++) { + d_taps[i][j] = tmp_taps[i + j*d_rate]; // add taps to channels in reverse order + } + + // Build a filter for each channel and add it's taps to it + d_filters[i]->set_taps(d_taps[i]); + } + + // Set the history to ensure enough input items for each filter + set_history (d_taps_per_filter); + + d_updated = true; +} + +void +gr_pfb_interpolator_ccf::print_taps() +{ + unsigned int i, j; + for(i = 0; i < d_rate; i++) { + printf("filter[%d]: [", i); + for(j = 0; j < d_taps_per_filter; j++) { + printf(" %.4e", d_taps[i][j]); + } + printf("]\n\n"); + } +} + +int +gr_pfb_interpolator_ccf::work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + gr_complex *in = (gr_complex *) input_items[0]; + gr_complex *out = (gr_complex *) output_items[0]; + + if (d_updated) { + d_updated = false; + return 0; // history requirements may have changed. + } + + int i = 0, count = 0; + + while(i < noutput_items) { + for(int j = 0; j < d_rate; j++) { + out[i] = d_filters[j]->filter(&in[count]); + i++; + } + count++; + } + + return i; +} diff --git a/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h b/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h new file mode 100644 index 00000000..50849d51 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h @@ -0,0 +1,129 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + + +#ifndef INCLUDED_GR_PFB_INTERPOLATOR_CCF_H +#define INCLUDED_GR_PFB_INTERPOLATOR_CCF_H + +#include + +class gr_pfb_interpolator_ccf; +typedef boost::shared_ptr gr_pfb_interpolator_ccf_sptr; +gr_pfb_interpolator_ccf_sptr gr_make_pfb_interpolator_ccf (unsigned int interp, + const std::vector &taps); + +class gr_fir_ccf; + +/*! + * \class gr_pfb_interpolator_ccf + * \brief Polyphase filterbank interpolator with gr_complex input, + * gr_complex output and float taps + * + * \ingroup filter_blk + * + * This block takes in a signal stream and performs interger up- + * sampling (interpolation) with a polyphase filterbank. The first + * input is the integer specifying how much to interpolate by. The + * second input is a vector (Python list) of floating-point taps of + * the prototype filter. + * + * The filter's taps should be based on the interpolation rate + * specified. That is, the bandwidth specified is relative to the + * bandwidth after interpolation. + * + * For example, using the GNU Radio's firdes utility to building + * filters, we build a low-pass filter with a sampling rate of + * fs, a 3-dB bandwidth of BW and a transition + * bandwidth of TB. We can also specify the out-of-band + * attenuation to use, ATT, and the filter window function (a + * Blackman-harris window in this case). The first input is the gain, + * which is also specified as the interpolation rate so that the + * output levels are the same as the input (this creates an overall + * increase in power). + * + * self._taps = gr.firdes.low_pass_2(interp, interp*fs, BW, TB, + * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS) + * + * The PFB interpolator code takes the taps generated above and builds + * a set of filters. The set contains interp number of + * filters and each filter contains ceil(taps.size()/interp) number of + * taps. Each tap from the filter prototype is sequentially inserted + * into the next filter. When all of the input taps are used, the + * remaining filters in the filterbank are filled out with 0's to make + * sure each filter has the same number of taps. + * + * The theory behind this block can be found in Chapter 7.1 of the + * following book. + * + * f. harris, Multirate Signal Processing for Communication + * Systems," Upper Saddle River, NJ: Prentice Hall, + * Inc. 2004. + */ + +class gr_pfb_interpolator_ccf : public gr_sync_interpolator +{ + private: + /*! + * Build the polyphase filterbank interpolator. + * \param interp (unsigned integer) Specifies the interpolation rate to use + * \param taps (vector/list of floats) The prototype filter to populate the filterbank. The taps + * should be generated at the interpolated sampling rate. + */ + friend gr_pfb_interpolator_ccf_sptr gr_make_pfb_interpolator_ccf (unsigned int interp, + const std::vector &taps); + + std::vector d_filters; + std::vector< std::vector > d_taps; + unsigned int d_rate; + unsigned int d_taps_per_filter; + bool d_updated; + + /*! + * Construct a Polyphase filterbank interpolator + * \param interp (unsigned integer) Specifies the interpolation rate to use + * \param taps (vector/list of floats) The prototype filter to populate the filterbank. The taps + * should be generated at the interpolated sampling rate. + */ + gr_pfb_interpolator_ccf (unsigned int interp, + const std::vector &taps); + +public: + ~gr_pfb_interpolator_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 &taps); + + /*! + * Print all of the filterbank taps to screen. + */ + void print_taps(); + + int work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); +}; + +#endif diff --git a/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.i b/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.i new file mode 100644 index 00000000..cf4302d4 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.i @@ -0,0 +1,39 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +GR_SWIG_BLOCK_MAGIC(gr,pfb_interpolator_ccf); + +gr_pfb_interpolator_ccf_sptr gr_make_pfb_interpolator_ccf (unsigned int interp, + const std::vector &taps); + +class gr_pfb_interpolator_ccf : public gr_sync_interpolator +{ + private: + gr_pfb_interpolator_ccf (unsigned int interp, + const std::vector &taps); + + public: + ~gr_pfb_interpolator_ccf (); + + void set_taps (const std::vector &taps); + void print_taps(); +}; diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am b/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am index f07abd4c..17be09cc 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am +++ b/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am @@ -47,6 +47,10 @@ grblkspython_PYTHON = \ ofdm_sync_pn.py \ ofdm_sync_pnac.py \ ofdm_sync_ml.py \ + pfb_arb_resampler.py \ + pfb_channelizer.py \ + pfb_decimator.py \ + pfb_interpolator.py \ pkt.py \ psk.py \ qam.py \ diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/pfb_arb_resampler.py b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_arb_resampler.py new file mode 100644 index 00000000..b1b3dfca --- /dev/null +++ b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_arb_resampler.py @@ -0,0 +1,50 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr + +class pfb_arb_resampler_ccf(gr.hier_block2): + ''' + Convinience wrapper for the polyphase filterbank arbitrary resampler. + + The block takes a single complex stream in and outputs a single complex + stream out. As such, it requires no extra glue to handle the input/output + streams. This block is provided to be consistent with the interface to the + other PFB block. + ''' + def __init__(self, rate, taps, flt_size=32): + gr.hier_block2.__init__(self, "pfb_arb_resampler_ccf", + gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature + gr.io_signature(1, 1, gr.sizeof_gr_complex)) # Output signature + + self._rate = rate + self._taps = taps + self._size = flt_size + + self.pfb = gr.pfb_arb_resampler_ccf(self._rate, self._taps, self._size) + + self.connect(self, self.pfb) + self.connect(self.pfb, self) + + + + diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py new file mode 100644 index 00000000..c45ae4d1 --- /dev/null +++ b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr + +class pfb_channelizer_ccf(gr.hier_block2): + ''' + Make a Polyphase Filter channelizer (complex in, complex out, floating-point taps) + + This simplifies the interface by allowing a single input stream to connect to this block. + It will then output a stream for each channel. + ''' + def __init__(self, numchans, taps): + gr.hier_block2.__init__(self, "pfb_channelizer_ccf", + gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature + gr.io_signature(numchans, numchans, gr.sizeof_gr_complex)) # Output signature + + self._numchans = numchans + self._taps = taps + + self.s2ss = gr.stream_to_streams(gr.sizeof_gr_complex, self._numchans) + self.pfb = gr.pfb_channelizer_ccf(self._numchans, self._taps) + self.v2s = gr.vector_to_streams(gr.sizeof_gr_complex, self._numchans) + + self.connect(self, self.s2ss) + + for i in xrange(self._numchans): + self.connect((self.s2ss,i), (self.pfb,i)) + + # Get independent streams from the filterbank and send them out + self.connect(self.pfb, self.v2s) + + for i in xrange(self._numchans): + self.connect((self.v2s,i), (self,i)) + + + + diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/pfb_decimator.py b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_decimator.py new file mode 100644 index 00000000..176d0473 --- /dev/null +++ b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_decimator.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr + +class pfb_decimator_ccf(gr.hier_block2): + ''' + Make a Polyphase Filter decimator (complex in, complex out, floating-point taps) + + This simplifies the interface by allowing a single input stream to connect to this block. + It will then output a stream that is the decimated output stream. + ''' + def __init__(self, decim, taps, channel=0): + gr.hier_block2.__init__(self, "pfb_decimator_ccf", + gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature + gr.io_signature(1, 1, gr.sizeof_gr_complex)) # Output signature + + self._decim = decim + self._taps = taps + self._channel = channel + + self.s2ss = gr.stream_to_streams(gr.sizeof_gr_complex, self._decim) + self.pfb = gr.pfb_decimator_ccf(self._decim, self._taps, self._channel) + + self.connect(self, self.s2ss) + + for i in xrange(self._decim): + self.connect((self.s2ss,i), (self.pfb,i)) + + self.connect(self.pfb, self) diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/pfb_interpolator.py b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_interpolator.py new file mode 100644 index 00000000..db294404 --- /dev/null +++ b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_interpolator.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr + +class pfb_interpolator_ccf(gr.hier_block2): + ''' + Make a Polyphase Filter interpolator (complex in, complex out, floating-point taps) + + The block takes a single complex stream in and outputs a single complex + stream out. As such, it requires no extra glue to handle the input/output + streams. This block is provided to be consistent with the interface to the + other PFB block. + ''' + def __init__(self, interp, taps): + gr.hier_block2.__init__(self, "pfb_interpolator_ccf", + gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature + gr.io_signature(1, 1, gr.sizeof_gr_complex)) # Output signature + + self._interp = interp + self._taps = taps + + self.pfb = gr.pfb_interpolator_ccf(self._interp, self._taps) + + self.connect(self, self.pfb) + self.connect(self.pfb, self) + + + + diff --git a/gnuradio-examples/python/Makefile.am b/gnuradio-examples/python/Makefile.am index 3a1acf51..ea03b438 100644 --- a/gnuradio-examples/python/Makefile.am +++ b/gnuradio-examples/python/Makefile.am @@ -32,5 +32,6 @@ SUBDIRS = \ multi_usrp \ network \ ofdm \ + pfb \ usrp \ usrp2 diff --git a/gnuradio-examples/python/pfb/Makefile.am b/gnuradio-examples/python/pfb/Makefile.am new file mode 100644 index 00000000..4aa9248e --- /dev/null +++ b/gnuradio-examples/python/pfb/Makefile.am @@ -0,0 +1,31 @@ +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +include $(top_srcdir)/Makefile.common + +ourdatadir = $(exampledir)/pfb + +dist_ourdata_SCRIPTS = \ + channelize.py \ + chirp_channelize.py \ + decimate.py \ + interpolate.py \ + fmtest.py diff --git a/gnuradio-examples/python/pfb/channelize.py b/gnuradio-examples/python/pfb/channelize.py new file mode 100755 index 00000000..bc83fae2 --- /dev/null +++ b/gnuradio-examples/python/pfb/channelize.py @@ -0,0 +1,177 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import os, time +import scipy, pylab +from scipy import fftpack +from pylab import mlab + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 2000000 # number of samples to use + self._fs = 9000 # initial sampling rate + self._M = 9 # Number of channels to channelize + + # Create a set of taps for the PFB channelizer + self._taps = gr.firdes.low_pass_2(1, self._fs, 475.50, 50, + attenuation_dB=10, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._M)) + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + # Create a set of signals at different frequencies + # freqs lists the frequencies of the signals that get stored + # in the list "signals", which then get summed together + self.signals = list() + self.add = gr.add_cc() + freqs = [-4070, -3050, -2030, -1010, 10, 1020, 2040, 3060, 4080] + for i in xrange(len(freqs)): + self.signals.append(gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freqs[i], 1)) + self.connect(self.signals[i], (self.add,i)) + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the channelizer filter + self.pfb = blks2.pfb_channelizer_ccf(self._M, self._taps) + + # Construct a vector sink for the input signal to the channelizer + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.add, self.head, self.pfb) + self.connect(self.add, self.snk_i) + + # Create a vector sink for each of M output channels of the filter and connect it + self.snks = list() + for i in xrange(self._M): + self.snks.append(gr.vector_sink_c()) + self.connect((self.pfb, i), self.snks[i]) + + +def main(): + tstart = time.time() + + tb = pfb_top_block() + tb.run() + + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig_in = pylab.figure(1, figsize=(16,9), facecolor="w") + fig1 = pylab.figure(2, figsize=(16,9), facecolor="w") + fig2 = pylab.figure(3, figsize=(16,9), facecolor="w") + + Ns = 1000 + Ne = 10000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._fs + + # Plot the input signal on its own figure + d = tb.snk_i.data()[Ns:Ne] + spin_f = fig_in.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + pin_f = spin_f.plot(f_in, X_in, "b") + spin_f.set_xlim([min(f_in), max(f_in)+1]) + spin_f.set_ylim([-200.0, 50.0]) + + spin_f.set_title("Input Signal", weight="bold") + spin_f.set_xlabel("Frequency (Hz)") + spin_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + spin_t = fig_in.add_subplot(2, 1, 2) + pin_t = spin_t.plot(t_in, x_in.real, "b") + pin_t = spin_t.plot(t_in, x_in.imag, "r") + + spin_t.set_xlabel("Time (s)") + spin_t.set_ylabel("Amplitude") + + Ncols = int(scipy.floor(scipy.sqrt(tb._M))) + Nrows = int(scipy.floor(tb._M / Ncols)) + if(tb._M % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = tb._fs / tb._M + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + for i in xrange(len(tb.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = tb.snks[i].data()[Ns:Ne] + + sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i) + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size)) + p2_f = sp1_f.plot(f_o, X_o, "b") + sp1_f.set_xlim([min(f_o), max(f_o)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title(("Channel %d" % i), weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i) + p2_o = sp2_o.plot(t_o, x_o.real, "b") + p2_o = sp2_o.plot(t_o, x_o.imag, "r") + sp2_o.set_xlim([min(t_o), max(t_o)+1]) + sp2_o.set_ylim([-2, 2]) + + sp2_o.set_title(("Channel %d" % i), weight="bold") + sp2_o.set_xlabel("Time (s)") + sp2_o.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-examples/python/pfb/chirp_channelize.py b/gnuradio-examples/python/pfb/chirp_channelize.py new file mode 100755 index 00000000..edebf5f5 --- /dev/null +++ b/gnuradio-examples/python/pfb/chirp_channelize.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import os, time +import scipy, pylab +from scipy import fftpack +from pylab import mlab + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 200000 # number of samples to use + self._fs = 9000 # initial sampling rate + self._M = 9 # Number of channels to channelize + + # Create a set of taps for the PFB channelizer + self._taps = gr.firdes.low_pass_2(1, self._fs, 500, 20, + attenuation_dB=10, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._M)) + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + repeated = True + if(repeated): + self.vco_input = gr.sig_source_f(self._fs, gr.GR_SIN_WAVE, 0.25, 110) + else: + amp = 100 + data = scipy.arange(0, amp, amp/float(self._N)) + self.vco_input = gr.vector_source_f(data, False) + + # Build a VCO controlled by either the sinusoid or single chirp tone + # Then convert this to a complex signal + self.vco = gr.vco_f(self._fs, 225, 1) + self.f2c = gr.float_to_complex() + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the channelizer filter + self.pfb = blks2.pfb_channelizer_ccf(self._M, self._taps) + + # Construct a vector sink for the input signal to the channelizer + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.vco_input, self.vco, self.f2c) + self.connect(self.f2c, self.head, self.pfb) + self.connect(self.f2c, self.snk_i) + + # Create a vector sink for each of M output channels of the filter and connect it + self.snks = list() + for i in xrange(self._M): + self.snks.append(gr.vector_sink_c()) + self.connect((self.pfb, i), self.snks[i]) + + +def main(): + tstart = time.time() + + tb = pfb_top_block() + tb.run() + + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig_in = pylab.figure(1, figsize=(16,9), facecolor="w") + fig1 = pylab.figure(2, figsize=(16,9), facecolor="w") + fig2 = pylab.figure(3, figsize=(16,9), facecolor="w") + fig3 = pylab.figure(4, figsize=(16,9), facecolor="w") + + Ns = 650 + Ne = 20000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._fs + + # Plot the input signal on its own figure + d = tb.snk_i.data()[Ns:Ne] + spin_f = fig_in.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + pin_f = spin_f.plot(f_in, X_in, "b") + spin_f.set_xlim([min(f_in), max(f_in)+1]) + spin_f.set_ylim([-200.0, 50.0]) + + spin_f.set_title("Input Signal", weight="bold") + spin_f.set_xlabel("Frequency (Hz)") + spin_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + spin_t = fig_in.add_subplot(2, 1, 2) + pin_t = spin_t.plot(t_in, x_in.real, "b") + pin_t = spin_t.plot(t_in, x_in.imag, "r") + + spin_t.set_xlabel("Time (s)") + spin_t.set_ylabel("Amplitude") + + Ncols = int(scipy.floor(scipy.sqrt(tb._M))) + Nrows = int(scipy.floor(tb._M / Ncols)) + if(tb._M % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = tb._fs / tb._M + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + for i in xrange(len(tb.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = tb.snks[i].data()[Ns:Ne] + + sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i) + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(X)) + f_o = freq + p2_f = sp1_f.plot(f_o, X_o, "b") + sp1_f.set_xlim([min(f_o), max(f_o)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title(("Channel %d" % i), weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i) + p2_o = sp2_o.plot(t_o, x_o.real, "b") + p2_o = sp2_o.plot(t_o, x_o.imag, "r") + sp2_o.set_xlim([min(t_o), max(t_o)+1]) + sp2_o.set_ylim([-2, 2]) + + sp2_o.set_title(("Channel %d" % i), weight="bold") + sp2_o.set_xlabel("Time (s)") + sp2_o.set_ylabel("Amplitude") + + + sp3 = fig3.add_subplot(1,1,1) + p3 = sp3.plot(t_o, x_o.real) + sp3.set_xlim([min(t_o), max(t_o)+1]) + sp3.set_ylim([-2, 2]) + + sp3.set_title("All Channels") + sp3.set_xlabel("Time (s)") + sp3.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-examples/python/pfb/decimate.py b/gnuradio-examples/python/pfb/decimate.py new file mode 100755 index 00000000..cb5d61b7 --- /dev/null +++ b/gnuradio-examples/python/pfb/decimate.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import os +import scipy, pylab +from scipy import fftpack +from pylab import mlab +import time + +#print os.getpid() +#raw_input() + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 10000000 # number of samples to use + self._fs = 10000 # initial sampling rate + self._decim = 20 # Decimation rate + + # Generate the prototype filter taps for the decimators with a 200 Hz bandwidth + self._taps = gr.firdes.low_pass_2(1, self._fs, 200, 150, + attenuation_dB=120, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._decim)) + print "Number of taps: ", len(self._taps) + print "Number of filters: ", self._decim + print "Taps per channel: ", tpc + + # Build the input signal source + # We create a list of freqs, and a sine wave is generated and added to the source + # for each one of these frequencies. + self.signals = list() + self.add = gr.add_cc() + freqs = [10, 20, 2040] + for i in xrange(len(freqs)): + self.signals.append(gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freqs[i], 1)) + self.connect(self.signals[i], (self.add,i)) + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct a PFB decimator filter + self.pfb = blks2.pfb_decimator_ccf(self._decim, self._taps, 0) + + # Construct a standard FIR decimating filter + self.dec = gr.fir_filter_ccf(self._decim, self._taps) + + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.add, self.head, self.pfb) + self.connect(self.add, self.snk_i) + + # Create the sink for the decimated siganl + self.snk = gr.vector_sink_c() + self.connect(self.pfb, self.snk) + + +def main(): + tb = pfb_top_block() + + tstart = time.time() + tb.run() + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig1 = pylab.figure(1, figsize=(16,9)) + fig2 = pylab.figure(2, figsize=(16,9)) + + Ns = 10000 + Ne = 10000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._fs + + # Plot the input to the decimator + + d = tb.snk_i.data()[Ns:Ns+Ne] + sp1_f = fig1.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + p1_f = sp1_f.plot(f_in, X_in, "b") + sp1_f.set_xlim([min(f_in), max(f_in)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title("Input Signal", weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + sp1_t = fig1.add_subplot(2, 1, 2) + p1_t = sp1_t.plot(t_in, x_in.real, "b") + p1_t = sp1_t.plot(t_in, x_in.imag, "r") + sp1_t.set_ylim([-tb._decim*1.1, tb._decim*1.1]) + + sp1_t.set_xlabel("Time (s)") + sp1_t.set_ylabel("Amplitude") + + + # Plot the output of the decimator + fs_o = tb._fs / tb._decim + + sp2_f = fig2.add_subplot(2, 1, 1) + d = tb.snk.data()[Ns:Ns+Ne] + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size)) + p2_f = sp2_f.plot(f_o, X_o, "b") + sp2_f.set_xlim([min(f_o), max(f_o)+1]) + sp2_f.set_ylim([-200.0, 50.0]) + + sp2_f.set_title("PFB Decimated Signal", weight="bold") + sp2_f.set_xlabel("Frequency (Hz)") + sp2_f.set_ylabel("Power (dBW)") + + + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_t = fig2.add_subplot(2, 1, 2) + p2_t = sp2_t.plot(t_o, x_o.real, "b-o") + p2_t = sp2_t.plot(t_o, x_o.imag, "r-o") + sp2_t.set_ylim([-2.5, 2.5]) + + sp2_t.set_xlabel("Time (s)") + sp2_t.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-examples/python/pfb/fmtest.py b/gnuradio-examples/python/pfb/fmtest.py new file mode 100755 index 00000000..97df0e0f --- /dev/null +++ b/gnuradio-examples/python/pfb/fmtest.py @@ -0,0 +1,197 @@ +#!/usr/bin/env python +# + + +from gnuradio import gr, eng_notation +from gnuradio import blks2 +from gnuradio.eng_option import eng_option +from optparse import OptionParser +import math, time, sys, scipy, pylab +from scipy import fftpack + +class fmtx(gr.hier_block2): + def __init__(self, lo_freq, audio_rate, if_rate): + + gr.hier_block2.__init__(self, "build_fm", + gr.io_signature(1, 1, gr.sizeof_float), # Input signature + gr.io_signature(1, 1, gr.sizeof_gr_complex)) # Output signature + + fmtx = blks2.nbfm_tx (audio_rate, if_rate, max_dev=5e3, tau=75e-6) + + # Local oscillator + lo = gr.sig_source_c (if_rate, # sample rate + gr.GR_SIN_WAVE, # waveform type + lo_freq, #frequency + 1.0, # amplitude + 0) # DC Offset + mixer = gr.multiply_cc () + + self.connect (self, fmtx, (mixer, 0)) + self.connect (lo, (mixer, 1)) + self.connect (mixer, self) + +class fmtest(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._nsamples = 1000000 + self._audio_rate = 8000 + + # Set up N channels with their own baseband and IF frequencies + self._N = 5 + chspacing = 16000 + freq = [10, 20, 30, 40, 50] + f_lo = [0, 1*chspacing, -1*chspacing, 2*chspacing, -2*chspacing] + + self._if_rate = 4*self._N*self._audio_rate + + # Create a signal source and frequency modulate it + self.sum = gr.add_cc () + for n in xrange(self._N): + sig = gr.sig_source_f(self._audio_rate, gr.GR_SIN_WAVE, freq[n], 0.5) + fm = fmtx(f_lo[n], self._audio_rate, self._if_rate) + self.connect(sig, fm) + self.connect(fm, (self.sum, n)) + + self.head = gr.head(gr.sizeof_gr_complex, self._nsamples) + self.snk_tx = gr.vector_sink_c() + self.channel = blks2.channel_model(0.1) + + self.connect(self.sum, self.head, self.channel, self.snk_tx) + + + # Design the channlizer + self._M = 10 + bw = chspacing/2.0 + t_bw = chspacing/10.0 + self._chan_rate = self._if_rate / self._M + self._taps = gr.firdes.low_pass_2(1, self._if_rate, bw, t_bw, + attenuation_dB=100, + window=gr.firdes.WIN_BLACKMAN_hARRIS) + tpc = math.ceil(float(len(self._taps)) / float(self._M)) + + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + self.pfb = blks2.pfb_channelizer_ccf(self._M, self._taps) + + self.connect(self.channel, self.pfb) + + # Create a file sink for each of M output channels of the filter and connect it + self.fmdet = list() + self.squelch = list() + self.snks = list() + for i in xrange(self._M): + self.fmdet.append(blks2.nbfm_rx(self._audio_rate, self._chan_rate)) + self.squelch.append(blks2.standard_squelch(self._audio_rate*10)) + self.snks.append(gr.vector_sink_f()) + self.connect((self.pfb, i), self.fmdet[i], self.squelch[i], self.snks[i]) + + def num_tx_channels(self): + return self._N + + def num_rx_channels(self): + return self._M + +def main(): + + fm = fmtest() + + tstart = time.time() + fm.run() + tend = time.time() + + if 1: + fig1 = pylab.figure(1, figsize=(12,10), facecolor="w") + fig2 = pylab.figure(2, figsize=(12,10), facecolor="w") + fig3 = pylab.figure(3, figsize=(12,10), facecolor="w") + + Ns = 10000 + Ne = 100000 + + fftlen = 8192 + winfunc = scipy.blackman + + # Plot transmitted signal + fs = fm._if_rate + + d = fm.snk_tx.data()[Ns:Ns+Ne] + sp1_f = fig1.add_subplot(2, 1, 1) + + X,freq = sp1_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + visible=False) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + p1_f = sp1_f.plot(f_in, X_in, "b") + sp1_f.set_xlim([min(f_in), max(f_in)+1]) + sp1_f.set_ylim([-120.0, 20.0]) + + sp1_f.set_title("Input Signal", weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + sp1_t = fig1.add_subplot(2, 1, 2) + p1_t = sp1_t.plot(t_in, x_in.real, "b-o") + #p1_t = sp1_t.plot(t_in, x_in.imag, "r-o") + sp1_t.set_ylim([-5, 5]) + + # Set up the number of rows and columns for plotting the subfigures + Ncols = int(scipy.floor(scipy.sqrt(fm.num_rx_channels()))) + Nrows = int(scipy.floor(fm.num_rx_channels() / Ncols)) + if(fm.num_rx_channels() % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = fm._audio_rate + for i in xrange(len(fm.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = fm.snks[i].data()[Ns:Ne] + + sp2_f = fig2.add_subplot(Nrows, Ncols, 1+i) + X,freq = sp2_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + visible=False) + #X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + X_o = 10.0*scipy.log10(abs(X)) + #f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size)) + f_o = scipy.arange(0, fs_o/2.0, fs_o/2.0/float(X_o.size)) + p2_f = sp2_f.plot(f_o, X_o, "b") + sp2_f.set_xlim([min(f_o), max(f_o)+0.1]) + sp2_f.set_ylim([-120.0, 20.0]) + sp2_f.grid(True) + + sp2_f.set_title(("Channel %d" % i), weight="bold") + sp2_f.set_xlabel("Frequency (kHz)") + sp2_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs_o + Tmax = len(d)*Ts + t_o = scipy.arange(0, Tmax, Ts) + + x_t = scipy.array(d) + sp2_t = fig3.add_subplot(Nrows, Ncols, 1+i) + p2_t = sp2_t.plot(t_o, x_t.real, "b") + p2_t = sp2_t.plot(t_o, x_t.imag, "r") + sp2_t.set_xlim([min(t_o), max(t_o)+1]) + sp2_t.set_ylim([-1, 1]) + + sp2_t.set_xlabel("Time (s)") + sp2_t.set_ylabel("Amplitude") + + + pylab.show() + + +if __name__ == "__main__": + main() diff --git a/gnuradio-examples/python/pfb/interpolate.py b/gnuradio-examples/python/pfb/interpolate.py new file mode 100755 index 00000000..a7a2522f --- /dev/null +++ b/gnuradio-examples/python/pfb/interpolate.py @@ -0,0 +1,226 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import os +import scipy, pylab +from scipy import fftpack +from pylab import mlab +import time + +#print os.getpid() +#raw_input() + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 100000 # number of samples to use + self._fs = 2000 # initial sampling rate + self._interp = 5 # Interpolation rate for PFB interpolator + self._ainterp = 5.5 # Resampling rate for the PFB arbitrary resampler + + # Frequencies of the signals we construct + freq1 = 100 + freq2 = 200 + + # Create a set of taps for the PFB interpolator + # This is based on the post-interpolation sample rate + self._taps = gr.firdes.low_pass_2(self._interp, self._interp*self._fs, freq2+50, 50, + attenuation_dB=120, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Create a set of taps for the PFB arbitrary resampler + # The filter size is the number of filters in the filterbank; 32 will give very low side-lobes, + # and larger numbers will reduce these even farther + # The taps in this filter are based on a sampling rate of the filter size since it acts + # internally as an interpolator. + flt_size = 32 + self._taps2 = gr.firdes.low_pass_2(flt_size, flt_size*self._fs, freq2+50, 150, + attenuation_dB=120, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._interp)) + print "Number of taps: ", len(self._taps) + print "Number of filters: ", self._interp + print "Taps per channel: ", tpc + + # Create a couple of signals at different frequencies + self.signal1 = gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freq1, 0.5) + self.signal2 = gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freq2, 0.5) + self.signal = gr.add_cc() + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the PFB interpolator filter + self.pfb = blks2.pfb_interpolator_ccf(self._interp, self._taps) + + # Construct the PFB arbitrary resampler filter + self.pfb_ar = blks2.pfb_arb_resampler_ccf(self._ainterp, self._taps2, flt_size) + self.snk_i = gr.vector_sink_c() + + #self.pfb_ar.pfb.print_taps() + #self.pfb.pfb.print_taps() + + # Connect the blocks + self.connect(self.signal1, self.head, (self.signal,0)) + self.connect(self.signal2, (self.signal,1)) + self.connect(self.signal, self.pfb) + self.connect(self.signal, self.pfb_ar) + self.connect(self.signal, self.snk_i) + + # Create the sink for the interpolated signals + self.snk1 = gr.vector_sink_c() + self.snk2 = gr.vector_sink_c() + self.connect(self.pfb, self.snk1) + self.connect(self.pfb_ar, self.snk2) + + +def main(): + tb = pfb_top_block() + + tstart = time.time() + tb.run() + tend = time.time() + print "Run time: %f" % (tend - tstart) + + + if 1: + fig1 = pylab.figure(1, figsize=(12,10), facecolor="w") + fig2 = pylab.figure(2, figsize=(12,10), facecolor="w") + fig3 = pylab.figure(3, figsize=(12,10), facecolor="w") + + Ns = 10000 + Ne = 10000 + + fftlen = 8192 + winfunc = scipy.blackman + + # Plot input signal + fs = tb._fs + + d = tb.snk_i.data()[Ns:Ns+Ne] + sp1_f = fig1.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + p1_f = sp1_f.plot(f_in, X_in, "b") + sp1_f.set_xlim([min(f_in), max(f_in)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + + sp1_f.set_title("Input Signal", weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + sp1_t = fig1.add_subplot(2, 1, 2) + p1_t = sp1_t.plot(t_in, x_in.real, "b-o") + #p1_t = sp1_t.plot(t_in, x_in.imag, "r-o") + sp1_t.set_ylim([-2.5, 2.5]) + + sp1_t.set_title("Input Signal", weight="bold") + sp1_t.set_xlabel("Time (s)") + sp1_t.set_ylabel("Amplitude") + + + # Plot output of PFB interpolator + fs_int = tb._fs*tb._interp + + sp2_f = fig2.add_subplot(2, 1, 1) + d = tb.snk1.data()[Ns:Ns+(tb._interp*Ne)] + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_o = scipy.arange(-fs_int/2.0, fs_int/2.0, fs_int/float(X_o.size)) + p2_f = sp2_f.plot(f_o, X_o, "b") + sp2_f.set_xlim([min(f_o), max(f_o)+1]) + sp2_f.set_ylim([-200.0, 50.0]) + + sp2_f.set_title("Output Signal from PFB Interpolator", weight="bold") + sp2_f.set_xlabel("Frequency (Hz)") + sp2_f.set_ylabel("Power (dBW)") + + Ts_int = 1.0/fs_int + Tmax = len(d)*Ts_int + + t_o = scipy.arange(0, Tmax, Ts_int) + x_o1 = scipy.array(d) + sp2_t = fig2.add_subplot(2, 1, 2) + p2_t = sp2_t.plot(t_o, x_o1.real, "b-o") + #p2_t = sp2_t.plot(t_o, x_o.imag, "r-o") + sp2_t.set_ylim([-2.5, 2.5]) + + sp2_t.set_title("Output Signal from PFB Interpolator", weight="bold") + sp2_t.set_xlabel("Time (s)") + sp2_t.set_ylabel("Amplitude") + + + # Plot output of PFB arbitrary resampler + fs_aint = tb._fs * tb._ainterp + + sp3_f = fig3.add_subplot(2, 1, 1) + d = tb.snk2.data()[Ns:Ns+(tb._interp*Ne)] + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_o = scipy.arange(-fs_aint/2.0, fs_aint/2.0, fs_aint/float(X_o.size)) + p3_f = sp3_f.plot(f_o, X_o, "b") + sp3_f.set_xlim([min(f_o), max(f_o)+1]) + sp3_f.set_ylim([-200.0, 50.0]) + + sp3_f.set_title("Output Signal from PFB Arbitrary Resampler", weight="bold") + sp3_f.set_xlabel("Frequency (Hz)") + sp3_f.set_ylabel("Power (dBW)") + + Ts_aint = 1.0/fs_aint + Tmax = len(d)*Ts_aint + + t_o = scipy.arange(0, Tmax, Ts_aint) + x_o2 = scipy.array(d) + sp3_f = fig3.add_subplot(2, 1, 2) + p3_f = sp3_f.plot(t_o, x_o2.real, "b-o") + p3_f = sp3_f.plot(t_o, x_o1.real, "m-o") + #p3_f = sp3_f.plot(t_o, x_o2.imag, "r-o") + sp3_f.set_ylim([-2.5, 2.5]) + + sp3_f.set_title("Output Signal from PFB Arbitrary Resampler", weight="bold") + sp3_f.set_xlabel("Time (s)") + sp3_f.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gr-utils/src/python/gr_plot_const.py b/gr-utils/src/python/gr_plot_const.py index 35530995..ec2272c7 100755 --- a/gr-utils/src/python/gr_plot_const.py +++ b/gr-utils/src/python/gr_plot_const.py @@ -35,9 +35,6 @@ except ImportError: from optparse import OptionParser -matplotlib.interactive(True) -matplotlib.use('TkAgg') - class draw_constellation: def __init__(self, filename, options): self.hfile = open(filename, "r") @@ -139,8 +136,9 @@ class draw_constellation: draw() def zoom(self, event): - newxlim = self.sp_iq.get_xlim() - if(newxlim.all() != self.xlim.all()): + newxlim = scipy.array(self.sp_iq.get_xlim()) + curxlim = scipy.array(self.xlim) + if(newxlim.all() != curxlim.all()): self.xlim = newxlim r = self.reals[int(ceil(self.xlim[0])) : int(ceil(self.xlim[1]))] i = self.imags[int(ceil(self.xlim[0])) : int(ceil(self.xlim[1]))] diff --git a/gr-utils/src/python/gr_plot_fft.py b/gr-utils/src/python/gr_plot_fft.py index 59a3f286..a9c1417f 100755 --- a/gr-utils/src/python/gr_plot_fft.py +++ b/gr-utils/src/python/gr_plot_fft.py @@ -20,14 +20,6 @@ # Boston, MA 02110-1301, USA. # -try: - import matplotlib - matplotlib.use('TkAgg') - matplotlib.interactive(True) -except ImportError: - print "Please install Matplotlib to run this script (http://matplotlib.sourceforge.net/)" - raise SystemExit, 1 - try: import scipy from scipy import fftpack @@ -164,8 +156,9 @@ class gr_plot_fft: draw() def zoom(self, event): - newxlim = self.sp_iq.get_xlim() - if(newxlim.all() != self.xlim.all()): + newxlim = scipy.array(self.sp_iq.get_xlim()) + curxlim = scipy.array(self.xlim) + if(newxlim.all() != curxlim.all()): self.xlim = newxlim xmin = max(0, int(ceil(self.sample_rate*(self.xlim[0] - self.position)))) xmax = min(int(ceil(self.sample_rate*(self.xlim[1] - self.position))), len(self.iq)) diff --git a/gr-utils/src/python/gr_plot_iq.py b/gr-utils/src/python/gr_plot_iq.py index 2a4142a8..371ce3b7 100755 --- a/gr-utils/src/python/gr_plot_iq.py +++ b/gr-utils/src/python/gr_plot_iq.py @@ -34,10 +34,7 @@ except ImportError: from optparse import OptionParser -matplotlib.interactive(True) -matplotlib.use('TkAgg') - -class draw_fft: +class draw_iq: def __init__(self, filename, options): self.hfile = open(filename, "r") self.block_length = options.block @@ -168,7 +165,7 @@ def main(): raise SystemExit, 1 filename = args[0] - dc = draw_fft(filename, options) + dc = draw_iq(filename, options) if __name__ == "__main__": try: diff --git a/gr-utils/src/python/gr_plot_psd.py b/gr-utils/src/python/gr_plot_psd.py index 669d7b57..0e3dbecd 100755 --- a/gr-utils/src/python/gr_plot_psd.py +++ b/gr-utils/src/python/gr_plot_psd.py @@ -20,14 +20,6 @@ # Boston, MA 02110-1301, USA. # -try: - import matplotlib - matplotlib.use('TkAgg') - matplotlib.interactive(True) -except ImportError: - print "Please install Matplotlib to run this script (http://matplotlib.sourceforge.net/)" - raise SystemExit, 1 - try: import scipy from scipy import fftpack @@ -187,8 +179,9 @@ class gr_plot_psd: draw() def zoom(self, event): - newxlim = self.sp_iq.get_xlim() - if(newxlim.all() != self.xlim.all()): + newxlim = scipy.array(self.sp_iq.get_xlim()) + curxlim = scipy.array(self.xlim) + if(newxlim.all() != curxlim.all()): self.xlim = newxlim xmin = max(0, int(ceil(self.sample_rate*(self.xlim[0] - self.position)))) xmax = min(int(ceil(self.sample_rate*(self.xlim[1] - self.position))), len(self.iq)) diff --git a/gr-utils/src/python/plot_data.py b/gr-utils/src/python/plot_data.py index 7c79e171..08cdd603 100644 --- a/gr-utils/src/python/plot_data.py +++ b/gr-utils/src/python/plot_data.py @@ -33,9 +33,6 @@ except ImportError: from optparse import OptionParser -matplotlib.interactive(True) -matplotlib.use('TkAgg') - class plot_data: def __init__(self, datatype, filenames, options): self.hfile = list() -- 2.30.2