[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Commit-gnuradio] r11584 - in gnuradio/branches/developers/n4hy/pfb_iir2
|
From: |
n4hy |
|
Subject: |
[Commit-gnuradio] r11584 - in gnuradio/branches/developers/n4hy/pfb_iir2: gnuradio-core/src/lib/filter gnuradio-examples/python/pfb |
|
Date: |
Wed, 12 Aug 2009 07:30:14 -0600 (MDT) |
Author: n4hy
Date: 2009-08-12 07:30:13 -0600 (Wed, 12 Aug 2009)
New Revision: 11584
Added:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/Makefile.am
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/channelize.py
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/chirp_channelize.py
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/decimate.py
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/interpolate.py
Log:
putting pfb to look at half band and different arbitrary resampling structure
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc
2009-08-12 13:30:13 UTC (rev 11584)
@@ -30,23 +30,38 @@
#include <gr_io_signature.h>
gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate,
- const
std::vector<float> &taps)
+ const
std::vector<float> &taps,
+ unsigned int
filter_size)
{
- return gr_pfb_arb_resampler_ccf_sptr (new gr_pfb_arb_resampler_ccf (rate,
taps));
+ return gr_pfb_arb_resampler_ccf_sptr (new gr_pfb_arb_resampler_ccf (rate,
taps,
+
filter_size));
}
gr_pfb_arb_resampler_ccf::gr_pfb_arb_resampler_ccf (float rate,
- const std::vector<float>
&taps)
+ const std::vector<float>
&taps,
+ unsigned int filter_size)
: gr_block ("pfb_arb_resampler_ccf",
gr_make_io_signature (1, 1, sizeof(gr_complex)),
gr_make_io_signature (1, 1, sizeof(gr_complex))),
d_updated (false)
{
- d_int_rate = 32;
+ /* The number of filters is specified by the user as the filter size;
+ this is also the interpolation rate of the filter. We use it and the
+ rate provided to determine the decimation rate. This acts as a
+ rational resampler. The flt_rate is calculated as the residual
+ between the integer decimation rate and the real decimation rate and
+ will be used to determine to interpolation point of the resampling
+ process.
+ */
+ d_int_rate = filter_size;
d_dec_rate = (unsigned int)floor(d_int_rate/rate);
d_flt_rate = (d_int_rate/rate) - d_dec_rate;
+
+ // The accumulator keeps track of overflow to increment the stride correctly.
d_acc = 0;
+
+ // Store the last filter between calls to work
d_last_filter = 0;
d_filters = std::vector<gr_fir_ccf*>(d_int_rate);
@@ -78,13 +93,21 @@
// Create d_numchan vectors to store each channel's taps
d_taps.resize(d_int_rate);
+
+ // Make a vector of the taps plus fill it out with 0's to fill
+ // each polyphase filter with exactly d_taps_per_filter
+ std::vector<float> tmp_taps;
+ tmp_taps = taps;
+ while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
+ tmp_taps.push_back(0.0);
+ }
// Partition the filter
for(i = 0; i < d_int_rate; i++) {
// Each channel uses all d_taps_per_filter with 0's if not enough taps to
fill out
d_taps[i] = std::vector<float>(d_taps_per_filter, 0);
for(j = 0; j < d_taps_per_filter; j++) {
- d_taps[i][j] = taps[i + j*d_int_rate]; // add taps to channels in
reverse order
+ d_taps[i][j] = tmp_taps[i + j*d_int_rate]; // add taps to channels in
reverse order
}
// Build a filter for each channel and add it's taps to it
@@ -124,34 +147,48 @@
return 0; // history requirements may have changed.
}
- int i = 0, j = 0, count = 0;
+ int i = 0, j, count = 0;
gr_complex o0, o1;
+ // Restore the last filter position
j = d_last_filter;
- while(i < noutput_items) {
+
+ // produce output as long as we can and there are enough input samples
+ while((i < noutput_items) && (count < ninput_items[0]-1)) {
+
// start j by wrapping around mod the number of channels
-
- while(j < d_int_rate) {
+ while((j < d_int_rate) && (i < noutput_items)) {
// Take the current filter output
o0 = d_filters[j]->filter(&in[count]);
// Take the next filter output; wrap around to 0 if necessary
if(j+1 == d_int_rate)
+ // Use the sample of the next input item through the first filter
o1 = d_filters[0]->filter(&in[count+1]);
- else
- o1 = d_filters[j]->filter(&in[count]);
+ else {
+ // Use the sample from the current input item through the nex filter
+ o1 = d_filters[j+1]->filter(&in[count]);
+ }
- out[i] = o0 + (o1 - o0)*d_flt_rate;
+ //out[i] = o0; // nearest-neighbor approach
+ out[i] = o0 + (o1 - o0)*d_acc; // linearly interpolate between
samples
i++;
-
- d_acc += d_dec_rate + d_flt_rate;
- j = (int)floor(d_acc);
+
+ // Accumulate the position in the stream for the interpolated point.
+ // If it goes above 1, roll around to zero and increment the stride
+ // length this time by the decimation rate plus 1 for the increment
+ // due to the acculated position.
+ d_acc += d_flt_rate;
+ j += d_dec_rate + (int)floor(d_acc);
+ d_acc = fmodf(d_acc, 1.0);
}
- count++;
- j = j % d_int_rate;
- d_acc = d_acc - floor(d_acc);
+ if(i < noutput_items) { // keep state for next entry
+ count++; // we have fully consumed another
input
+ j = j % d_int_rate; // roll filter around
+ }
}
+ // Store the current filter position
d_last_filter = j;
consume_each(count);
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h
2009-08-12 13:30:13 UTC (rev 11584)
@@ -29,41 +29,124 @@
class gr_pfb_arb_resampler_ccf;
typedef boost::shared_ptr<gr_pfb_arb_resampler_ccf>
gr_pfb_arb_resampler_ccf_sptr;
gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate,
- const
std::vector<float> &taps);
+ const
std::vector<float> &taps,
+ unsigned int
filter_size=32);
class gr_fir_ccf;
/*!
- * \brief FIR filter with gr_complex input, gr_complex output and float taps
+ * \class gr_pfb_arb_resampler_ccf
+ *
+ * \brief Polyphase filterbank arbitrary resampler with
+ * gr_complex input, gr_complex output and float taps
+ *
* \ingroup filter_blk
+ *
+ * This block takes in a signal stream and performs arbitrary
+ * resampling. The resampling rate can be any real
+ * number <EM>r</EM>. The resampling is done by constructing
+ * <EM>N</EM> filters where <EM>N</EM> is the interpolation rate. We
+ * then calculate <EM>D</EM> where <EM>D = floor(N/r)</EM>.
+ *
+ * Using <EM>N</EM> and <EM>D</EM>, we can perform rational resampling
+ * where <EM>N/D</EM> is a rational number close to the input rate
+ * <EM>r</EM> where we have <EM>N</EM> filters and we cycle through
+ * them as a polyphase filterbank with a stride of <EM>D</EM> so that
+ * <EM>i+1 = (i + D) % N</EM>.
+ *
+ * To get the arbitrary rate, we want to interpolate between two
+ * points. For each value out, we take an output from the current
+ * filter, <EM>i</EM>, and the next filter <EM>i+1</EM> and then
+ * linearly interpolate between the two based on the real resampling
+ * rate we want.
+ *
+ * The linear interpolation only provides us with an approximation to
+ * the real sampling rate specified. The error is a quantization error
+ * between the two filters we used as our interpolation points. To
+ * this end, the number of filters, <EM>N</EM>, used determines the
+ * quantization error; the larger <EM>N</EM>, the smaller the
+ * noise. You can design for a specified noise floor by setting the
+ * filter size (parameters <EM>filter_size</EM>). The size defaults to
+ * 32 filters, which is about as good as most implementations need.
+ *
+ * The trick with designing this filter is in how to specify the taps
+ * of the prototype filter. Like the PFB interpolator, the taps are
+ * specified using the interpolated filter rate. In this case, that
+ * rate is the input sample rate multiplied by the number of filters
+ * in the filterbank, which is also the interpolation rate. All other
+ * values should be relative to this rate.
+ *
+ * For example, for a 32-filter arbitrary resampler and using the
+ * GNU Radio's firdes utility to build the filter, we build a low-pass
+ * filter with a sampling rate of <EM>fs</EM>, a 3-dB bandwidth of
+ * <EM>BW</EM> and a transition bandwidth of <EM>TB</EM>. We can also
+ * specify the out-of-band attenuation to use, <EM>ATT</EM>, and the
+ * filter window function (a Blackman-harris window in this case). The
+ * first input is the gain of the filter, which we specify here as the
+ * interpolation rate (<EM>32</EM>).
+ *
+ * <B><EM>self._taps = gr.firdes.low_pass_2(32, 32*fs, BW, TB,
+ * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS)</EM></B>
+ *
+ * The theory behind this block can be found in Chapter 7.5 of
+ * the following book.
+ *
+ * <B><EM>f. harris, Multirate Signal Processing for Communication
+ * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004.</EM></B>
*/
+
class gr_pfb_arb_resampler_ccf : public gr_block
{
private:
+ /*!
+ * Build the polyphase filterbank arbitray resampler.
+ * \param rate (float) Specifies the resampling rate to use
+ * \param taps (vector/list of floats) The prototype filter to populate the
filterbank. The taps
+ * should be generated at the
filter_size sampling rate.
+ * \param filter_size (unsigned int) The number of filters in the filter
bank. This is directly
+ related to quantization noise
introduced during the resampling.
+ Defaults to 32 filters.
+ */
friend gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float
rate,
- const
std::vector<float> &taps);
+ const
std::vector<float> &taps,
+ unsigned
int filter_size);
std::vector<gr_fir_ccf*> d_filters;
std::vector< std::vector<float> > d_taps;
- unsigned int d_int_rate;
- unsigned int d_dec_rate;
- float d_flt_rate;
+ unsigned int d_int_rate; // the number of filters
(interpolation rate)
+ unsigned int d_dec_rate; // the stride through the
filters (decimation rate)
+ float d_flt_rate; // residual rate for the
linear interpolation
float d_acc;
unsigned int d_last_filter;
unsigned int d_taps_per_filter;
bool d_updated;
/*!
- * Construct a Polyphase filterbank for channelization with the given
- * number of channels and taps
+ * Build the polyphase filterbank arbitray resampler.
+ * \param rate (float) Specifies the resampling rate to use
+ * \param taps (vector/list of floats) The prototype filter to populate the
filterbank. The taps
+ * should be generated at the
filter_size sampling rate.
+ * \param filter_size (unsigned int) The number of filters in the filter
bank. This is directly
+ related to quantization noise
introduced during the resampling.
+ Defaults to 32 filters.
*/
gr_pfb_arb_resampler_ccf (float rate,
- const std::vector<float> &taps);
+ const std::vector<float> &taps,
+ unsigned int filter_size);
public:
~gr_pfb_arb_resampler_ccf ();
+ /*!
+ * Resets the filterbank's filter taps with the new prototype filter
+ * \param taps (vector/list of floats) The prototype filter to populate
the filterbank. The taps
+ * should be generated at the
interpolated sampling rate.
+ */
void set_taps (const std::vector<float> &taps);
+
+ /*!
+ * Print all of the filterbank taps to screen.
+ */
void print_taps();
int general_work (int noutput_items,
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i
2009-08-12 13:30:13 UTC (rev 11584)
@@ -23,13 +23,15 @@
GR_SWIG_BLOCK_MAGIC(gr,pfb_arb_resampler_ccf);
gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate,
- const
std::vector<float> &taps);
+ const
std::vector<float> &taps,
+ unsigned int
filter_size=32);
class gr_pfb_arb_resampler_ccf : public gr_block
{
private:
gr_pfb_arb_resampler_ccf (float rate,
- const std::vector<float> &taps);
+ const std::vector<float> &taps,
+ unsigned int filter_size);
public:
~gr_pfb_arb_resampler_ccf ();
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc
2009-08-12 13:30:13 UTC (rev 11584)
@@ -76,18 +76,26 @@
d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_numchans);
// Create d_numchan vectors to store each channel's taps
- std::vector< std::vector<float> > vtaps(d_numchans);
+ 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<float> 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
- vtaps[i] = std::vector<float>(d_taps_per_filter, 0);
+ d_taps[i] = std::vector<float>(d_taps_per_filter, 0);
for(j = 0; j < d_taps_per_filter; j++) {
- vtaps[i][j] = taps[i + j*d_numchans]; // add taps to channels in
reverse order
+ 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(vtaps[i]);
+ d_filters[i]->set_taps(d_taps[i]);
}
// Set the history to ensure enough input items for each filter
@@ -96,6 +104,20 @@
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,
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h
2009-08-12 13:30:13 UTC (rev 11584)
@@ -34,25 +34,90 @@
class gr_fir_ccf;
class gri_fft_complex;
+
/*!
- * \brief FIR filter with gr_complex input, gr_complex output and float taps
+ * \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 <EM>M</EM>
+ * channels of equal bandwidth. Each of the resulting channels is
+ * decimated to the new rate that is the input sampling rate
+ * <EM>fs</EM> divided by the number of channels, <EM>M</EM>.
+ *
+ * The PFB channelizer code takes the taps generated above and builds
+ * a set of filters. The set contains <EM>M</EM> 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 <EM>i</EM> and performs the inner
+ * product calculation to <EM>i+(n-1)</EM> where <EM>n</EM> 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 <EM>M</EM> 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 <EM>i</EM> in
+ * the vector is the next sample from the <EM>i</EM>th channel. This
+ * is most easily handled by sending the output to a
+ * gr_vector_to_streams block to handle the conversion and passing
+ * <EM>M</EM> 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
+ * <EM>M</EM> 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
+ * <EM>fs</EM>, a 3-dB bandwidth of <EM>BW</EM> and a transition
+ * bandwidth of <EM>TB</EM>. We can also specify the out-of-band
+ * attenuation to use, <EM>ATT</EM>, and the filter window
+ * function (a Blackman-harris window in this case). The first input
+ * is the gain of the filter, which we specify here as unity.
+ *
+ * <B><EM>self._taps = gr.firdes.low_pass_2(1, fs, BW, TB,
+ * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS)</EM></B>
+ *
+ * The theory behind this block can be found in Chapter 6 of
+ * the following book.
+ *
+ * <B><EM>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
<EM>M</EM>
+ * \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<float> &taps);
std::vector<gr_fir_ccf*> d_filters;
+ std::vector< std::vector<float> > d_taps;
gri_fft_complex *d_fft;
unsigned int d_numchans;
unsigned int d_taps_per_filter;
bool d_updated;
/*!
- * Construct a Polyphase filterbank for channelization with the given
- * number of channels and taps
+ * Build the polyphase filterbank decimator.
+ * \param numchans (unsigned integer) Specifies the number of channels
<EM>M</EM>
+ * \param taps (vector/list of floats) The prototype filter to populate
the filterbank.
*/
gr_pfb_channelizer_ccf (unsigned int numchans,
const std::vector<float> &taps);
@@ -60,7 +125,16 @@
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<float> &taps);
+
+ /*!
+ * Print all of the filterbank taps to screen.
+ */
+ void print_taps();
int work (int noutput_items,
gr_vector_const_void_star &input_items,
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.cc
2009-08-12 13:30:13 UTC (rev 11584)
@@ -82,18 +82,26 @@
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<float> > 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<float> 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
- vtaps[i] = std::vector<float>(d_taps_per_filter, 0);
+ d_taps[i] = std::vector<float>(d_taps_per_filter, 0);
for(j = 0; j < d_taps_per_filter; j++) {
- vtaps[i][j] = taps[i + j*d_rate]; // add taps to channels in reverse
order
+ 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(vtaps[i]);
+ d_filters[i]->set_taps(d_taps[i]);
}
// Set the history to ensure enough input items for each filter
@@ -102,7 +110,21 @@
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,
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_decimator_ccf.h
2009-08-12 13:30:13 UTC (rev 11584)
@@ -36,17 +36,77 @@
class gri_fft_complex;
/*!
- * \brief FIR filter with gr_complex input, gr_complex output and float taps
+ * \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 <EM>channel</EM> parameter specifies which channel to use since
+ * this class is capable of bandpass decimation. Given a complex input
+ * stream at a sampling rate of <EM>fs</EM> and a decimation rate of
+ * <EM>decim</EM>, the input frequency domain is split into
+ * <EM>decim</EM> 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
+ * <EM>fs</EM>, a 3-dB bandwidth of <EM>BW</EM> and a transition
+ * bandwidth of <EM>TB</EM>. We can also specify the out-of-band
+ * attenuation to use, <EM>ATT</EM>, and the filter window
+ * function (a Blackman-harris window in this case). The first input
+ * is the gain of the filter, which we specify here as unity.
+ *
+ * <B><EM>self._taps = gr.firdes.low_pass_2(1, fs, BW, TB,
+ * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS)</EM></B>
+ *
+ * The PFB decimator code takes the taps generated above and builds a
+ * set of filters. The set contains <EM>decim</EM> 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.
+ *
+ * <B><EM>f. harris, Multirate Signal Processing for Communication
+ * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004.</EM></B>
*/
+
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<float> &taps,
unsigned int
channel);
std::vector<gr_fir_ccf*> d_filters;
+ std::vector< std::vector<float> > d_taps;
gri_fft_complex *d_fft;
unsigned int d_rate;
unsigned int d_chan;
@@ -55,8 +115,10 @@
gr_complex *d_rotator;
/*!
- * Construct a Polyphase filterbank for channelization with the given
- * number of channels and taps
+ * 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<float> &taps,
@@ -65,9 +127,19 @@
public:
~gr_pfb_decimator_ccf ();
- void set_taps (const std::vector<float> &taps);
- //void set_channel (unsigned int channel);
+ /*!
+ * Resets the filterbank's filter taps with the new prototype filter
+ * \param taps (vector/list of floats) The prototype filter to populate
the filterbank.
+ */
+ void set_taps (const std::vector<float> &taps);
+ /*!
+ * 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);
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.cc
2009-08-12 13:30:13 UTC (rev 11584)
@@ -75,13 +75,21 @@
// Create d_numchan vectors to store each channel's taps
//std::vector< std::vector<float> > 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<float> 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<float>(d_taps_per_filter, 0);
for(j = 0; j < d_taps_per_filter; j++) {
- d_taps[i][j] = taps[i + j*d_rate]; // add taps to channels in reverse
order
+ 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
@@ -103,7 +111,7 @@
for(j = 0; j < d_taps_per_filter; j++) {
printf(" %.4e", d_taps[i][j]);
}
- printf("]\n");
+ printf("]\n\n");
}
}
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-core/src/lib/filter/gr_pfb_interpolator_ccf.h
2009-08-12 13:30:13 UTC (rev 11584)
@@ -34,12 +34,60 @@
class gr_fir_ccf;
/*!
- * \brief FIR filter with gr_complex input, gr_complex output and float taps
+ * \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
+ * <EM>fs</EM>, a 3-dB bandwidth of <EM>BW</EM> and a transition
+ * bandwidth of <EM>TB</EM>. We can also specify the out-of-band
+ * attenuation to use, 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).
+ *
+ * <B><EM>self._taps = gr.firdes.low_pass_2(interp, interp*fs, BW, TB,
+ * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS)</EM></B>
+ *
+ * The PFB interpolator code takes the taps generated above and builds
+ * a set of filters. The set contains <EM>interp</EM> 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.
+ *
+ * <B><EM>f. harris, <EM>Multirate Signal Processing for Communication
+ * Systems</EM>," Upper Saddle River, NJ: Prentice Hall,
+ * Inc. 2004.</EM></B>
*/
+
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<float> &taps);
@@ -50,8 +98,10 @@
bool d_updated;
/*!
- * Construct a Polyphase filterbank for channelization with the given
- * number of channels and taps
+ * 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<float> &taps);
@@ -59,7 +109,16 @@
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<float> &taps);
+
+ /*!
+ * Print all of the filterbank taps to screen.
+ */
void print_taps();
int work (int noutput_items,
Added:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/Makefile.am
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/Makefile.am
(rev 0)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/Makefile.am
2009-08-12 13:30:13 UTC (rev 11584)
@@ -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
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/channelize.py
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/channelize.py
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/channelize.py
2009-08-12 13:30:13 UTC (rev 11584)
@@ -1,49 +1,72 @@
#!/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 math
-import os
+import os, time
import scipy, pylab
from scipy import fftpack
-import time
+from pylab import mlab
-#print os.getpid()
-#raw_input()
-
class pfb_top_block(gr.top_block):
def __init__(self):
gr.top_block.__init__(self)
- self._N = 2000000
- self._fs = 9000
- self._Tmax = self._N * (1.0/self._fs)
- self._M = 9
- self._taps = gr.firdes.low_pass_2(1, self._fs, 475,50, 100,5)
+ self._N = 2000000 # number of samples to use
+ self._fs = 9000 # initial sampling rate
+ self._M = 9 # Number of channels to channelize
- fc = 200
+ # 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)
- tpc = math.ceil(float(len(self._taps)) / float(self._M))
-
+ # 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, 4300]
+ 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 file sink for each of M output channels of the filter and
connect it
+ # 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())
@@ -75,9 +98,9 @@
d = tb.snk_i.data()[Ns:Ne]
spin_f = fig_in.add_subplot(2, 1, 1)
- X,freq = spin_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs,
- window = lambda d: d*winfunc(fftlen),
- visible=False)
+ 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")
@@ -88,7 +111,11 @@
spin_f.set_xlabel("Frequency (Hz)")
spin_f.set_ylabel("Power (dBW)")
- t_in = scipy.arange(0, tb._Tmax, tb._Tmax/float(len(d)))
+
+ 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")
@@ -105,6 +132,8 @@
# 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
@@ -112,9 +141,9 @@
d = tb.snks[i].data()[Ns:Ne]
sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i)
- X,freq = sp1_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o,
- window = lambda d: d*winfunc(fftlen),
- visible=False)
+ 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")
@@ -125,15 +154,15 @@
sp1_f.set_xlabel("Frequency (Hz)")
sp1_f.set_ylabel("Power (dBW)")
-
x_o = scipy.array(d)
- t_o = scipy.arange(0, tb._Tmax, tb._Tmax/float(x_o.size))
+ 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")
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/chirp_channelize.py
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/chirp_channelize.py
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/chirp_channelize.py
2009-08-12 13:30:13 UTC (rev 11584)
@@ -1,49 +1,76 @@
#!/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 math
-import os
+import os, time
import scipy, pylab
from scipy import fftpack
-import time
+from pylab import mlab
-#print os.getpid()
-#raw_input()
-
class pfb_top_block(gr.top_block):
def __init__(self):
gr.top_block.__init__(self)
- self._N = 200000
- self._fs = 9000
- self._Tmax = self._N * (1.0/self._fs)
- self._M = 9
- self._taps = gr.firdes.low_pass_2(1, self._fs, 800,100 ,90,5)
- fc = 200
+ self._N = 200000 # number of samples to use
+ self._fs = 9000 # initial sampling rate
+ self._M = 9 # Number of channels to channelize
- tpc = math.ceil(float(len(self._taps)) / float(self._M))
+ # 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
- data = scipy.arange(0, 1000, 100.0/float(self._N))
-
- #self.vco_input = gr.sig_source_f(self._fs, gr.GR_SIN_WAVE, 0.25, 100)
- self.vco_input = gr.vector_source_f(data, False)
+ 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 rest
+ # 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 file sink for each of M output channels of the filter and
connect it
+ # 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())
@@ -60,42 +87,100 @@
print "Run time: %f" % (tend - tstart)
if 1:
- fig1 = pylab.figure(1, figsize=(16,9))
- fig2 = pylab.figure(2, figsize=(16,9))
+ 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 = 100
- Ne = tb._N
+ 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]
- f_in = scipy.arange(-tb._fs/2.0, tb._fs/2.0, tb._fs/float(len(d)))
- X_in = 10.0*scipy.log10(fftpack.fftshift(fftpack.fft(d, f_in.size)))
- sp1_in = fig1.add_subplot(4, 3, 1)
- p1_in = sp1_in.plot(f_in, X_in)
+ 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, tb._Tmax, tb._Tmax/float(len(d)))
+ t_in = scipy.arange(0, Tmax, Ts)
x_in = scipy.array(d)
- sp2_in = fig2.add_subplot(4, 3, 1)
- p2_in = sp2_in.plot(t_in, x_in.real, "b")
- p2_in = sp2_in.plot(t_in, x_in.imag, "r")
-
+ 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]
- f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(len(d)))
- x_o = 10.0*scipy.log10(fftpack.fftshift(fftpack.fft(d)))
- sp1_o = fig1.add_subplot(4, 3, 2+i)
- p1_o = sp1_o.plot(f_o, x_o)
- sp1_o.set_ylim([-50.0, 30.0])
-
+
+ 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, tb._Tmax, tb._Tmax/float(x_o.size))
- sp2_o = fig2.add_subplot(4, 3, 2+i)
+ 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_ylim([-1.0, 1.0])
-
+ 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()
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/decimate.py
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/decimate.py
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/decimate.py
2009-08-12 13:30:13 UTC (rev 11584)
@@ -1,10 +1,30 @@
#!/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 math
import os
import scipy, pylab
from scipy import fftpack
+from pylab import mlab
import time
#print os.getpid()
@@ -14,28 +34,38 @@
def __init__(self):
gr.top_block.__init__(self)
- self._N = 10000000
- self._fs = 10000
- self._Tmax = self._N * (1.0/self._fs)
- self._decim = 20
- self._taps = gr.firdes.low_pass_2(1, self._fs, 200, 50,120,5)
+ 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)
- tpc = math.ceil(float(len(self._taps)) / float(self._decim))
-
+ # 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, 2040]
+ 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)
- #self.pfb = gr.fir_filter_ccf(self._decim, self._taps)
+
+ # 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
@@ -71,9 +101,9 @@
d = tb.snk_i.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,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")
@@ -84,7 +114,10 @@
sp1_f.set_xlabel("Frequency (Hz)")
sp1_f.set_ylabel("Power (dBW)")
- t_in = scipy.arange(0, tb._Tmax, tb._Tmax/float(len(d)))
+ 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")
@@ -100,9 +133,9 @@
sp2_f = fig2.add_subplot(2, 1, 1)
d = tb.snk.data()[Ns:Ns+Ne]
- X,freq = sp2_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o,
- window = lambda d: d*winfunc(fftlen),
- visible=False)
+ 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")
@@ -114,12 +147,15 @@
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, tb._Tmax, tb._Tmax/float(x_o.size))
+ 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")
- p2_t = sp2_t.plot(t_o, x_o.imag, "r")
- sp2_t.set_ylim([-1.5, 1.5])
+ 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")
Modified:
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/interpolate.py
===================================================================
---
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/interpolate.py
2009-08-12 03:39:03 UTC (rev 11583)
+++
gnuradio/branches/developers/n4hy/pfb_iir2/gnuradio-examples/python/pfb/interpolate.py
2009-08-12 13:30:13 UTC (rev 11584)
@@ -1,10 +1,30 @@
#!/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 math
import os
import scipy, pylab
from scipy import fftpack
+from pylab import mlab
import time
#print os.getpid()
@@ -14,36 +34,52 @@
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 = 500
+ freq2 = 200
- self._N = 100000
- self._fs = 2000
- self._Tmax = self._N * (1.0/self._fs)
- self._interp = 8
- self._ainterp = 8
+ # 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)
- self._taps = gr.firdes.low_pass_2(self._interp, self._interp*self._fs,
freq2+50, 150,120,5)
- self._taps2 = gr.firdes.low_pass_2(32.0, 1,
- 0.4 / (32.0), 0.1/(32.0),100,5)
+ # 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)
- tpc = math.ceil(float(len(self._taps)) / float(self._interp))
-
+ # 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)
- self.pfb_ar = blks2.pfb_arb_resampler_ccf(self._ainterp, self._taps2)
+
+ # 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))
@@ -51,7 +87,7 @@
self.connect(self.signal, self.pfb_ar)
self.connect(self.signal, self.snk_i)
- # Create the sink for the interpolated siganl
+ # 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)
@@ -84,9 +120,9 @@
d = tb.snk_i.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,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")
@@ -118,9 +154,9 @@
sp2_f = fig2.add_subplot(2, 1, 1)
d = tb.snk1.data()[Ns:Ns+(tb._interp*Ne)]
- X,freq = sp2_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs,
- window = lambda d: d*winfunc(fftlen),
- visible=False)
+ 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")
@@ -151,20 +187,15 @@
sp3_f = fig3.add_subplot(2, 1, 1)
d = tb.snk2.data()[Ns:Ns+(tb._interp*Ne)]
- X,freq = sp3_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs,
- window = lambda d: d*winfunc(fftlen),
- visible=False)
+ 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])
- #f_o = scipy.arange(-fs_int/2.0, fs_int/2.0, fs_int/float(len(d)))
- #X_o = 10.0*scipy.log10(fftpack.fftshift(fftpack.fft(d)))
- #p3_f = sp3_f.plot(f_o, X_o)
- #sp3_f.set_ylim([-50.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)")
@@ -173,10 +204,11 @@
Tmax = len(d)*Ts_aint
t_o = scipy.arange(0, Tmax, Ts_aint)
- x_o1 = scipy.array(d)
+ x_o2 = scipy.array(d)
sp3_f = fig3.add_subplot(2, 1, 2)
- p3_f = sp3_f.plot(t_o, x_o1.real, "b-o")
- #p3_f = sp3_f.plot(t_o, x_o.imag, "r-o")
+ 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")
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Commit-gnuradio] r11584 - in gnuradio/branches/developers/n4hy/pfb_iir2: gnuradio-core/src/lib/filter gnuradio-examples/python/pfb,
n4hy <=