diff options
Diffstat (limited to 'gnuradio-core/src')
57 files changed, 2845 insertions, 800 deletions
diff --git a/gnuradio-core/src/lib/filter/Makefile.am b/gnuradio-core/src/lib/filter/Makefile.am index d5afd571b..23c1dadc3 100644 --- a/gnuradio-core/src/lib/filter/Makefile.am +++ b/gnuradio-core/src/lib/filter/Makefile.am @@ -184,6 +184,8 @@ libfilter_la_common_SOURCES = \ $(GENERATED_CC) \ gr_adaptive_fir_ccf.cc \ gr_cma_equalizer_cc.cc \ + gri_fft_filter_fff_generic.cc \ + gri_fft_filter_ccc_generic.cc \ gr_fft_filter_ccc.cc \ gr_fft_filter_fff.cc \ gr_goertzel_fc.cc \ @@ -206,7 +208,8 @@ libfilter_la_common_SOURCES = \ gr_pfb_decimator_ccf.cc \ gr_pfb_interpolator_ccf.cc \ gr_pfb_arb_resampler_ccf.cc \ - gr_pfb_clock_sync_ccf.cc + gr_pfb_clock_sync_ccf.cc \ + gr_pfb_clock_sync_fff.cc libfilter_qa_la_common_SOURCES = \ qa_filter.cc \ @@ -258,6 +261,8 @@ grinclude_HEADERS = \ gr_altivec.h \ gr_cma_equalizer_cc.h \ gr_cpu.h \ + gri_fft_filter_fff_generic.h \ + gri_fft_filter_ccc_generic.h \ gr_fft_filter_ccc.h \ gr_fft_filter_fff.h \ gr_filter_delay_fc.h \ @@ -286,7 +291,8 @@ grinclude_HEADERS = \ gr_pfb_decimator_ccf.h \ gr_pfb_interpolator_ccf.h \ gr_pfb_arb_resampler_ccf.h \ - gr_pfb_clock_sync_ccf.h + gr_pfb_clock_sync_ccf.h \ + gr_pfb_clock_sync_fff.h noinst_HEADERS = \ assembly.h \ @@ -342,6 +348,7 @@ swiginclude_HEADERS = \ gr_pfb_interpolator_ccf.i \ gr_pfb_arb_resampler_ccf.i \ gr_pfb_clock_sync_ccf.i \ + gr_pfb_clock_sync_fff.i \ $(GENERATED_I) endif diff --git a/gnuradio-core/src/lib/filter/filter.i b/gnuradio-core/src/lib/filter/filter.i index 91f55c514..bdfb8fa8d 100644 --- a/gnuradio-core/src/lib/filter/filter.i +++ b/gnuradio-core/src/lib/filter/filter.i @@ -37,6 +37,7 @@ #include <gr_pfb_interpolator_ccf.h> #include <gr_pfb_arb_resampler_ccf.h> #include <gr_pfb_clock_sync_ccf.h> +#include <gr_pfb_clock_sync_fff.h> %} %include "gr_iir_filter_ffd.i" @@ -58,5 +59,6 @@ %include "gr_pfb_interpolator_ccf.i" %include "gr_pfb_arb_resampler_ccf.i" %include "gr_pfb_clock_sync_ccf.i" +%include "gr_pfb_clock_sync_fff.i" %include "filter_generated.i" diff --git a/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.cc b/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.cc index 3dd40d56d..4540c6e4a 100644 --- a/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.cc +++ b/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.cc @@ -30,6 +30,8 @@ #endif #include <gr_fft_filter_ccc.h> +//#include <gri_fft_filter_ccc_sse.h> +#include <gri_fft_filter_ccc_generic.h> #include <gr_io_signature.h> #include <gri_fft.h> #include <math.h> @@ -52,32 +54,23 @@ gr_fft_filter_ccc::gr_fft_filter_ccc (int decimation, const std::vector<gr_compl gr_make_io_signature (1, 1, sizeof (gr_complex)), gr_make_io_signature (1, 1, sizeof (gr_complex)), decimation), - d_fftsize(-1), d_fwdfft(0), d_invfft(0), d_updated(false) + d_updated(false) { - // if (decimation != 1) - // throw std::invalid_argument("gr_fft_filter_ccc: decimation must be 1"); - set_history(1); - actual_set_taps(taps); +#if 1 // don't enable the sse version until handling it is worked out + d_filter = new gri_fft_filter_ccc_generic(decimation, taps); +#else + d_filter = new gri_fft_filter_ccc_sse(decimation, taps); +#endif + d_nsamples = d_filter->set_taps(taps); + set_output_multiple(d_nsamples); } gr_fft_filter_ccc::~gr_fft_filter_ccc () { - delete d_fwdfft; - delete d_invfft; + delete d_filter; } -#if 0 -static void -print_vector_complex(const std::string label, const std::vector<gr_complex> &x) -{ - std::cout << label; - for (unsigned i = 0; i < x.size(); i++) - std::cout << x[i] << " "; - std::cout << "\n"; -} -#endif - void gr_fft_filter_ccc::set_taps (const std::vector<gr_complex> &taps) { @@ -85,130 +78,26 @@ gr_fft_filter_ccc::set_taps (const std::vector<gr_complex> &taps) d_updated = true; } -/* - * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps - */ -void -gr_fft_filter_ccc::actual_set_taps (const std::vector<gr_complex> &taps) -{ - int i = 0; - compute_sizes(taps.size()); - - d_tail.resize(tailsize()); - for (i = 0; i < tailsize(); i++) - d_tail[i] = 0; - - gr_complex *in = d_fwdfft->get_inbuf(); - gr_complex *out = d_fwdfft->get_outbuf(); - - float scale = 1.0 / d_fftsize; - - // Compute forward xform of taps. - // Copy taps into first ntaps slots, then pad with zeros - for (i = 0; i < d_ntaps; i++) - in[i] = taps[i] * scale; - - for (; i < d_fftsize; i++) - in[i] = 0; - - d_fwdfft->execute(); // do the xform - - // now copy output to d_xformed_taps - for (i = 0; i < d_fftsize; i++) - d_xformed_taps[i] = out[i]; - - //print_vector_complex("transformed taps:", d_xformed_taps); -} - -// determine and set d_ntaps, d_nsamples, d_fftsize - -void -gr_fft_filter_ccc::compute_sizes(int ntaps) -{ - int old_fftsize = d_fftsize; - d_ntaps = ntaps; - d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2)))); - d_nsamples = d_fftsize - d_ntaps + 1; - - if (0) - fprintf(stderr, "gr_fft_filter: ntaps = %d, fftsize = %d, nsamples = %d\n", - d_ntaps, d_fftsize, d_nsamples); - - assert(d_fftsize == d_ntaps + d_nsamples -1 ); - - if (d_fftsize != old_fftsize){ // compute new plans - delete d_fwdfft; - delete d_invfft; - d_fwdfft = new gri_fft_complex(d_fftsize, true); - d_invfft = new gri_fft_complex(d_fftsize, false); - d_xformed_taps.resize(d_fftsize); - } - - set_output_multiple(d_nsamples); -} - int gr_fft_filter_ccc::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]; + const gr_complex *in = (const gr_complex *) input_items[0]; gr_complex *out = (gr_complex *) output_items[0]; if (d_updated){ - actual_set_taps(d_new_taps); + d_nsamples = d_filter->set_taps(d_new_taps); d_updated = false; + set_output_multiple(d_nsamples); return 0; // output multiple may have changed } assert(noutput_items % d_nsamples == 0); - int dec_ctr = 0; - int j = 0; - int ninput_items = noutput_items * decimation(); - - for (int i = 0; i < ninput_items; i += d_nsamples){ - - memcpy(d_fwdfft->get_inbuf(), &in[i], d_nsamples * sizeof(gr_complex)); - - for (j = d_nsamples; j < d_fftsize; j++) - d_fwdfft->get_inbuf()[j] = 0; - - d_fwdfft->execute(); // compute fwd xform - - gr_complex *a = d_fwdfft->get_outbuf(); - gr_complex *b = &d_xformed_taps[0]; - gr_complex *c = d_invfft->get_inbuf(); - - for (j = 0; j < d_fftsize; j++) // filter in the freq domain - c[j] = a[j] * b[j]; - - d_invfft->execute(); // compute inv xform - - // add in the overlapping tail - - for (j = 0; j < tailsize(); j++) - d_invfft->get_outbuf()[j] += d_tail[j]; - - // copy nsamples to output - - //memcpy(out, d_invfft->get_outbuf(), d_nsamples * sizeof(gr_complex)); - //out += d_nsamples; - - j = dec_ctr; - while (j < d_nsamples) { - *out++ = d_invfft->get_outbuf()[j]; - j += decimation(); - } - dec_ctr = (j - d_nsamples); - - // stash the tail - memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples, - tailsize() * sizeof(gr_complex)); - } + d_filter->filter(noutput_items, in, out); - assert((out - (gr_complex *) output_items[0]) == noutput_items); - assert(dec_ctr == 0); + //assert((out - (gr_complex *) output_items[0]) == noutput_items); return noutput_items; } diff --git a/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.h b/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.h index c5363dcbb..68b19e775 100644 --- a/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.h +++ b/gnuradio-core/src/lib/filter/gr_fft_filter_ccc.h @@ -28,8 +28,8 @@ class gr_fft_filter_ccc; typedef boost::shared_ptr<gr_fft_filter_ccc> gr_fft_filter_ccc_sptr; gr_fft_filter_ccc_sptr gr_make_fft_filter_ccc (int decimation, const std::vector<gr_complex> &taps); -class gr_fir_ccc; -class gri_fft_complex; +//class gri_fft_filter_ccc_sse; +class gri_fft_filter_ccc_generic; /*! * \brief Fast FFT filter with gr_complex input, gr_complex output and gr_complex taps @@ -40,15 +40,14 @@ class gr_fft_filter_ccc : public gr_sync_decimator private: friend gr_fft_filter_ccc_sptr gr_make_fft_filter_ccc (int decimation, const std::vector<gr_complex> &taps); - int d_ntaps; int d_nsamples; - int d_fftsize; // fftsize = ntaps + nsamples - 1 - gri_fft_complex *d_fwdfft; // forward "plan" - gri_fft_complex *d_invfft; // inverse "plan" - std::vector<gr_complex> d_tail; // state carried between blocks for overlap-add - std::vector<gr_complex> d_xformed_taps; // Fourier xformed taps - std::vector<gr_complex> d_new_taps; bool d_updated; +#if 1 // don't enable the sse version until handling it is worked out + gri_fft_filter_ccc_generic *d_filter; +#else + gri_fft_filter_ccc_sse *d_filter; +#endif + std::vector<gr_complex> d_new_taps; /*! * Construct a FFT filter with the given taps @@ -58,10 +57,6 @@ class gr_fft_filter_ccc : public gr_sync_decimator */ gr_fft_filter_ccc (int decimation, const std::vector<gr_complex> &taps); - void compute_sizes(int ntaps); - int tailsize() const { return d_ntaps - 1; } - void actual_set_taps (const std::vector<gr_complex> &taps); - public: ~gr_fft_filter_ccc (); diff --git a/gnuradio-core/src/lib/filter/gr_fft_filter_fff.cc b/gnuradio-core/src/lib/filter/gr_fft_filter_fff.cc index 57232f3fb..e8857fe8c 100644 --- a/gnuradio-core/src/lib/filter/gr_fft_filter_fff.cc +++ b/gnuradio-core/src/lib/filter/gr_fft_filter_fff.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2005 Free Software Foundation, Inc. + * Copyright 2005,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -25,13 +25,11 @@ #endif #include <gr_fft_filter_fff.h> +#include <gri_fft_filter_fff_generic.h> +//#include <gri_fft_filter_fff_sse.h> #include <gr_io_signature.h> -#include <gri_fft.h> -#include <math.h> #include <assert.h> #include <stdexcept> -#include <gr_firdes.h> - #include <cstdio> #include <iostream> @@ -48,37 +46,24 @@ gr_fft_filter_fff::gr_fft_filter_fff (int decimation, const std::vector<float> & gr_make_io_signature (1, 1, sizeof (float)), gr_make_io_signature (1, 1, sizeof (float)), decimation), - d_fftsize(-1), d_fwdfft(0), d_invfft(0), d_updated(false) + d_updated(false) { set_history(1); - actual_set_taps(taps); -} - -gr_fft_filter_fff::~gr_fft_filter_fff () -{ - delete d_fwdfft; - delete d_invfft; -} + +#if 1 // don't enable the sse version until handling it is worked out + d_filter = new gri_fft_filter_fff_generic(decimation, taps); +#else + d_filter = new gri_fft_filter_fff_sse(decimation, taps); +#endif -#if 0 -static void -print_vector_complex(const std::string label, const std::vector<gr_complex> &x) -{ - std::cout << label; - for (unsigned i = 0; i < x.size(); i++) - std::cout << x[i] << " "; - std::cout << "\n"; + d_nsamples = d_filter->set_taps(taps); + set_output_multiple(d_nsamples); } -static void -print_vector_float(const std::string label, const std::vector<float> &x) +gr_fft_filter_fff::~gr_fft_filter_fff () { - std::cout << label; - for (unsigned i = 0; i < x.size(); i++) - std::cout << x[i] << " "; - std::cout << "\n"; + delete d_filter; } -#endif void gr_fft_filter_fff::set_taps (const std::vector<float> &taps) @@ -87,68 +72,6 @@ gr_fft_filter_fff::set_taps (const std::vector<float> &taps) d_updated = true; } -/* - * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps - */ -void -gr_fft_filter_fff::actual_set_taps (const std::vector<float> &taps) -{ - int i = 0; - compute_sizes(taps.size()); - - d_tail.resize(tailsize()); - for (i = 0; i < tailsize(); i++) - d_tail[i] = 0; - - float *in = d_fwdfft->get_inbuf(); - gr_complex *out = d_fwdfft->get_outbuf(); - - float scale = 1.0 / d_fftsize; - - // Compute forward xform of taps. - // Copy taps into first ntaps slots, then pad with zeros - for (i = 0; i < d_ntaps; i++) - in[i] = taps[i] * scale; - - for (; i < d_fftsize; i++) - in[i] = 0; - - d_fwdfft->execute(); // do the xform - - // now copy output to d_xformed_taps - for (i = 0; i < d_fftsize/2+1; i++) - d_xformed_taps[i] = out[i]; - - //print_vector_complex("transformed taps:", d_xformed_taps); -} - -// determine and set d_ntaps, d_nsamples, d_fftsize - -void -gr_fft_filter_fff::compute_sizes(int ntaps) -{ - int old_fftsize = d_fftsize; - d_ntaps = ntaps; - d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2)))); - d_nsamples = d_fftsize - d_ntaps + 1; - - if (0) - fprintf(stderr, "gr_fft_filter: ntaps = %d, fftsize = %d, nsamples = %d\n", - d_ntaps, d_fftsize, d_nsamples); - - assert(d_fftsize == d_ntaps + d_nsamples -1 ); - - if (d_fftsize != old_fftsize){ // compute new plans - delete d_fwdfft; - delete d_invfft; - d_fwdfft = new gri_fft_real_fwd(d_fftsize); - d_invfft = new gri_fft_real_rev(d_fftsize); - d_xformed_taps.resize(d_fftsize/2+1); - } - - set_output_multiple(d_nsamples); -} - int gr_fft_filter_fff::work (int noutput_items, gr_vector_const_void_star &input_items, @@ -158,59 +81,17 @@ gr_fft_filter_fff::work (int noutput_items, float *out = (float *) output_items[0]; if (d_updated){ - actual_set_taps(d_new_taps); + d_nsamples = d_filter->set_taps(d_new_taps); d_updated = false; + set_output_multiple(d_nsamples); return 0; // output multiple may have changed } assert(noutput_items % d_nsamples == 0); + + d_filter->filter(noutput_items, in, out); - int dec_ctr = 0; - int j = 0; - int ninput_items = noutput_items * decimation(); - - for (int i = 0; i < ninput_items; i += d_nsamples){ - - memcpy(d_fwdfft->get_inbuf(), &in[i], d_nsamples * sizeof(float)); - - for (j = d_nsamples; j < d_fftsize; j++) - d_fwdfft->get_inbuf()[j] = 0; - - d_fwdfft->execute(); // compute fwd xform - - gr_complex *a = d_fwdfft->get_outbuf(); - gr_complex *b = &d_xformed_taps[0]; - gr_complex *c = d_invfft->get_inbuf(); - - for (j = 0; j < d_fftsize/2+1; j++) // filter in the freq domain - c[j] = a[j] * b[j]; - - d_invfft->execute(); // compute inv xform - - // add in the overlapping tail - - for (j = 0; j < tailsize(); j++) - d_invfft->get_outbuf()[j] += d_tail[j]; - - // copy nsamples to output - - //memcpy(out, d_invfft->get_outbuf(), d_nsamples * sizeof(float)); - //out += d_nsamples; - - j = dec_ctr; - while (j < d_nsamples) { - *out++ = d_invfft->get_outbuf()[j]; - j += decimation(); - } - dec_ctr = (j - d_nsamples); - - // stash the tail - memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples, - tailsize() * sizeof(float)); - } - - assert((out - (float *) output_items[0]) == noutput_items); - assert(dec_ctr == 0); + //assert((out - (float *) output_items[0]) == noutput_items); return noutput_items; } diff --git a/gnuradio-core/src/lib/filter/gr_fft_filter_fff.h b/gnuradio-core/src/lib/filter/gr_fft_filter_fff.h index b26361107..6eaa21500 100644 --- a/gnuradio-core/src/lib/filter/gr_fft_filter_fff.h +++ b/gnuradio-core/src/lib/filter/gr_fft_filter_fff.h @@ -28,9 +28,8 @@ class gr_fft_filter_fff; typedef boost::shared_ptr<gr_fft_filter_fff> gr_fft_filter_fff_sptr; gr_fft_filter_fff_sptr gr_make_fft_filter_fff (int decimation, const std::vector<float> &taps); -class gr_fir_fff; -class gri_fft_real_fwd; -class gri_fft_real_rev; +class gri_fft_filter_fff_generic; +//class gri_fft_filter_fff_sse; /*! * \brief Fast FFT filter with float input, float output and float taps @@ -41,15 +40,14 @@ class gr_fft_filter_fff : public gr_sync_decimator private: friend gr_fft_filter_fff_sptr gr_make_fft_filter_fff (int decimation, const std::vector<float> &taps); - int d_ntaps; int d_nsamples; - int d_fftsize; // fftsize = ntaps + nsamples - 1 - gri_fft_real_fwd *d_fwdfft; // forward "plan" - gri_fft_real_rev *d_invfft; // inverse "plan" - std::vector<float> d_tail; // state carried between blocks for overlap-add - std::vector<gr_complex> d_xformed_taps; // Fourier xformed taps - std::vector<float> d_new_taps; bool d_updated; +#if 1 // don't enable the sse version until handling it is worked out + gri_fft_filter_fff_generic *d_filter; +#else + gri_fft_filter_fff_sse *d_filter; +#endif + std::vector<float> d_new_taps; /*! * Construct a FFT filter with the given taps @@ -58,10 +56,6 @@ class gr_fft_filter_fff : public gr_sync_decimator * \param taps float filter taps */ gr_fft_filter_fff (int decimation, const std::vector<float> &taps); - - void compute_sizes(int ntaps); - int tailsize() const { return d_ntaps - 1; } - void actual_set_taps (const std::vector<float> &taps); public: ~gr_fft_filter_fff (); 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 index 8971d3d39..5a6e753ab 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc +++ b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009 Free Software Foundation, Inc. + * Copyright 2009,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -47,6 +47,8 @@ gr_pfb_arb_resampler_ccf::gr_pfb_arb_resampler_ccf (float rate, gr_make_io_signature (1, 1, sizeof(gr_complex))), d_updated (false) { + d_acc = 0; // start accumulator at 0 + /* The number of filters is specified by the user as the filter size; this is also the interpolation rate of the filter. We use it and the rate provided to determine the decimation rate. This acts as a @@ -59,22 +61,26 @@ gr_pfb_arb_resampler_ccf::gr_pfb_arb_resampler_ccf (float rate, 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_start_index = 0; d_filters = std::vector<gr_fir_ccf*>(d_int_rate); + d_diff_filters = std::vector<gr_fir_ccf*>(d_int_rate); // Create an FIR filter for each channel and zero out the taps std::vector<float> vtaps(0, d_int_rate); - for(unsigned int i = 0; i < d_int_rate; i++) { + for(int i = 0; i < d_int_rate; i++) { d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps); + d_diff_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps); } // Now, actually set the filters' taps - set_taps(taps); + std::vector<float> dtaps; + create_diff_taps(taps, dtaps); + create_taps(taps, d_taps, d_filters); + create_taps(dtaps, d_dtaps, d_diff_filters); } gr_pfb_arb_resampler_ccf::~gr_pfb_arb_resampler_ccf () @@ -85,20 +91,22 @@ gr_pfb_arb_resampler_ccf::~gr_pfb_arb_resampler_ccf () } void -gr_pfb_arb_resampler_ccf::set_taps (const std::vector<float> &taps) +gr_pfb_arb_resampler_ccf::create_taps (const std::vector<float> &newtaps, + std::vector< std::vector<float> > &ourtaps, + std::vector<gr_fir_ccf*> &ourfilter) { - unsigned int i,j; + int i,j; - unsigned int ntaps = taps.size(); + unsigned int ntaps = newtaps.size(); d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate); // Create d_numchan vectors to store each channel's taps - d_taps.resize(d_int_rate); - + ourtaps.resize(d_int_rate); + // Make a vector of the taps plus fill it out with 0's to fill // each polyphase filter with exactly d_taps_per_filter std::vector<float> tmp_taps; - tmp_taps = taps; + tmp_taps = newtaps; while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) { tmp_taps.push_back(0.0); } @@ -106,22 +114,37 @@ gr_pfb_arb_resampler_ccf::set_taps (const std::vector<float> &taps) // 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); + ourtaps[d_int_rate-1-i] = std::vector<float>(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 + ourtaps[d_int_rate - 1 - i][j] = tmp_taps[i + j*d_int_rate]; } // Build a filter for each channel and add it's taps to it - d_filters[i]->set_taps(d_taps[i]); + ourfilter[i]->set_taps(ourtaps[d_int_rate-1-i]); } // Set the history to ensure enough input items for each filter - set_history (d_taps_per_filter); + set_history (d_taps_per_filter + 1); d_updated = true; } void +gr_pfb_arb_resampler_ccf::create_diff_taps(const std::vector<float> &newtaps, + std::vector<float> &difftaps) +{ + // Calculate the differential taps (derivative filter) by taking the difference + // between two taps. Duplicate the last one to make both filters the same length. + float tap; + difftaps.clear(); + for(unsigned int i = 0; i < newtaps.size()-1; i++) { + tap = newtaps[i+1] - newtaps[i]; + difftaps.push_back(tap); + } + difftaps.push_back(tap); +} + +void gr_pfb_arb_resampler_ccf::print_taps() { unsigned int i, j; @@ -148,7 +171,7 @@ gr_pfb_arb_resampler_ccf::general_work (int noutput_items, return 0; // history requirements may have changed. } - int i = 0, j, count = 0; + int i = 0, j, count = d_start_index; gr_complex o0, o1; // Restore the last filter position @@ -159,39 +182,30 @@ gr_pfb_arb_resampler_ccf::general_work (int noutput_items, // start j by wrapping around mod the number of channels while((j < d_int_rate) && (i < noutput_items)) { - // Take the current filter output + // Take the current filter and derivative filter output o0 = d_filters[j]->filter(&in[count]); + o1 = d_diff_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 + out[i] = o0 + o1*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. + + // Adjust accumulator and index into filterbank d_acc += d_flt_rate; j += d_dec_rate + (int)floor(d_acc); d_acc = fmodf(d_acc, 1.0); } if(i < noutput_items) { // keep state for next entry - count++; // we have fully consumed another input + float ss = (int)(j / d_int_rate); // number of items to skip ahead by + count += ss; // we have fully consumed another input j = j % d_int_rate; // roll filter around } } - // Store the current filter position + // Store the current filter position and start of next sample d_last_filter = j; + d_start_index = std::max(0, count - ninput_items[0]); - consume_each(count); + // consume all we've processed but no more than we can + consume_each(std::min(count, ninput_items[0])); 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 index d4c886ec3..cf5a79d4e 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h +++ b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.h @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009 Free Software Foundation, Inc. + * Copyright 2009,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -112,12 +112,15 @@ class gr_pfb_arb_resampler_ccf : public gr_block unsigned int filter_size); std::vector<gr_fir_ccf*> d_filters; + std::vector<gr_fir_ccf*> d_diff_filters; std::vector< std::vector<float> > d_taps; + std::vector< std::vector<float> > d_dtaps; unsigned int d_int_rate; // the number of filters (interpolation rate) unsigned int d_dec_rate; // the stride through the filters (decimation rate) float d_flt_rate; // residual rate for the linear interpolation float d_acc; unsigned int d_last_filter; + int d_start_index; unsigned int d_taps_per_filter; bool d_updated; @@ -133,16 +136,26 @@ class gr_pfb_arb_resampler_ccf : public gr_block gr_pfb_arb_resampler_ccf (float rate, const std::vector<float> &taps, unsigned int filter_size); - -public: - ~gr_pfb_arb_resampler_ccf (); - + + void create_diff_taps(const std::vector<float> &newtaps, + std::vector<float> &difftaps); + /*! * Resets the filterbank's filter taps with the new prototype filter - * \param taps (vector/list of floats) The prototype filter to populate the filterbank. The taps - * should be generated at the interpolated sampling rate. + * \param newtaps (vector of floats) The prototype filter to populate the filterbank. + * The taps should be generated at the interpolated sampling rate. + * \param ourtaps (vector of floats) Reference to our internal member of holding the taps. + * \param ourfilter (vector of filters) Reference to our internal filter to set the taps for. */ - void set_taps (const std::vector<float> &taps); + void create_taps (const std::vector<float> &newtaps, + std::vector< std::vector<float> > &ourtaps, + std::vector<gr_fir_ccf*> &ourfilter); + + +public: + ~gr_pfb_arb_resampler_ccf (); + + // FIXME: See about a set_taps function during runtime. /*! * Print all of the filterbank taps to screen. 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 index e365e0314..4f07af861 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i +++ b/gnuradio-core/src/lib/filter/gr_pfb_arb_resampler_ccf.i @@ -36,6 +36,6 @@ class gr_pfb_arb_resampler_ccf : public gr_block public: ~gr_pfb_arb_resampler_ccf (); - void set_taps (const std::vector<float> &taps); + //void set_taps (const std::vector<float> &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 index 7e34551c8..5fda47880 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc +++ b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009 Free Software Foundation, Inc. + * Copyright 2009,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -33,20 +33,33 @@ #include <cstring> gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, - const std::vector<float> &taps) + const std::vector<float> &taps, + float oversample_rate) { - return gr_pfb_channelizer_ccf_sptr (new gr_pfb_channelizer_ccf (numchans, taps)); + return gr_pfb_channelizer_ccf_sptr (new gr_pfb_channelizer_ccf (numchans, taps, + oversample_rate)); } gr_pfb_channelizer_ccf::gr_pfb_channelizer_ccf (unsigned int numchans, - const std::vector<float> &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) + const std::vector<float> &taps, + float oversample_rate) + : gr_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_oversample_rate(oversample_rate) { - d_numchans = numchans; + // The over sampling rate must be rationally related to the number of channels + // in that it must be N/i for i in [1,N], which gives an outputsample rate + // of [fs/N, fs] where fs is the input sample rate. + // This tests the specified input sample rate to see if it conforms to this + // requirement within a few significant figures. + double intp = 0; + double x = (10000.0*rint(numchans / oversample_rate)) / 10000.0; + double fltp = modf(numchans / oversample_rate, &intp); + if(fltp != 0.0) + throw std::invalid_argument("gr_pfb_channelizer: oversample rate must be N/i for i in [1, N]"); + d_filters = std::vector<gr_fir_ccf*>(d_numchans); // Create an FIR filter for each channel and zero out the taps @@ -60,10 +73,28 @@ gr_pfb_channelizer_ccf::gr_pfb_channelizer_ccf (unsigned int numchans, // Create the FFT to handle the output de-spinning of the channels d_fft = new gri_fft_complex (d_numchans, false); + + // Although the filters change, we use this look up table + // to set the index of the FFT input buffer, which equivalently + // performs the FFT shift operation on every other turn. + d_rate_ratio = (int)rintf(d_numchans / d_oversample_rate); + d_idxlut = new int[d_numchans]; + for(unsigned int i = 0; i < d_numchans; i++) { + d_idxlut[i] = d_numchans - ((i + d_rate_ratio) % d_numchans) - 1; + } + + // Calculate the number of filtering rounds to do to evenly + // align the input vectors with the output channels + d_output_multiple = 1; + while((d_output_multiple * d_rate_ratio) % d_numchans != 0) + d_output_multiple++; + set_output_multiple(d_output_multiple); } gr_pfb_channelizer_ccf::~gr_pfb_channelizer_ccf () { + delete [] d_idxlut; + for(unsigned int i = 0; i < d_numchans; i++) { delete d_filters[i]; } @@ -101,7 +132,7 @@ gr_pfb_channelizer_ccf::set_taps (const std::vector<float> &taps) } // Set the history to ensure enough input items for each filter - set_history (d_taps_per_filter); + set_history (d_taps_per_filter+1); d_updated = true; } @@ -121,9 +152,10 @@ gr_pfb_channelizer_ccf::print_taps() int -gr_pfb_channelizer_ccf::work (int noutput_items, - gr_vector_const_void_star &input_items, - gr_vector_void_star &output_items) +gr_pfb_channelizer_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]; @@ -133,20 +165,35 @@ gr_pfb_channelizer_ccf::work (int noutput_items, 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]; + int n=1, i=-1, j=0, last; + int toconsume = (int)rintf(noutput_items/d_oversample_rate); + while(n <= toconsume) { + j = 0; + i = (i + d_rate_ratio) % d_numchans; + last = i; + while(i >= 0) { + in = (gr_complex*)input_items[j]; + d_fft->get_inbuf()[d_idxlut[j]] = d_filters[i]->filter(&in[n]); + j++; + i--; + } - // Filter current input stream from bottom filter to top - d_fft->get_inbuf()[j] = d_filters[j]->filter(&in[i]); + i = d_numchans-1; + while(i > last) { + in = (gr_complex*)input_items[j]; + d_fft->get_inbuf()[d_idxlut[j]] = d_filters[i]->filter(&in[n-1]); + j++; + i--; } + n += (i+d_rate_ratio) >= (int)d_numchans; + // despin through FFT d_fft->execute(); - memcpy(&out[d_numchans*i], d_fft->get_outbuf(), d_numchans*sizeof(gr_complex)); + memcpy(out, d_fft->get_outbuf(), d_numchans*sizeof(gr_complex)); + out += d_numchans; } - + + consume_each(toconsume); 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 index b2e67e817..751673bc7 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h +++ b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.h @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009 Free Software Foundation, Inc. + * Copyright 2009,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -24,12 +24,13 @@ #ifndef INCLUDED_GR_PFB_CHANNELIZER_CCF_H #define INCLUDED_GR_PFB_CHANNELIZER_CCF_H -#include <gr_sync_block.h> +#include <gr_block.h> class gr_pfb_channelizer_ccf; typedef boost::shared_ptr<gr_pfb_channelizer_ccf> gr_pfb_channelizer_ccf_sptr; gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, - const std::vector<float> &taps); + const std::vector<float> &taps, + float oversample_rate=1); class gr_fir_ccf; class gri_fft_complex; @@ -88,6 +89,19 @@ class gri_fft_complex; * <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 filter output can also be overs ampled. The over sampling rate + * is the ratio of the the actual output sampling rate to the normal + * output sampling rate. It must be rationally related to the number + * of channels as N/i for i in [1,N], which gives an outputsample rate + * of [fs/N, fs] where fs is the input sample rate and N is the number + * of channels. + * + * For example, for 6 channels with fs = 6000 Hz, the normal rate is + * 6000/6 = 1000 Hz. Allowable oversampling rates are 6/6, 6/5, 6/4, + * 6/3, 6/2, and 6/1 where the output sample rate of a 6/1 oversample + * ratio is 6000 Hz, or 6 times the normal 1000 Hz. A rate of 6/5 = 1.2, + * so the output rate would be 1200 Hz. + * * The theory behind this block can be found in Chapter 6 of * the following book. * @@ -96,31 +110,50 @@ class gri_fft_complex; * */ -class gr_pfb_channelizer_ccf : public gr_sync_block +class gr_pfb_channelizer_ccf : public gr_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. + * \param oversample_rate (float) The over sampling rate is the ratio of the the actual + * output sampling rate to the normal output sampling rate. + * It must be rationally related to the number of channels + * as N/i for i in [1,N], which gives an outputsample rate + * of [fs/N, fs] where fs is the input sample rate and N is + * the number of channels. + * + * For example, for 6 channels with fs = 6000 Hz, the normal + * rate is 6000/6 = 1000 Hz. Allowable oversampling rates + * are 6/6, 6/5, 6/4, 6/3, 6/2, and 6/1 where the output + * sample rate of a 6/1 oversample ratio is 6000 Hz, or + * 6 times the normal 1000 Hz. */ friend gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, - const std::vector<float> &taps); + const std::vector<float> &taps, + float oversample_rate); + bool d_updated; + unsigned int d_numchans; + float d_oversample_rate; 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; + gri_fft_complex *d_fft; + int *d_idxlut; + int d_rate_ratio; + int d_output_multiple; /*! * 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. + * \param oversample_rate (float) The output over sampling rate. */ gr_pfb_channelizer_ccf (unsigned int numchans, - const std::vector<float> &taps); + const std::vector<float> &taps, + float oversample_rate); public: ~gr_pfb_channelizer_ccf (); @@ -136,9 +169,10 @@ public: */ void print_taps(); - int work (int noutput_items, - gr_vector_const_void_star &input_items, - gr_vector_void_star &output_items); + 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_channelizer_ccf.i b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.i index 4bef90e22..63e3e0fe6 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.i +++ b/gnuradio-core/src/lib/filter/gr_pfb_channelizer_ccf.i @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009 Free Software Foundation, Inc. + * Copyright 2009,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -23,13 +23,15 @@ GR_SWIG_BLOCK_MAGIC(gr,pfb_channelizer_ccf); gr_pfb_channelizer_ccf_sptr gr_make_pfb_channelizer_ccf (unsigned int numchans, - const std::vector<float> &taps); + const std::vector<float> &taps, + float oversample_rate=1); -class gr_pfb_channelizer_ccf : public gr_sync_block +class gr_pfb_channelizer_ccf : public gr_block { private: gr_pfb_channelizer_ccf (unsigned int numchans, - const std::vector<float> &taps); + const std::vector<float> &taps, + float oversample_rate); public: ~gr_pfb_channelizer_ccf (); diff --git a/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.cc b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.cc index 433b7d613..ff4fb70a3 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.cc +++ b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009 Free Software Foundation, Inc. + * Copyright 2009,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -45,8 +45,8 @@ gr_pfb_clock_sync_ccf_sptr gr_make_pfb_clock_sync_ccf (double sps, float gain, max_rate_deviation)); } -int ios[] = {sizeof(gr_complex), sizeof(float), sizeof(float), sizeof(float)}; -std::vector<int> iosig(ios, ios+sizeof(ios)/sizeof(int)); +static int ios[] = {sizeof(gr_complex), sizeof(float), sizeof(float), sizeof(float)}; +static std::vector<int> iosig(ios, ios+sizeof(ios)/sizeof(int)); gr_pfb_clock_sync_ccf::gr_pfb_clock_sync_ccf (double sps, float gain, const std::vector<float> &taps, unsigned int filter_size, @@ -93,9 +93,16 @@ gr_pfb_clock_sync_ccf::~gr_pfb_clock_sync_ccf () { for(int i = 0; i < d_nfilters; i++) { delete d_filters[i]; + delete d_diff_filters[i]; } } +bool +gr_pfb_clock_sync_ccf::check_topology(int ninputs, int noutputs) +{ + return noutputs == 1 || noutputs == 4; +} + void gr_pfb_clock_sync_ccf::set_taps (const std::vector<float> &newtaps, std::vector< std::vector<float> > &ourtaps, @@ -219,8 +226,8 @@ gr_pfb_clock_sync_ccf::general_work (int noutput_items, gr_complex *in = (gr_complex *) input_items[0]; gr_complex *out = (gr_complex *) output_items[0]; - float *err, *outrate, *outk; - if(output_items.size() > 2) { + float *err = 0, *outrate = 0, *outk = 0; + if(output_items.size() == 4) { err = (float *) output_items[1]; outrate = (float*)output_items[2]; outk = (float*)output_items[3]; @@ -271,7 +278,7 @@ gr_pfb_clock_sync_ccf::general_work (int noutput_items, i++; count += (int)floor(d_sps); - if(output_items.size() > 2) { + if(output_items.size() == 4) { err[i] = error; outrate[i] = d_rate_f; outk[i] = d_k; diff --git a/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.h b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.h index a07192a7f..4e6ef5fc4 100644 --- a/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.h +++ b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.h @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009 Free Software Foundation, Inc. + * Copyright 2009,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -43,6 +43,71 @@ class gr_fir_ccf; * * \ingroup filter_blk * + * This block performs timing synchronization for PAM signals by minimizing the + * derivative of the filtered signal, which in turn maximizes the SNR and + * minimizes ISI. + * + * This approach works by setting up two filterbanks; one filterbanke contains the + * signal's pulse shaping matched filter (such as a root raised cosine filter), + * where each branch of the filterbank contains a different phase of the filter. + * The second filterbank contains the derivatives of the filters in the first + * filterbank. Thinking of this in the time domain, the first filterbank contains + * filters that have a sinc shape to them. We want to align the output signal to + * be sampled at exactly the peak of the sinc shape. The derivative of the sinc + * contains a zero at the maximum point of the sinc (sinc(0) = 1, sinc(0)' = 0). + * Furthermore, the region around the zero point is relatively linear. We make + * use of this fact to generate the error signal. + * + * If the signal out of the derivative filters is d_i[n] for the ith filter, and + * the output of the matched filter is x_i[n], we calculate the error as: + * e[n] = (Re{x_i[n]} * Re{d_i[n]} + Im{x_i[n]} * Im{d_i[n]}) / 2.0 + * This equation averages the error in the real and imaginary parts. There are two + * reasons we multiply by the signal itself. First, if the symbol could be positive + * or negative going, but we want the error term to always tell us to go in the + * same direction depending on which side of the zero point we are on. The sign of + * x_i[n] adjusts the error term to do this. Second, the magnitude of x_i[n] scales + * the error term depending on the symbol's amplitude, so larger signals give us + * a stronger error term because we have more confidence in that symbol's value. + * Using the magnitude of x_i[n] instead of just the sign is especially good for + * signals with low SNR. + * + * The error signal, e[n], gives us a value proportional to how far away from the zero + * point we are in the derivative signal. We want to drive this value to zero, so we + * set up a second order loop. We have two variables for this loop; d_k is the filter + * number in the filterbank we are on and d_rate is the rate which we travel through + * the filters in the steady state. That is, due to the natural clock differences between + * the transmitter and receiver, d_rate represents that difference and would traverse + * the filter phase paths to keep the receiver locked. Thinking of this as a second-order + * PLL, the d_rate is the frequency and d_k is the phase. So we update d_rate and d_k + * using the standard loop equations based on two error signals, d_alpha and d_beta. + * We have these two values set based on each other for a critically damped system, so in + * the block constructor, we just ask for "gain," which is d_alpha while d_beta is + * equal to (gain^2)/4. + * + * The clock sync block needs to know the number of samples per second (sps), because it + * only returns a single point representing the sample. The sps can be any positive real + * number and does not need to be an integer. The filter taps must also be specified. The + * taps are generated by first conceiving of the prototype filter that would be the signal's + * matched filter. Then interpolate this by the number of filters in the filterbank. These + * are then distributed among all of the filters. So if the prototype filter was to have + * 45 taps in it, then each path of the filterbank will also have 45 taps. This is easily + * done by building the filter with the sample rate multiplied by the number of filters + * to use. + * + * The number of filters can also be set and defaults to 32. With 32 filters, you get a + * good enough resolution in the phase to produce very small, almost unnoticeable, ISI. + * Going to 64 filters can reduce this more, but after that there is very little gained + * for the extra complexity. + * + * The initial phase is another settable parameter and refers to the filter path the + * algorithm initially looks at (i.e., d_k starts at init_phase). This value defaults + * to zero, but it might be useful to start at a different phase offset, such as the mid- + * point of the filters. + * + * The final parameter is the max_rate_devitation, which defaults to 1.5. This is how far + * we allow d_rate to swing, positive or negative, from 0. Constraining the rate can help + * keep the algorithm from walking too far away to lock during times when there is no signal. + * */ class gr_pfb_clock_sync_ccf : public gr_block @@ -50,6 +115,14 @@ class gr_pfb_clock_sync_ccf : public gr_block private: /*! * Build the polyphase filterbank timing synchronizer. + * \param sps (double) The number of samples per second in the incoming signal + * \param gain (float) The alpha gain of the control loop; beta = (gain^2)/4 by default. + * \param taps (vector<int>) The filter taps. + * \param filter_size (uint) The number of filters in the filterbank (default = 32). + * \param init_phase (float) The initial phase to look at, or which filter to start + * with (default = 0). + * \param max_rate_deviation (float) Distance from 0 d_rate can get (default = 1.5). + * */ friend gr_pfb_clock_sync_ccf_sptr gr_make_pfb_clock_sync_ccf (double sps, float gain, const std::vector<float> &taps, @@ -96,29 +169,53 @@ public: void set_taps (const std::vector<float> &taps, std::vector< std::vector<float> > &ourtaps, std::vector<gr_fir_ccf*> &ourfilter); + + /*! + * Returns the taps of the matched filter + */ std::vector<float> channel_taps(int channel); + + /*! + * Returns the taps in the derivative filter + */ std::vector<float> diff_channel_taps(int channel); /*! * Print all of the filterbank taps to screen. */ void print_taps(); + + /*! + * Print all of the filterbank taps of the derivative filter to screen. + */ void print_diff_taps(); + /*! + * Set the gain value alpha for the control loop + */ void set_alpha(float alpha) { d_alpha = alpha; } + + /*! + * Set the gain value beta for the control loop + */ void set_beta(float beta) { d_beta = beta; } + /*! + * Set the maximum deviation from 0 d_rate can have + */ void set_max_rate_deviation(float m) { d_max_dev = m; } + bool check_topology(int ninputs, int noutputs); + int general_work (int noutput_items, gr_vector_int &ninput_items, gr_vector_const_void_star &input_items, diff --git a/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.cc b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.cc new file mode 100644 index 000000000..86de3b5a1 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.cc @@ -0,0 +1,288 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009,2010 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 <cstdio> +#include <cmath> + +#include <gr_pfb_clock_sync_fff.h> +#include <gr_fir_fff.h> +#include <gr_fir_util.h> +#include <gr_io_signature.h> +#include <gr_math.h> + +gr_pfb_clock_sync_fff_sptr gr_make_pfb_clock_sync_fff (double sps, float gain, + const std::vector<float> &taps, + unsigned int filter_size, + float init_phase, + float max_rate_deviation) +{ + return gr_pfb_clock_sync_fff_sptr (new gr_pfb_clock_sync_fff (sps, gain, taps, + filter_size, + init_phase, + max_rate_deviation)); +} + +static int ios[] = {sizeof(float), sizeof(float), sizeof(float), sizeof(float)}; +static std::vector<int> iosig(ios, ios+sizeof(ios)/sizeof(int)); +gr_pfb_clock_sync_fff::gr_pfb_clock_sync_fff (double sps, float gain, + const std::vector<float> &taps, + unsigned int filter_size, + float init_phase, + float max_rate_deviation) + : gr_block ("pfb_clock_sync_fff", + gr_make_io_signature (1, 1, sizeof(float)), + gr_make_io_signaturev (1, 4, iosig)), + d_updated (false), d_nfilters(filter_size), + d_max_dev(max_rate_deviation) +{ + d_nfilters = filter_size; + d_sps = floor(sps); + + // Store the last filter between calls to work + // The accumulator keeps track of overflow to increment the stride correctly. + // set it here to the fractional difference based on the initial phaes + set_alpha(gain); + set_beta(0.25*gain*gain); + d_k = init_phase; + d_rate = (sps-floor(sps))*(double)d_nfilters; + d_rate_i = (int)floor(d_rate); + d_rate_f = d_rate - (float)d_rate_i; + d_filtnum = (int)floor(d_k); + + d_filters = std::vector<gr_fir_fff*>(d_nfilters); + d_diff_filters = std::vector<gr_fir_fff*>(d_nfilters); + + // Create an FIR filter for each channel and zero out the taps + std::vector<float> vtaps(0, d_nfilters); + for(int i = 0; i < d_nfilters; i++) { + d_filters[i] = gr_fir_util::create_gr_fir_fff(vtaps); + d_diff_filters[i] = gr_fir_util::create_gr_fir_fff(vtaps); + } + + // Now, actually set the filters' taps + std::vector<float> dtaps; + create_diff_taps(taps, dtaps); + set_taps(taps, d_taps, d_filters); + set_taps(dtaps, d_dtaps, d_diff_filters); +} + +gr_pfb_clock_sync_fff::~gr_pfb_clock_sync_fff () +{ + for(int i = 0; i < d_nfilters; i++) { + delete d_filters[i]; + delete d_diff_filters[i]; + } +} + +bool +gr_pfb_clock_sync_fff::check_topology(int ninputs, int noutputs) +{ + return noutputs == 1 || noutputs == 4; +} + +void +gr_pfb_clock_sync_fff::set_taps (const std::vector<float> &newtaps, + std::vector< std::vector<float> > &ourtaps, + std::vector<gr_fir_fff*> &ourfilter) +{ + int i,j; + + unsigned int ntaps = newtaps.size(); + d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_nfilters); + + // Create d_numchan vectors to store each channel's taps + ourtaps.resize(d_nfilters); + + // Make a vector of the taps plus fill it out with 0's to fill + // each polyphase filter with exactly d_taps_per_filter + std::vector<float> tmp_taps; + tmp_taps = newtaps; + while((float)(tmp_taps.size()) < d_nfilters*d_taps_per_filter) { + tmp_taps.push_back(0.0); + } + + // Partition the filter + for(i = 0; i < d_nfilters; i++) { + // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out + ourtaps[d_nfilters-1-i] = std::vector<float>(d_taps_per_filter, 0); + for(j = 0; j < d_taps_per_filter; j++) { + ourtaps[d_nfilters - 1 - i][j] = tmp_taps[i + j*d_nfilters]; + } + + // Build a filter for each channel and add it's taps to it + ourfilter[i]->set_taps(ourtaps[d_nfilters-1-i]); + } + + // Set the history to ensure enough input items for each filter + set_history (d_taps_per_filter + d_sps); + + d_updated = true; +} + +void +gr_pfb_clock_sync_fff::create_diff_taps(const std::vector<float> &newtaps, + std::vector<float> &difftaps) +{ + float maxtap = 1e-20; + difftaps.clear(); + difftaps.push_back(0); //newtaps[0]); + for(unsigned int i = 1; i < newtaps.size()-1; i++) { + float tap = newtaps[i+1] - newtaps[i-1]; + difftaps.push_back(tap); + if(tap > maxtap) { + maxtap = tap; + } + } + difftaps.push_back(0);//-newtaps[newtaps.size()-1]); + + // Scale the differential taps; helps scale error term to better update state + // FIXME: should this be scaled this way or use the same gain as the taps? + for(unsigned int i = 0; i < difftaps.size(); i++) { + difftaps[i] /= maxtap; + } +} + +void +gr_pfb_clock_sync_fff::print_taps() +{ + int i, j; + printf("[ "); + for(i = 0; i < d_nfilters; i++) { + printf("[%.4e, ", d_taps[i][0]); + for(j = 1; j < d_taps_per_filter-1; j++) { + printf("%.4e,", d_taps[i][j]); + } + printf("%.4e],", d_taps[i][j]); + } + printf(" ]\n"); +} + +void +gr_pfb_clock_sync_fff::print_diff_taps() +{ + int i, j; + printf("[ "); + for(i = 0; i < d_nfilters; i++) { + printf("[%.4e, ", d_dtaps[i][0]); + for(j = 1; j < d_taps_per_filter-1; j++) { + printf("%.4e,", d_dtaps[i][j]); + } + printf("%.4e],", d_dtaps[i][j]); + } + printf(" ]\n"); +} + + +std::vector<float> +gr_pfb_clock_sync_fff::channel_taps(int channel) +{ + std::vector<float> taps; + for(int i = 0; i < d_taps_per_filter; i++) { + taps.push_back(d_taps[channel][i]); + } + return taps; +} + +std::vector<float> +gr_pfb_clock_sync_fff::diff_channel_taps(int channel) +{ + std::vector<float> taps; + for(int i = 0; i < d_taps_per_filter; i++) { + taps.push_back(d_dtaps[channel][i]); + } + return taps; +} + + +int +gr_pfb_clock_sync_fff::general_work (int noutput_items, + gr_vector_int &ninput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + float *in = (float *) input_items[0]; + float *out = (float *) output_items[0]; + + float *err = 0, *outrate = 0, *outk = 0; + if(output_items.size() == 4) { + err = (float *) output_items[1]; + outrate = (float*)output_items[2]; + outk = (float*)output_items[3]; + } + + if (d_updated) { + d_updated = false; + return 0; // history requirements may have changed. + } + + // We need this many to process one output + int nrequired = ninput_items[0] - d_taps_per_filter; + + int i = 0, count = 0; + float error; + + // produce output as long as we can and there are enough input samples + while((i < noutput_items) && (count < nrequired)) { + d_filtnum = (int)floor(d_k); + + // Keep the current filter number in [0, d_nfilters] + // If we've run beyond the last filter, wrap around and go to next sample + // If we've go below 0, wrap around and go to previous sample + while(d_filtnum >= d_nfilters) { + d_k -= d_nfilters; + d_filtnum -= d_nfilters; + count += 1; + } + while(d_filtnum < 0) { + d_k += d_nfilters; + d_filtnum += d_nfilters; + count -= 1; + } + + out[i] = d_filters[d_filtnum]->filter(&in[count]); + float diff = d_diff_filters[d_filtnum]->filter(&in[count]); + error = out[i] * diff; + + // Run the control loop to update the current phase (k) and tracking rate + d_k = d_k + d_alpha*error + d_rate_i + d_rate_f; + d_rate_f = d_rate_f + d_beta*error; + + // Keep our rate within a good range + d_rate_f = gr_branchless_clip(d_rate_f, d_max_dev); + + i++; + count += (int)floor(d_sps); + + if(output_items.size() == 4) { + err[i] = error; + outrate[i] = d_rate_f; + outk[i] = d_k; + } + } + consume_each(count); + + return i; +} diff --git a/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.h b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.h new file mode 100644 index 000000000..fa1279a7c --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.h @@ -0,0 +1,225 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009,2010 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_CLOCK_SYNC_FFF_H +#define INCLUDED_GR_PFB_CLOCK_SYNC_FFF_H + +#include <gr_block.h> + +class gr_pfb_clock_sync_fff; +typedef boost::shared_ptr<gr_pfb_clock_sync_fff> gr_pfb_clock_sync_fff_sptr; +gr_pfb_clock_sync_fff_sptr gr_make_pfb_clock_sync_fff (double sps, float gain, + const std::vector<float> &taps, + unsigned int filter_size=32, + float init_phase=0, + float max_rate_deviation=1.5); + +class gr_fir_fff; + +/*! + * \class gr_pfb_clock_sync_fff + * + * \brief Timing synchronizer using polyphase filterbanks + * + * \ingroup filter_blk + * + * This block performs timing synchronization for PAM signals by minimizing the + * derivative of the filtered signal, which in turn maximizes the SNR and + * minimizes ISI. + * + * This approach works by setting up two filterbanks; one filterbanke contains the + * signal's pulse shaping matched filter (such as a root raised cosine filter), + * where each branch of the filterbank contains a different phase of the filter. + * The second filterbank contains the derivatives of the filters in the first + * filterbank. Thinking of this in the time domain, the first filterbank contains + * filters that have a sinc shape to them. We want to align the output signal to + * be sampled at exactly the peak of the sinc shape. The derivative of the sinc + * contains a zero at the maximum point of the sinc (sinc(0) = 1, sinc(0)' = 0). + * Furthermore, the region around the zero point is relatively linear. We make + * use of this fact to generate the error signal. + * + * If the signal out of the derivative filters is d_i[n] for the ith filter, and + * the output of the matched filter is x_i[n], we calculate the error as: + * e[n] = (Re{x_i[n]} * Re{d_i[n]} + Im{x_i[n]} * Im{d_i[n]}) / 2.0 + * This equation averages the error in the real and imaginary parts. There are two + * reasons we multiply by the signal itself. First, if the symbol could be positive + * or negative going, but we want the error term to always tell us to go in the + * same direction depending on which side of the zero point we are on. The sign of + * x_i[n] adjusts the error term to do this. Second, the magnitude of x_i[n] scales + * the error term depending on the symbol's amplitude, so larger signals give us + * a stronger error term because we have more confidence in that symbol's value. + * Using the magnitude of x_i[n] instead of just the sign is especially good for + * signals with low SNR. + * + * The error signal, e[n], gives us a value proportional to how far away from the zero + * point we are in the derivative signal. We want to drive this value to zero, so we + * set up a second order loop. We have two variables for this loop; d_k is the filter + * number in the filterbank we are on and d_rate is the rate which we travel through + * the filters in the steady state. That is, due to the natural clock differences between + * the transmitter and receiver, d_rate represents that difference and would traverse + * the filter phase paths to keep the receiver locked. Thinking of this as a second-order + * PLL, the d_rate is the frequency and d_k is the phase. So we update d_rate and d_k + * using the standard loop equations based on two error signals, d_alpha and d_beta. + * We have these two values set based on each other for a critically damped system, so in + * the block constructor, we just ask for "gain," which is d_alpha while d_beta is + * equal to (gain^2)/4. + * + * The clock sync block needs to know the number of samples per second (sps), because it + * only returns a single point representing the sample. The sps can be any positive real + * number and does not need to be an integer. The filter taps must also be specified. The + * taps are generated by first conceiving of the prototype filter that would be the signal's + * matched filter. Then interpolate this by the number of filters in the filterbank. These + * are then distributed among all of the filters. So if the prototype filter was to have + * 45 taps in it, then each path of the filterbank will also have 45 taps. This is easily + * done by building the filter with the sample rate multiplied by the number of filters + * to use. + * + * The number of filters can also be set and defaults to 32. With 32 filters, you get a + * good enough resolution in the phase to produce very small, almost unnoticeable, ISI. + * Going to 64 filters can reduce this more, but after that there is very little gained + * for the extra complexity. + * + * The initial phase is another settable parameter and refers to the filter path the + * algorithm initially looks at (i.e., d_k starts at init_phase). This value defaults + * to zero, but it might be useful to start at a different phase offset, such as the mid- + * point of the filters. + * + * The final parameter is the max_rate_devitation, which defaults to 1.5. This is how far + * we allow d_rate to swing, positive or negative, from 0. Constraining the rate can help + * keep the algorithm from walking too far away to lock during times when there is no signal. + * + */ + +class gr_pfb_clock_sync_fff : public gr_block +{ + private: + /*! + * Build the polyphase filterbank timing synchronizer. + * \param sps (double) The number of samples per second in the incoming signal + * \param gain (float) The alpha gain of the control loop; beta = (gain^2)/4 by default. + * \param taps (vector<int>) The filter taps. + * \param filter_size (uint) The number of filters in the filterbank (default = 32). + * \param init_phase (float) The initial phase to look at, or which filter to start + * with (default = 0). + * \param max_rate_deviation (float) Distance from 0 d_rate can get (default = 1.5). + * + */ + friend gr_pfb_clock_sync_fff_sptr gr_make_pfb_clock_sync_fff (double sps, float gain, + const std::vector<float> &taps, + unsigned int filter_size, + float init_phase, + float max_rate_deviation); + + bool d_updated; + double d_sps; + double d_sample_num; + float d_alpha; + float d_beta; + int d_nfilters; + std::vector<gr_fir_fff*> d_filters; + std::vector<gr_fir_fff*> d_diff_filters; + std::vector< std::vector<float> > d_taps; + std::vector< std::vector<float> > d_dtaps; + float d_k; + float d_rate; + float d_rate_i; + float d_rate_f; + float d_max_dev; + int d_filtnum; + int d_taps_per_filter; + + /*! + * Build the polyphase filterbank timing synchronizer. + */ + gr_pfb_clock_sync_fff (double sps, float gain, + const std::vector<float> &taps, + unsigned int filter_size, + float init_phase, + float max_rate_deviation); + + void create_diff_taps(const std::vector<float> &newtaps, + std::vector<float> &difftaps); + +public: + ~gr_pfb_clock_sync_fff (); + + /*! + * Resets the filterbank's filter taps with the new prototype filter + */ + void set_taps (const std::vector<float> &taps, + std::vector< std::vector<float> > &ourtaps, + std::vector<gr_fir_fff*> &ourfilter); + + /*! + * Returns the taps of the matched filter + */ + std::vector<float> channel_taps(int channel); + + /*! + * Returns the taps in the derivative filter + */ + std::vector<float> diff_channel_taps(int channel); + + /*! + * Print all of the filterbank taps to screen. + */ + void print_taps(); + + /*! + * Print all of the filterbank taps of the derivative filter to screen. + */ + void print_diff_taps(); + + /*! + * Set the gain value alpha for the control loop + */ + void set_alpha(float alpha) + { + d_alpha = alpha; + } + + /*! + * Set the gain value beta for the control loop + */ + void set_beta(float beta) + { + d_beta = beta; + } + + /*! + * Set the maximum deviation from 0 d_rate can have + */ + void set_max_rate_deviation(float m) + { + d_max_dev = m; + } + + bool check_topology(int ninputs, int noutputs); + + 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_clock_sync_fff.i b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.i new file mode 100644 index 000000000..d6bb7873c --- /dev/null +++ b/gnuradio-core/src/lib/filter/gr_pfb_clock_sync_fff.i @@ -0,0 +1,54 @@ +/* -*- 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_clock_sync_fff); + +gr_pfb_clock_sync_fff_sptr gr_make_pfb_clock_sync_fff (double sps, float gain, + const std::vector<float> &taps, + unsigned int filter_size=32, + float init_phase=0, + float max_rate_deviation=1.5); + +class gr_pfb_clock_sync_fff : public gr_block +{ + private: + gr_pfb_clock_sync_fff (double sps, float gain, + const std::vector<float> &taps, + unsigned int filter_size, + float init_phase, + float max_rate_deviation); + + public: + ~gr_pfb_clock_sync_fff (); + + void set_taps (const std::vector<float> &taps, + std::vector< std::vector<float> > &ourtaps, + std::vector<gr_fir_fff*> &ourfilter); + + std::vector<float> channel_taps(int channel); + std::vector<float> diff_channel_taps(int channel); + void print_taps(); + void print_diff_taps(); + void set_alpha(float alpha); + void set_beta(float beta); + void set_max_rate_deviation(float m); +}; diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_generic.cc b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_generic.cc new file mode 100644 index 000000000..1bf4a6f4b --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_generic.cc @@ -0,0 +1,167 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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 <gri_fft_filter_ccc_generic.h> +#include <gri_fft.h> +#include <assert.h> +#include <stdexcept> +#include <cstdio> +#include <cstring> +#include <fftw3.h> + +gri_fft_filter_ccc_generic::gri_fft_filter_ccc_generic (int decimation, + const std::vector<gr_complex> &taps) + : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), d_invfft(0) +{ + set_taps(taps); +} + +gri_fft_filter_ccc_generic::~gri_fft_filter_ccc_generic () +{ + delete d_fwdfft; + delete d_invfft; +} + +#if 0 +static void +print_vector_complex(const std::string label, const std::vector<gr_complex> &x) +{ + std::cout << label; + for (unsigned i = 0; i < x.size(); i++) + std::cout << x[i] << " "; + std::cout << "\n"; +} +#endif + + +/* + * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps + */ +int +gri_fft_filter_ccc_generic::set_taps (const std::vector<gr_complex> &taps) +{ + int i = 0; + compute_sizes(taps.size()); + + d_tail.resize(tailsize()); + for (i = 0; i < tailsize(); i++) + d_tail[i] = 0; + + gr_complex *in = d_fwdfft->get_inbuf(); + gr_complex *out = d_fwdfft->get_outbuf(); + + float scale = 1.0 / d_fftsize; + + // Compute forward xform of taps. + // Copy taps into first ntaps slots, then pad with zeros + for (i = 0; i < d_ntaps; i++) + in[i] = taps[i] * scale; + + for (; i < d_fftsize; i++) + in[i] = 0; + + d_fwdfft->execute(); // do the xform + + // now copy output to d_xformed_taps + for (i = 0; i < d_fftsize; i++) + d_xformed_taps[i] = out[i]; + + return d_nsamples; +} + +// determine and set d_ntaps, d_nsamples, d_fftsize + +void +gri_fft_filter_ccc_generic::compute_sizes(int ntaps) +{ + int old_fftsize = d_fftsize; + d_ntaps = ntaps; + d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2)))); + d_nsamples = d_fftsize - d_ntaps + 1; + + if (0) + fprintf(stderr, "gri_fft_filter_ccc_generic: ntaps = %d, fftsize = %d, nsamples = %d\n", + d_ntaps, d_fftsize, d_nsamples); + + assert(d_fftsize == d_ntaps + d_nsamples -1 ); + + if (d_fftsize != old_fftsize){ // compute new plans + delete d_fwdfft; + delete d_invfft; + d_fwdfft = new gri_fft_complex(d_fftsize, true); + d_invfft = new gri_fft_complex(d_fftsize, false); + d_xformed_taps.resize(d_fftsize); + } +} + +int +gri_fft_filter_ccc_generic::filter (int nitems, const gr_complex *input, gr_complex *output) +{ + int dec_ctr = 0; + int j = 0; + int ninput_items = nitems * d_decimation; + + for (int i = 0; i < ninput_items; i += d_nsamples){ + + memcpy(d_fwdfft->get_inbuf(), &input[i], d_nsamples * sizeof(gr_complex)); + + for (j = d_nsamples; j < d_fftsize; j++) + d_fwdfft->get_inbuf()[j] = 0; + + d_fwdfft->execute(); // compute fwd xform + + gr_complex *a = d_fwdfft->get_outbuf(); + gr_complex *b = &d_xformed_taps[0]; + gr_complex *c = d_invfft->get_inbuf(); + + for (j = 0; j < d_fftsize; j+=1) { // filter in the freq domain + c[j] = a[j] * b[j]; + } + + d_invfft->execute(); // compute inv xform + + // add in the overlapping tail + + for (j = 0; j < tailsize(); j++) + d_invfft->get_outbuf()[j] += d_tail[j]; + + // copy nsamples to output + j = dec_ctr; + while (j < d_nsamples) { + *output++ = d_invfft->get_outbuf()[j]; + j += d_decimation; + } + dec_ctr = (j - d_nsamples); + + // stash the tail + memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples, + tailsize() * sizeof(gr_complex)); + } + + assert(dec_ctr == 0); + + return nitems; +} diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_generic.h b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_generic.h new file mode 100644 index 000000000..3cd9105c7 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_generic.h @@ -0,0 +1,82 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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_GRI_FFT_FILTER_CCC_GENERIC_H +#define INCLUDED_GRI_FFT_FILTER_CCC_GENERIC_H + +#include <gr_complex.h> +#include <vector> + +class gri_fft_complex; + +/*! + * \brief Fast FFT filter with gr_complex input, gr_complex output and gr_complex taps + * \ingroup filter_blk + */ +class gri_fft_filter_ccc_generic +{ + private: + int d_ntaps; + int d_nsamples; + int d_fftsize; // fftsize = ntaps + nsamples - 1 + int d_decimation; + gri_fft_complex *d_fwdfft; // forward "plan" + gri_fft_complex *d_invfft; // inverse "plan" + std::vector<gr_complex> d_tail; // state carried between blocks for overlap-add + std::vector<gr_complex> d_xformed_taps; // Fourier xformed taps + std::vector<gr_complex> d_new_taps; + + void compute_sizes(int ntaps); + int tailsize() const { return d_ntaps - 1; } + + public: + /*! + * \brief Construct an FFT filter for complex vectors with the given taps and decimation rate. + * + * This is the basic implementation for performing FFT filter for fast convolution + * in other blocks for complex vectors (such as gr_fft_filter_ccc). + * \param decimation The decimation rate of the filter (int) + * \param taps The filter taps (complex) + */ + gri_fft_filter_ccc_generic (int decimation, const std::vector<gr_complex> &taps); + ~gri_fft_filter_ccc_generic (); + + /*! + * \brief Set new taps for the filter. + * + * Sets new taps and resets the class properties to handle different sizes + * \param taps The filter taps (complex) + */ + int set_taps (const std::vector<gr_complex> &taps); + + /*! + * \brief Perform the filter operation + * + * \param nitems The number of items to produce + * \param input The input vector to be filtered + * \param output The result of the filter operation + */ + int filter (int nitems, const gr_complex *input, gr_complex *output); + +}; + +#endif /* INCLUDED_GRI_FFT_FILTER_CCC_GENERIC_H */ diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_sse.cc b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_sse.cc new file mode 100644 index 000000000..b7d925ff3 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_sse.cc @@ -0,0 +1,186 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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 <gri_fft_filter_ccc_sse.h> +#include <gri_fft.h> +#include <assert.h> +#include <stdexcept> +#include <cstdio> +#include <xmmintrin.h> +#include <fftw3.h> + +gri_fft_filter_ccc_sse::gri_fft_filter_ccc_sse (int decimation, + const std::vector<gr_complex> &taps) + : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), d_invfft(0) +{ + d_xformed_taps = (gr_complex*)fftwf_malloc(1*sizeof(gr_complex)); + set_taps(taps); +} + +gri_fft_filter_ccc_sse::~gri_fft_filter_ccc_sse () +{ + fftwf_free(d_xformed_taps); + delete d_fwdfft; + delete d_invfft; +} + +#if 0 +static void +print_vector_complex(const std::string label, const std::vector<gr_complex> &x) +{ + std::cout << label; + for (unsigned i = 0; i < x.size(); i++) + std::cout << x[i] << " "; + std::cout << "\n"; +} +#endif + + +/* + * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps + */ +int +gri_fft_filter_ccc_sse::set_taps (const std::vector<gr_complex> &taps) +{ + int i = 0; + compute_sizes(taps.size()); + + d_tail.resize(tailsize()); + for (i = 0; i < tailsize(); i++) + d_tail[i] = 0; + + gr_complex *in = d_fwdfft->get_inbuf(); + gr_complex *out = d_fwdfft->get_outbuf(); + + float scale = 1.0 / d_fftsize; + + // Compute forward xform of taps. + // Copy taps into first ntaps slots, then pad with zeros + for (i = 0; i < d_ntaps; i++) + in[i] = taps[i] * scale; + + for (; i < d_fftsize; i++) + in[i] = 0; + + d_fwdfft->execute(); // do the xform + + // now copy output to d_xformed_taps + for (i = 0; i < d_fftsize; i++) + d_xformed_taps[i] = out[i]; + + return d_nsamples; +} + +// determine and set d_ntaps, d_nsamples, d_fftsize + +void +gri_fft_filter_ccc_sse::compute_sizes(int ntaps) +{ + int old_fftsize = d_fftsize; + d_ntaps = ntaps; + d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2)))); + d_nsamples = d_fftsize - d_ntaps + 1; + + if (0) + fprintf(stderr, "gri_fft_filter_ccc_sse: ntaps = %d, fftsize = %d, nsamples = %d\n", + d_ntaps, d_fftsize, d_nsamples); + + assert(d_fftsize == d_ntaps + d_nsamples -1 ); + + if (d_fftsize != old_fftsize){ // compute new plans + delete d_fwdfft; + delete d_invfft; + d_fwdfft = new gri_fft_complex(d_fftsize, true); + d_invfft = new gri_fft_complex(d_fftsize, false); + + fftwf_free(d_xformed_taps); + d_xformed_taps = (gr_complex*)fftwf_malloc((d_fftsize)*sizeof(gr_complex)); + } +} + +int +gri_fft_filter_ccc_sse::filter (int nitems, const gr_complex *input, gr_complex *output) +{ + int dec_ctr = 0; + int j = 0; + int ninput_items = nitems * d_decimation; + + for (int i = 0; i < ninput_items; i += d_nsamples){ + + memcpy(d_fwdfft->get_inbuf(), &input[i], d_nsamples * sizeof(gr_complex)); + + for (j = d_nsamples; j < d_fftsize; j++) + d_fwdfft->get_inbuf()[j] = 0; + + d_fwdfft->execute(); // compute fwd xform + + float *a = (float*)(d_fwdfft->get_outbuf()); + float *b = (float*)(&d_xformed_taps[0]); + float *c = (float*)(d_invfft->get_inbuf()); + + __m128 x0, x1, x2, t0, t1, m; + m = _mm_set_ps(-1, 1, -1, 1); + for (j = 0; j < 2*d_fftsize; j+=4) { // filter in the freq domain + x0 = _mm_load_ps(&a[j]); + t0 = _mm_load_ps(&b[j]); + + t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(3, 3, 1, 1)); + t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 2, 0, 0)); + t1 = _mm_mul_ps(t1, m); + + x1 = _mm_mul_ps(x0, t0); + x2 = _mm_mul_ps(x0, t1); + + x2 = _mm_shuffle_ps(x2, x2, _MM_SHUFFLE(2, 3, 0, 1)); + x2 = _mm_add_ps(x1, x2); + + _mm_store_ps(&c[j], x2); + } + + d_invfft->execute(); // compute inv xform + + // add in the overlapping tail + + for (j = 0; j < tailsize(); j++) + d_invfft->get_outbuf()[j] += d_tail[j]; + + // copy nsamples to output + j = dec_ctr; + while (j < d_nsamples) { + *output++ = d_invfft->get_outbuf()[j]; + j += d_decimation; + } + dec_ctr = (j - d_nsamples); + + // stash the tail + memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples, + tailsize() * sizeof(gr_complex)); + } + + assert(dec_ctr == 0); + + return nitems; +} diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_sse.h b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_sse.h new file mode 100644 index 000000000..d1c54f01f --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_ccc_sse.h @@ -0,0 +1,82 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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_GRI_FFT_FILTER_CCC_SSE_H +#define INCLUDED_GRI_FFT_FILTER_CCC_SSE_H + +#include <gr_complex.h> +#include <vector> + +class gri_fft_complex; + +/*! + * \brief Fast FFT filter with gr_complex input, gr_complex output and gr_complex taps + * \ingroup filter_blk + */ +class gri_fft_filter_ccc_sse +{ + private: + int d_ntaps; + int d_nsamples; + int d_fftsize; // fftsize = ntaps + nsamples - 1 + int d_decimation; + gri_fft_complex *d_fwdfft; // forward "plan" + gri_fft_complex *d_invfft; // inverse "plan" + std::vector<gr_complex> d_tail; // state carried between blocks for overlap-add + gr_complex *d_xformed_taps; + std::vector<gr_complex> d_new_taps; + + void compute_sizes(int ntaps); + int tailsize() const { return d_ntaps - 1; } + + public: + /*! + * \brief Construct an FFT filter for complex vectors with the given taps and decimation rate. + * + * This is the basic implementation for performing FFT filter for fast convolution + * in other blocks for complex vectors (such as gr_fft_filter_ccc). + * \param decimation The decimation rate of the filter (int) + * \param taps The filter taps (complex) + */ + gri_fft_filter_ccc_sse (int decimation, const std::vector<gr_complex> &taps); + ~gri_fft_filter_ccc_sse (); + + /*! + * \brief Set new taps for the filter. + * + * Sets new taps and resets the class properties to handle different sizes + * \param taps The filter taps (complex) + */ + int set_taps (const std::vector<gr_complex> &taps); + + /*! + * \brief Perform the filter operation + * + * \param nitems The number of items to produce + * \param input The input vector to be filtered + * \param output The result of the filter operation + */ + int filter (int nitems, const gr_complex *input, gr_complex *output); + +}; + +#endif /* INCLUDED_GRI_FFT_FILTER_CCC_SSE_H */ diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_fff_generic.cc b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_generic.cc new file mode 100644 index 000000000..74058fa93 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_generic.cc @@ -0,0 +1,158 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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 <gri_fft_filter_fff_generic.h> +#include <gri_fft.h> +#include <assert.h> +#include <stdexcept> +#include <cstdio> +#include <cstring> + +gri_fft_filter_fff_generic::gri_fft_filter_fff_generic (int decimation, + const std::vector<float> &taps) + : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), d_invfft(0) +{ + set_taps(taps); +} + +gri_fft_filter_fff_generic::~gri_fft_filter_fff_generic () +{ + delete d_fwdfft; + delete d_invfft; +} + +/* + * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps + */ +int +gri_fft_filter_fff_generic::set_taps (const std::vector<float> &taps) +{ + int i = 0; + compute_sizes(taps.size()); + + d_tail.resize(tailsize()); + for (i = 0; i < tailsize(); i++) + d_tail[i] = 0; + + float *in = d_fwdfft->get_inbuf(); + gr_complex *out = d_fwdfft->get_outbuf(); + + float scale = 1.0 / d_fftsize; + + // Compute forward xform of taps. + // Copy taps into first ntaps slots, then pad with zeros + for (i = 0; i < d_ntaps; i++) + in[i] = taps[i] * scale; + + for (; i < d_fftsize; i++) + in[i] = 0; + + d_fwdfft->execute(); // do the xform + + // now copy output to d_xformed_taps + for (i = 0; i < d_fftsize/2+1; i++) + d_xformed_taps[i] = out[i]; + + return d_nsamples; +} + +// determine and set d_ntaps, d_nsamples, d_fftsize + +void +gri_fft_filter_fff_generic::compute_sizes(int ntaps) +{ + int old_fftsize = d_fftsize; + d_ntaps = ntaps; + d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2)))); + d_nsamples = d_fftsize - d_ntaps + 1; + + if (0) + fprintf(stderr, "gri_fft_filter_fff_generic: ntaps = %d, fftsize = %d, nsamples = %d\n", + d_ntaps, d_fftsize, d_nsamples); + + assert(d_fftsize == d_ntaps + d_nsamples -1 ); + + if (d_fftsize != old_fftsize){ // compute new plans + delete d_fwdfft; + delete d_invfft; + d_fwdfft = new gri_fft_real_fwd(d_fftsize); + d_invfft = new gri_fft_real_rev(d_fftsize); + d_xformed_taps.resize(d_fftsize/2+1); + } +} + +int +gri_fft_filter_fff_generic::filter (int nitems, const float *input, float *output) +{ + int dec_ctr = 0; + int j = 0; + int ninput_items = nitems * d_decimation; + + for (int i = 0; i < ninput_items; i += d_nsamples){ + + memcpy(d_fwdfft->get_inbuf(), &input[i], d_nsamples * sizeof(float)); + + for (j = d_nsamples; j < d_fftsize; j++) + d_fwdfft->get_inbuf()[j] = 0; + + d_fwdfft->execute(); // compute fwd xform + + gr_complex *a = d_fwdfft->get_outbuf(); + gr_complex *b = &d_xformed_taps[0]; + gr_complex *c = d_invfft->get_inbuf(); + + for (j = 0; j < d_fftsize/2+1; j++) { // filter in the freq domain + c[j] = a[j] * b[j]; + } + + d_invfft->execute(); // compute inv xform + + // add in the overlapping tail + + for (j = 0; j < tailsize(); j++) + d_invfft->get_outbuf()[j] += d_tail[j]; + + // copy nsamples to output + + //memcpy(out, d_invfft->get_outbuf(), d_nsamples * sizeof(float)); + //out += d_nsamples; + + j = dec_ctr; + while (j < d_nsamples) { + *output++ = d_invfft->get_outbuf()[j]; + j += d_decimation; + } + dec_ctr = (j - d_nsamples); + + // stash the tail + memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples, + tailsize() * sizeof(float)); + } + + assert(dec_ctr == 0); + + return nitems; +} diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_fff_generic.h b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_generic.h new file mode 100644 index 000000000..6c31632d5 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_generic.h @@ -0,0 +1,80 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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_GRI_FFT_FILTER_FFF_GENERIC_H +#define INCLUDED_GRI_FFT_FILTER_FFF_GENERIC_H + +#include <gr_complex.h> +#include <vector> + +class gri_fft_real_fwd; +class gri_fft_real_rev; + +class gri_fft_filter_fff_generic +{ + private: + int d_ntaps; + int d_nsamples; + int d_fftsize; // fftsize = ntaps + nsamples - 1 + int d_decimation; + gri_fft_real_fwd *d_fwdfft; // forward "plan" + gri_fft_real_rev *d_invfft; // inverse "plan" + std::vector<float> d_tail; // state carried between blocks for overlap-add + std::vector<gr_complex> d_xformed_taps; // Fourier xformed taps + std::vector<float> d_new_taps; + + + void compute_sizes(int ntaps); + int tailsize() const { return d_ntaps - 1; } + + public: + /*! + * \brief Construct a FFT filter for float vectors with the given taps and decimation rate. + * + * This is the basic implementation for performing FFT filter for fast convolution + * in other blocks for floating point vectors (such as gr_fft_filter_fff). + * \param decimation The decimation rate of the filter (int) + * \param taps The filter taps (float) + */ + gri_fft_filter_fff_generic (int decimation, const std::vector<float> &taps); + ~gri_fft_filter_fff_generic (); + + /*! + * \brief Set new taps for the filter. + * + * Sets new taps and resets the class properties to handle different sizes + * \param taps The filter taps (float) + */ + int set_taps (const std::vector<float> &taps); + + /*! + * \brief Perform the filter operation + * + * \param nitems The number of items to produce + * \param input The input vector to be filtered + * \param output The result of the filter operation + */ + int filter (int nitems, const float *input, float *output); + +}; + +#endif /* INCLUDED_GRI_FFT_FILTER_FFF_GENERIC_H */ diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_fff_sse.cc b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_sse.cc new file mode 100644 index 000000000..2680e6594 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_sse.cc @@ -0,0 +1,184 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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 <gri_fft_filter_fff_sse.h> +#include <gri_fft.h> +#include <assert.h> +#include <stdexcept> +#include <cstdio> +#include <xmmintrin.h> +#include <fftw3.h> + +gri_fft_filter_fff_sse::gri_fft_filter_fff_sse (int decimation, + const std::vector<float> &taps) + : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), d_invfft(0) +{ + d_xformed_taps = (gr_complex*)fftwf_malloc(1*sizeof(gr_complex)); + set_taps(taps); +} + +gri_fft_filter_fff_sse::~gri_fft_filter_fff_sse () +{ + fftwf_free(d_xformed_taps); + delete d_fwdfft; + delete d_invfft; +} + +/* + * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps + */ +int +gri_fft_filter_fff_sse::set_taps (const std::vector<float> &taps) +{ + int i = 0; + compute_sizes(taps.size()); + + d_tail.resize(tailsize()); + for (i = 0; i < tailsize(); i++) + d_tail[i] = 0; + + float *in = d_fwdfft->get_inbuf(); + gr_complex *out = d_fwdfft->get_outbuf(); + + float scale = 1.0 / d_fftsize; + + // Compute forward xform of taps. + // Copy taps into first ntaps slots, then pad with zeros + for (i = 0; i < d_ntaps; i++) + in[i] = taps[i] * scale; + + for (; i < d_fftsize; i++) + in[i] = 0; + + d_fwdfft->execute(); // do the xform + + // now copy output to d_xformed_taps + for (i = 0; i < d_fftsize/2+1; i++) + d_xformed_taps[i] = out[i]; + + return d_nsamples; +} + +// determine and set d_ntaps, d_nsamples, d_fftsize + +void +gri_fft_filter_fff_sse::compute_sizes(int ntaps) +{ + int old_fftsize = d_fftsize; + d_ntaps = ntaps; + d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2)))); + d_nsamples = d_fftsize - d_ntaps + 1; + + if (0) + fprintf(stderr, "gri_fft_filter_fff_sse: ntaps = %d, fftsize = %d, nsamples = %d\n", + d_ntaps, d_fftsize, d_nsamples); + + assert(d_fftsize == d_ntaps + d_nsamples -1 ); + + if (d_fftsize != old_fftsize){ // compute new plans + delete d_fwdfft; + delete d_invfft; + d_fwdfft = new gri_fft_real_fwd(d_fftsize); + d_invfft = new gri_fft_real_rev(d_fftsize); + //d_xformed_taps.resize(d_fftsize/2+1); + + fftwf_free(d_xformed_taps); + d_xformed_taps = (gr_complex*)fftwf_malloc((d_fftsize/2+1)*sizeof(gr_complex)); + } +} + +int +gri_fft_filter_fff_sse::filter (int nitems, const float *input, float *output) +{ + int dec_ctr = 0; + int j = 0; + int ninput_items = nitems * d_decimation; + + for (int i = 0; i < ninput_items; i += d_nsamples){ + + memcpy(d_fwdfft->get_inbuf(), &input[i], d_nsamples * sizeof(float)); + + for (j = d_nsamples; j < d_fftsize; j++) + d_fwdfft->get_inbuf()[j] = 0; + + d_fwdfft->execute(); // compute fwd xform + + float *a = (float*)(d_fwdfft->get_outbuf()); + float *b = (float*)(&d_xformed_taps[0]); + float *c = (float*)(d_invfft->get_inbuf()); + + __m128 x0, x1, x2, t0, t1, m; + m = _mm_set_ps(-1, 1, -1, 1); + for (j = 0; j < d_fftsize; j+=4) { // filter in the freq domain + x0 = _mm_load_ps(&a[j]); + t0 = _mm_load_ps(&b[j]); + + t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(3, 3, 1, 1)); + t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 2, 0, 0)); + t1 = _mm_mul_ps(t1, m); + + x1 = _mm_mul_ps(x0, t0); + x2 = _mm_mul_ps(x0, t1); + + x2 = _mm_shuffle_ps(x2, x2, _MM_SHUFFLE(2, 3, 0, 1)); + x2 = _mm_add_ps(x1, x2); + + _mm_store_ps(&c[j], x2); + } + + // Finish off the last one; do the complex multiply as floats + j = d_fftsize/2; + c[j] = (a[j] * b[j]) - (a[j+1] * b[j+1]); + c[j+1] = (a[j] * b[j+1]) + (a[j+1] * b[j]); + + d_invfft->execute(); // compute inv xform + + // add in the overlapping tail + + for (j = 0; j < tailsize(); j++) + d_invfft->get_outbuf()[j] += d_tail[j]; + + // copy nsamples to output + + //memcpy(out, d_invfft->get_outbuf(), d_nsamples * sizeof(float)); + //out += d_nsamples; + + j = dec_ctr; + while (j < d_nsamples) { + *output++ = d_invfft->get_outbuf()[j]; + j += d_decimation; + } + dec_ctr = (j - d_nsamples); + + // stash the tail + memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples, + tailsize() * sizeof(float)); + } + + assert(dec_ctr == 0); + + return nitems; +} diff --git a/gnuradio-core/src/lib/filter/gri_fft_filter_fff_sse.h b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_sse.h new file mode 100644 index 000000000..8258bb824 --- /dev/null +++ b/gnuradio-core/src/lib/filter/gri_fft_filter_fff_sse.h @@ -0,0 +1,81 @@ +/* -*- c++ -*- */ +/* + * Copyright 2010 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_GRI_FFT_FILTER_FFF_SSE_H +#define INCLUDED_GRI_FFT_FILTER_FFF_SSE_H + +#include <gr_complex.h> +#include <vector> + +class gri_fft_real_fwd; +class gri_fft_real_rev; + +class gri_fft_filter_fff_sse +{ + private: + int d_ntaps; + int d_nsamples; + int d_fftsize; // fftsize = ntaps + nsamples - 1 + int d_decimation; + gri_fft_real_fwd *d_fwdfft; // forward "plan" + gri_fft_real_rev *d_invfft; // inverse "plan" + std::vector<float> d_tail; // state carried between blocks for overlap-add + //std::vector<gr_complex> d_xformed_taps; // Fourier xformed taps + gr_complex *d_xformed_taps; + std::vector<float> d_new_taps; + + + void compute_sizes(int ntaps); + int tailsize() const { return d_ntaps - 1; } + + public: + /*! + * \brief Construct a FFT filter for float vectors with the given taps and decimation rate. + * + * This is the basic implementation for performing FFT filter for fast convolution + * in other blocks for floating point vectors (such as gr_fft_filter_fff). + * \param decimation The decimation rate of the filter (int) + * \param taps The filter taps (float) + */ + gri_fft_filter_fff_sse (int decimation, const std::vector<float> &taps); + ~gri_fft_filter_fff_sse (); + + /*! + * \brief Set new taps for the filter. + * + * Sets new taps and resets the class properties to handle different sizes + * \param taps The filter taps (float) + */ + int set_taps (const std::vector<float> &taps); + + /*! + * \brief Perform the filter operation + * + * \param nitems The number of items to produce + * \param input The input vector to be filtered + * \param output The result of the filter operation + */ + int filter (int nitems, const float *input, float *output); + +}; + +#endif /* INCLUDED_GRI_FFT_FILTER_FFF_SSE_H */ diff --git a/gnuradio-core/src/lib/general/Makefile.am b/gnuradio-core/src/lib/general/Makefile.am index cf6ff1e65..b5f5c346b 100644 --- a/gnuradio-core/src/lib/general/Makefile.am +++ b/gnuradio-core/src/lib/general/Makefile.am @@ -34,6 +34,7 @@ EXTRA_DIST = \ gr_constants.cc.in libgeneral_la_SOURCES = \ + gr_additive_scrambler_bb.cc \ gr_agc_cc.cc \ gr_agc_ff.cc \ gr_agc2_cc.cc \ @@ -76,6 +77,7 @@ libgeneral_la_SOURCES = \ gr_fft_vcc_fftw.cc \ gr_fft_vfc.cc \ gr_firdes.cc \ + gr_fll_band_edge_cc.cc \ gr_float_to_char.cc \ gr_float_to_complex.cc \ gr_float_to_short.cc \ @@ -187,6 +189,7 @@ libgeneral_qa_la_SOURCES = \ qa_gri_lfsr.cc grinclude_HEADERS = \ + gr_additive_scrambler_bb.h \ gr_agc_cc.h \ gr_agc_ff.h \ gr_agc2_cc.h \ @@ -229,6 +232,7 @@ grinclude_HEADERS = \ gr_fft_vcc_fftw.h \ gr_fft_vfc.h \ gr_firdes.h \ + gr_fll_band_edge_cc.h \ gr_float_to_char.h \ gr_float_to_complex.h \ gr_float_to_short.h \ @@ -358,6 +362,7 @@ noinst_HEADERS = \ if PYTHON swiginclude_HEADERS = \ general.i \ + gr_additive_scrambler_bb.i \ gr_agc_cc.i \ gr_agc_ff.i \ gr_agc2_cc.i \ @@ -396,6 +401,7 @@ swiginclude_HEADERS = \ gr_fft_vcc.i \ gr_fft_vfc.i \ gr_firdes.i \ + gr_fll_band_edge_cc.i \ gr_float_to_char.i \ gr_float_to_complex.i \ gr_float_to_short.i \ diff --git a/gnuradio-core/src/lib/general/general.i b/gnuradio-core/src/lib/general/general.i index e1161c8eb..6929f1e6e 100644 --- a/gnuradio-core/src/lib/general/general.i +++ b/gnuradio-core/src/lib/general/general.i @@ -140,6 +140,8 @@ #include <gr_wavelet_ff.h> #include <gr_wvps_ff.h> #include <gr_copy.h> +#include <gr_fll_band_edge_cc.h> +#include <gr_additive_scrambler_bb.h> %} %include "gr_nop.i" @@ -260,3 +262,5 @@ %include "gr_wavelet_ff.i" %include "gr_wvps_ff.i" %include "gr_copy.i" +%include "gr_fll_band_edge_cc.i" +%include "gr_additive_scrambler_bb.i" diff --git a/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.cc b/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.cc new file mode 100644 index 000000000..91e02c2d3 --- /dev/null +++ b/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.cc @@ -0,0 +1,65 @@ +/* -*- c++ -*- */ +/* + * Copyright 2008,2010 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 <gr_additive_scrambler_bb.h> +#include <gr_io_signature.h> + +gr_additive_scrambler_bb_sptr +gr_make_additive_scrambler_bb(int mask, int seed, int len, int count) +{ + return gr_additive_scrambler_bb_sptr(new gr_additive_scrambler_bb(mask, seed, len, count)); +} + +gr_additive_scrambler_bb::gr_additive_scrambler_bb(int mask, int seed, int len, int count) + : gr_sync_block("additive_scrambler_bb", + gr_make_io_signature (1, 1, sizeof (unsigned char)), + gr_make_io_signature (1, 1, sizeof (unsigned char))), + d_lfsr(mask, seed, len), + d_count(count), + d_bits(0) +{ +} + +int +gr_additive_scrambler_bb::work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + const unsigned char *in = (const unsigned char *) input_items[0]; + unsigned char *out = (unsigned char *) output_items[0]; + + for (int i = 0; i < noutput_items; i++) { + out[i] = in[i]^d_lfsr.next_bit(); + if (d_count > 0) { + if (++d_bits == d_count) { + d_lfsr.reset(); + d_bits = 0; + } + } + } + + return noutput_items; +} diff --git a/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.h b/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.h new file mode 100644 index 000000000..6c9493050 --- /dev/null +++ b/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.h @@ -0,0 +1,67 @@ +/* -*- c++ -*- */ +/* + * Copyright 2008,2010 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_ADDITIVE_SCRAMBLER_BB_H +#define INCLUDED_GR_ADDITIVE_SCRAMBLER_BB_H + +#include <gr_sync_block.h> +#include "gri_lfsr.h" + +class gr_additive_scrambler_bb; +typedef boost::shared_ptr<gr_additive_scrambler_bb> gr_additive_scrambler_bb_sptr; + +gr_additive_scrambler_bb_sptr gr_make_additive_scrambler_bb(int mask, int seed, int len, int count=0); + +/*! + * Scramble an input stream using an LFSR. This block works on the LSB only + * of the input data stream, i.e., on an "unpacked binary" stream, and + * produces the same format on its output. + * + * \param mask Polynomial mask for LFSR + * \param seed Initial shift register contents + * \param len Shift register length + * \param count Number of bits after which shift register is reset, 0=never + * + * The scrambler works by XORing the incoming bit stream by the output of + * the LFSR. Optionally, after 'count' bits have been processed, the shift + * register is reset to the seed value. This allows processing fixed length + * vectors of samples. + * + * \ingroup coding_blk + */ + +class gr_additive_scrambler_bb : public gr_sync_block +{ + friend gr_additive_scrambler_bb_sptr gr_make_additive_scrambler_bb(int mask, int seed, int len, int count); + + gri_lfsr d_lfsr; + int d_count; + int d_bits; + + gr_additive_scrambler_bb(int mask, int seed, int len, int count); + +public: + int work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); +}; + +#endif /* INCLUDED_GR_ADDITIVE_SCRAMBLER_BB_H */ diff --git a/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.i b/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.i new file mode 100644 index 000000000..0ca9c1cd7 --- /dev/null +++ b/gnuradio-core/src/lib/general/gr_additive_scrambler_bb.i @@ -0,0 +1,31 @@ +/* -*- c++ -*- */ +/* + * Copyright 2008,2010 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,additive_scrambler_bb); + +gr_additive_scrambler_bb_sptr gr_make_additive_scrambler_bb(int mask, int seed, int len, int count=0); + +class gr_additive_scrambler_bb : public gr_sync_block +{ +private: + gr_additive_scrambler_bb(int mask, int seed, int len, int count); +}; diff --git a/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.cc b/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.cc new file mode 100644 index 000000000..7f2c468b7 --- /dev/null +++ b/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.cc @@ -0,0 +1,213 @@ +/* -*- 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 <gr_fll_band_edge_cc.h> +#include <gr_fir_ccc.h> +#include <gr_fir_util.h> +#include <gri_fft.h> +#include <gr_io_signature.h> +#include <gr_expj.h> +#include <gr_math.h> +#include <cstdio> + +#define M_TWOPI (2*M_PI) + +float sinc(float x) +{ + if(x == 0) + return 1; + else + return sin(M_PI*x)/(M_PI*x); +} + + + +gr_fll_band_edge_cc_sptr gr_make_fll_band_edge_cc (float samps_per_sym, float rolloff, + int filter_size, float gain_alpha, float gain_beta) +{ + return gr_fll_band_edge_cc_sptr (new gr_fll_band_edge_cc (samps_per_sym, rolloff, + filter_size, gain_alpha, gain_beta)); +} + + +static int ios[] = {sizeof(gr_complex), sizeof(float), sizeof(float), sizeof(gr_complex)}; +static std::vector<int> iosig(ios, ios+sizeof(ios)/sizeof(int)); +gr_fll_band_edge_cc::gr_fll_band_edge_cc (float samps_per_sym, float rolloff, + int filter_size, float alpha, float beta) + : gr_sync_block ("fll_band_edge_cc", + gr_make_io_signature (1, 1, sizeof(gr_complex)), + gr_make_io_signaturev (1, 4, iosig)), + d_alpha(alpha), d_beta(beta), d_updated (false) +{ + // base this on the number of samples per symbol + d_max_freq = M_TWOPI * (2.0/samps_per_sym); + d_min_freq = -M_TWOPI * (2.0/samps_per_sym); + + d_freq = 0; + d_phase = 0; + + set_alpha(alpha); + + design_filter(samps_per_sym, rolloff, filter_size); +} + +gr_fll_band_edge_cc::~gr_fll_band_edge_cc () +{ + delete d_filter_lower; + delete d_filter_upper; +} + +void +gr_fll_band_edge_cc::set_alpha(float alpha) +{ + //float eta = sqrt(2.0)/2.0; + //float theta = alpha; + //d_alpha = (4*eta*theta) / (1.0 + 2.0*eta*theta + theta*theta); + //d_beta = (4*theta*theta) / (1.0 + 2.0*eta*theta + theta*theta); + d_alpha = alpha; +} + +void +gr_fll_band_edge_cc::design_filter(float samps_per_sym, float rolloff, int filter_size) +{ + int M = rint(filter_size / samps_per_sym); + float power = 0; + std::vector<float> bb_taps; + for(int i = 0; i < filter_size; i++) { + float k = -M + i*2.0/samps_per_sym; + float tap = sinc(rolloff*k - 0.5) + sinc(rolloff*k + 0.5); + power += tap; + + bb_taps.push_back(tap); + } + + int N = (bb_taps.size() - 1.0)/2.0; + std::vector<gr_complex> taps_lower; + std::vector<gr_complex> taps_upper; + for(unsigned int i = 0; i < bb_taps.size(); i++) { + float tap = bb_taps[i] / power; + + float k = (-N + (int)i)/(2.0*samps_per_sym); + + gr_complex t1 = tap * gr_expj(-2*M_PI*(1+rolloff)*k); + gr_complex t2 = tap * gr_expj(2*M_PI*(1+rolloff)*k); + + taps_lower.push_back(t1); + taps_upper.push_back(t2); + } + + std::vector<gr_complex> vtaps(0, taps_lower.size()); + d_filter_upper = gr_fir_util::create_gr_fir_ccc(vtaps); + d_filter_lower = gr_fir_util::create_gr_fir_ccc(vtaps); + + d_filter_lower->set_taps(taps_lower); + d_filter_upper->set_taps(taps_upper); + + d_updated = true; + + // Set the history to ensure enough input items for each filter + set_history(filter_size+1); + +} + +void +gr_fll_band_edge_cc::print_taps() +{ + unsigned int i; + std::vector<gr_complex> taps_upper = d_filter_upper->get_taps(); + std::vector<gr_complex> taps_lower = d_filter_lower->get_taps(); + + printf("Upper Band-edge: ["); + for(i = 0; i < taps_upper.size(); i++) { + printf(" %.4e + %.4ej,", taps_upper[i].real(), taps_upper[i].imag()); + } + printf("]\n\n"); + + printf("Lower Band-edge: ["); + for(i = 0; i < taps_lower.size(); i++) { + printf(" %.4e + %.4ej,", taps_lower[i].real(), taps_lower[i].imag()); + } + printf("]\n\n"); +} + +int +gr_fll_band_edge_cc::work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + const gr_complex *in = (const gr_complex *) input_items[0]; + gr_complex *out = (gr_complex *) output_items[0]; + + float *frq, *phs; + gr_complex *err; + if(output_items.size() > 2) { + frq = (float *) output_items[1]; + phs = (float *) output_items[2]; + err = (gr_complex *) output_items[3]; + } + + if (d_updated) { + d_updated = false; + return 0; // history requirements may have changed. + } + + int i; + gr_complex nco_out; + gr_complex out_upper, out_lower; + float error; + float avg_k = 0.1; + for(i = 0; i < noutput_items; i++) { + nco_out = gr_expj(d_phase); + out[i] = in[i] * nco_out; + + out_upper = (d_filter_upper->filter(&out[i])); + out_lower = (d_filter_lower->filter(&out[i])); + error = -real((out_upper + out_lower) * conj(out_upper - out_lower)); + d_error = avg_k*error + avg_k*d_error; // average error + + d_freq = d_freq + d_beta * d_error; + d_phase = d_phase + d_freq + d_alpha * d_error; + + if(d_phase > M_PI) + d_phase -= M_TWOPI; + else if(d_phase < -M_PI) + d_phase += M_TWOPI; + + if (d_freq > d_max_freq) + d_freq = d_max_freq; + else if (d_freq < d_min_freq) + d_freq = d_min_freq; + + if(output_items.size() > 2) { + frq[i] = d_freq; + phs[i] = d_phase; + err[i] = d_error; + } + } + + + return noutput_items; +} diff --git a/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.h b/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.h new file mode 100644 index 000000000..db060793e --- /dev/null +++ b/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.h @@ -0,0 +1,139 @@ +/* -*- 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_FLL_BAND_EDGE_CC_H +#define INCLUDED_GR_FLL_BAND_EDGE_CC_H + +#include <gr_sync_block.h> + +class gr_fll_band_edge_cc; +typedef boost::shared_ptr<gr_fll_band_edge_cc> gr_fll_band_edge_cc_sptr; +gr_fll_band_edge_cc_sptr gr_make_fll_band_edge_cc (float samps_per_sym, float rolloff, + int filter_size, float alpha, float beta); + +class gr_fir_ccc; +class gri_fft_complex; + +/*! + * \class gr_fll_band_edge_cc + * \brief Frequency Lock Loop using band-edge filters + * + * \ingroup general + * + * The frequency lock loop derives a band-edge filter that covers the upper and lower bandwidths + * of a digitally-modulated signal. The bandwidth range is determined by the excess bandwidth + * (e.g., rolloff factor) of the modulated signal. The placement in frequency of the band-edges + * is determined by the oversampling ratio (number of samples per symbol) and the excess bandwidth. + * The size of the filters should be fairly large so as to average over a number of symbols. + * + * The FLL works by filtering the upper and lower band edges into x_u(t) and x_l(t), respectively. + * These are combined to form cc(t) = x_u(t) + x_l(t) and ss(t) = x_u(t) - x_l(t). Combining + * these to form the signal e(t) = Re{cc(t) \\times ss(t)^*} (where ^* is the complex conjugate) + * provides an error signal at the DC term that is directly proportional to the carrier frequency. + * We then make a second-order loop using the error signal that is the running average of e(t). + * + * In theory, the band-edge filter is the derivative of the matched filter in frequency, + * (H_be(f) = \\frac{H(f)}{df}. In practice, this comes down to a quarter sine wave at the point + * of the matched filter's rolloff (if it's a raised-cosine, the derivative of a cosine is a sine). + * Extend this sine by another quarter wave to make a half wave around the band-edges is equivalent + * in time to the sum of two sinc functions. The baseband filter fot the band edges is therefore + * derived from this sum of sincs. The band edge filters are then just the baseband signal + * modulated to the correct place in frequency. All of these calculations are done in the + * 'design_filter' function. + * + * Note: We use FIR filters here because the filters have to have a flat phase response over the + * entire frequency range to allow their comparisons to be valid. + */ + +class gr_fll_band_edge_cc : public gr_sync_block +{ + private: + /*! + * Build the FLL + * \param samps_per_sym (float) Number of samples per symbol of signal + * \param rolloff (float) Rolloff factor of signal + * \param filter_size (int) Size (in taps) of the filter + * \param alpha (float) Loop gain 1 + * \param beta (float) Loop gain 2 + */ + friend gr_fll_band_edge_cc_sptr gr_make_fll_band_edge_cc (float samps_per_sym, float rolloff, + int filter_size, float alpha, float beta); + + float d_alpha; + float d_beta; + float d_max_freq; + float d_min_freq; + + gr_fir_ccc* d_filter_upper; + gr_fir_ccc* d_filter_lower; + bool d_updated; + float d_error; + float d_freq; + float d_phase; + + /*! + * Build the FLL + * \param samps_per_sym (float) number of samples per symbol + * \param rolloff (float) Rolloff (excess bandwidth) of signal filter + * \param filter_size (int) number of filter taps to generate + * \param alpha (float) Alpha gain in the control loop + * \param beta (float) Beta gain in the control loop + */ + gr_fll_band_edge_cc(float samps_per_sym, float rolloff, + int filter_size, float alpha, float beta); + +public: + ~gr_fll_band_edge_cc (); + + /*! + * Design the band-edge filter based on the number of samples per symbol, + * filter rolloff factor, and the filter size + * \param samps_per_sym (float) Number of samples per symbol of signal + * \param rolloff (float) Rolloff factor of signal + * \param filter_size (int) Size (in taps) of the filter + */ + void design_filter(float samps_per_sym, float rolloff, int filter_size); + + /*! + * Set the alpha gainvalue + * \param alpha (float) new gain value + */ + void set_alpha(float alpha); + + /*! + * Set the beta gain value + * \param beta (float) new gain value + */ + void set_beta(float beta) { d_beta = beta; } + + /*! + * Print the 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/general/gr_fll_band_edge_cc.i b/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.i new file mode 100644 index 000000000..c9c792c8a --- /dev/null +++ b/gnuradio-core/src/lib/general/gr_fll_band_edge_cc.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,fll_band_edge_cc); + +gr_fll_band_edge_cc_sptr gr_make_fll_band_edge_cc (float samps_per_sym, float rolloff, + int filter_size, float alpha, float beta); + +class gr_fll_band_edge_cc : public gr_sync_block +{ + private: + gr_fll_band_edge_cc (float samps_per_sym, float rolloff, + int filter_size, float alpha, float beta); + + public: + ~gr_fll_band_edge_cc (); + + void set_alpha (float alpha); + void set_beta (float beta); + void design_filter(float samps_per_sym, float rolloff, int filter_size); + void print_taps(); +}; diff --git a/gnuradio-core/src/lib/general/gr_head.cc b/gnuradio-core/src/lib/general/gr_head.cc index 01035ffcd..b52735c06 100644 --- a/gnuradio-core/src/lib/general/gr_head.cc +++ b/gnuradio-core/src/lib/general/gr_head.cc @@ -27,7 +27,7 @@ #include <gr_io_signature.h> #include <string.h> -gr_head::gr_head (size_t sizeof_stream_item, int nitems) +gr_head::gr_head (size_t sizeof_stream_item, unsigned long long nitems) : gr_sync_block ("head", gr_make_io_signature (1, 1, sizeof_stream_item), gr_make_io_signature (1, 1, sizeof_stream_item)), @@ -36,7 +36,7 @@ gr_head::gr_head (size_t sizeof_stream_item, int nitems) } gr_head_sptr -gr_make_head (size_t sizeof_stream_item, int nitems) +gr_make_head (size_t sizeof_stream_item, unsigned long long nitems) { return gnuradio::get_initial_sptr(new gr_head (sizeof_stream_item, nitems)); } @@ -49,7 +49,7 @@ gr_head::work (int noutput_items, if (d_ncopied_items >= d_nitems) return -1; // Done! - unsigned n = std::min (d_nitems - d_ncopied_items, noutput_items); + unsigned n = std::min (d_nitems - d_ncopied_items, (unsigned long long) noutput_items); if (n == 0) return 0; diff --git a/gnuradio-core/src/lib/general/gr_head.h b/gnuradio-core/src/lib/general/gr_head.h index 430d5f8b9..f7eee1064 100644 --- a/gnuradio-core/src/lib/general/gr_head.h +++ b/gnuradio-core/src/lib/general/gr_head.h @@ -38,11 +38,11 @@ typedef boost::shared_ptr<gr_head> gr_head_sptr; class gr_head : public gr_sync_block { - friend gr_head_sptr gr_make_head (size_t sizeof_stream_item, int nitems); - gr_head (size_t sizeof_stream_item, int nitems); + friend gr_head_sptr gr_make_head (size_t sizeof_stream_item, unsigned long long nitems); + gr_head (size_t sizeof_stream_item, unsigned long long nitems); - int d_nitems; - int d_ncopied_items; + unsigned long long d_nitems; + unsigned long long d_ncopied_items; public: int work (int noutput_items, @@ -53,7 +53,7 @@ class gr_head : public gr_sync_block }; gr_head_sptr -gr_make_head (size_t sizeof_stream_item, int nitems); +gr_make_head (size_t sizeof_stream_item, unsigned long long nitems); #endif /* INCLUDED_GR_HEAD_H */ diff --git a/gnuradio-core/src/lib/general/gr_head.i b/gnuradio-core/src/lib/general/gr_head.i index 2a88b885f..3aece9601 100644 --- a/gnuradio-core/src/lib/general/gr_head.i +++ b/gnuradio-core/src/lib/general/gr_head.i @@ -22,7 +22,7 @@ GR_SWIG_BLOCK_MAGIC(gr,head); -gr_head_sptr gr_make_head(size_t sizeof_stream_item, int nitems); +gr_head_sptr gr_make_head(size_t sizeof_stream_item, unsigned long long nitems); class gr_head : public gr_block { gr_head(); diff --git a/gnuradio-core/src/lib/general/gr_ofdm_sampler.cc b/gnuradio-core/src/lib/general/gr_ofdm_sampler.cc index 74bd65a50..7f6b2b01c 100644 --- a/gnuradio-core/src/lib/general/gr_ofdm_sampler.cc +++ b/gnuradio-core/src/lib/general/gr_ofdm_sampler.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2007,2008 Free Software Foundation, Inc. + * Copyright 2007,2008,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -45,6 +45,7 @@ gr_ofdm_sampler::gr_ofdm_sampler (unsigned int fft_length, gr_make_io_signature2 (2, 2, sizeof (gr_complex)*fft_length, sizeof(char)*fft_length)), d_state(STATE_NO_SIG), d_timeout_max(timeout), d_fft_length(fft_length), d_symbol_length(symbol_length) { + set_relative_rate(1.0/(double) fft_length); // buffer allocator hint } void diff --git a/gnuradio-core/src/lib/general/gri_lfsr.h b/gnuradio-core/src/lib/general/gri_lfsr.h index 715da78a9..f691e36ec 100644 --- a/gnuradio-core/src/lib/general/gri_lfsr.h +++ b/gnuradio-core/src/lib/general/gri_lfsr.h @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2008 Free Software Foundation, Inc. + * Copyright 2008,2010 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -86,6 +86,7 @@ class gri_lfsr private: uint32_t d_shift_register; uint32_t d_mask; + uint32_t d_seed; uint32_t d_shift_register_length; // less than 32 static uint32_t @@ -99,7 +100,10 @@ class gri_lfsr public: gri_lfsr(uint32_t mask, uint32_t seed, uint32_t reg_len) - : d_shift_register(seed), d_mask(mask), d_shift_register_length(reg_len) + : d_shift_register(seed), + d_mask(mask), + d_seed(seed), + d_shift_register_length(reg_len) { if (reg_len > 31) throw std::invalid_argument("reg_len must be <= 31"); @@ -126,6 +130,10 @@ class gri_lfsr return output; } + /*! + * Reset shift register to initial seed value + */ + void reset() { d_shift_register = d_seed; } /*! * Rotate the register through x number of bits diff --git a/gnuradio-core/src/lib/io/Makefile.am b/gnuradio-core/src/lib/io/Makefile.am index 9eacd137d..c52554645 100644 --- a/gnuradio-core/src/lib/io/Makefile.am +++ b/gnuradio-core/src/lib/io/Makefile.am @@ -39,7 +39,6 @@ libio_la_SOURCES = \ gr_oscope_guts.cc \ gr_oscope_sink_f.cc \ gr_oscope_sink_x.cc \ - gri_logger.cc \ i2c.cc \ i2c_bitbang.cc \ i2c_bbio.cc \ @@ -72,7 +71,6 @@ grinclude_HEADERS = \ gr_oscope_sink_f.h \ gr_oscope_sink_x.h \ gr_trigger_mode.h \ - gri_logger.h \ i2c.h \ i2c_bitbang.h \ i2c_bbio.h \ diff --git a/gnuradio-core/src/lib/io/gri_logger.cc b/gnuradio-core/src/lib/io/gri_logger.cc deleted file mode 100644 index 473a7c5ed..000000000 --- a/gnuradio-core/src/lib/io/gri_logger.cc +++ /dev/null @@ -1,178 +0,0 @@ -/* -*- c++ -*- */ -/* - * Copyright 2006,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 - -#if 0 // This needs reimplementation with boost threads and synchronization - -#include <gri_logger.h> -#include <stdio.h> -#include <stdarg.h> -#include <stdexcept> -#include <boost/weak_ptr.hpp> -#include <string.h> - - -/* - * This class creates the thread that reads from the ringbuffer and - * and writes to the file. This is opaque to the user. - */ -class gri_log_poster : public omni_thread -{ - FILE *d_fp; - gr_buffer_sptr d_writer; - gr_buffer_reader_sptr d_reader; - omni_semaphore d_ringbuffer_ready; - volatile bool d_time_to_die; - volatile bool d_writer_overrun; - - virtual void* run_undetached(void * arg); - -public: - gri_log_poster(const char *filename); - ~gri_log_poster(); - - void kill() { d_time_to_die = true; post(); } - gr_buffer_sptr writer() const { return d_writer; } - void post() { d_ringbuffer_ready.post(); } - void note_writer_overrun() { d_writer_overrun = true; } -}; - -gri_log_poster::gri_log_poster(const char *filename) - : omni_thread(), - d_ringbuffer_ready(1, 1), // binary semaphore - d_time_to_die(false), - d_writer_overrun(false) -{ - if ((d_fp = fopen(filename, "w")) == 0){ - perror (filename); - throw std::runtime_error("can't open file"); - } - - // Create a 1MB buffer. - d_writer = gr_make_buffer(1 * 1024 * 1024, sizeof(unsigned char)); - d_reader = gr_buffer_add_reader(d_writer, 0); - - start_undetached(); // start the thread -} - -gri_log_poster::~gri_log_poster() -{ - if (d_fp != 0){ - fclose(d_fp); - d_fp = 0; - } -} - -/* - * This is the body of the logging thread. - */ -void * -gri_log_poster::run_undetached(void *arg) -{ - int nbytes; - - //fprintf(stderr, "Enter: run_undetached!\n"); - - while (!d_time_to_die){ - while ((nbytes = d_reader->items_available()) > 0){ - fwrite(d_reader->read_pointer(), 1, nbytes, d_fp); - d_reader->update_read_pointer(nbytes); - } - fflush(d_fp); - d_ringbuffer_ready.wait(); - - if (d_writer_overrun){ - fputs(">>>>> gri_logger: writer overrun. Info lost <<<<<\n", d_fp); - d_writer_overrun = false; - } - } - - // fprintf(stderr, "Exit: run_undetached!\n"); - return 0; -} - -// ------------------------------------------------------------------------ - -static boost::weak_ptr<gri_logger> s_singleton; // weak pointer IQ test ;-) -static omni_mutex s_singleton_mutex; - -gri_logger_sptr -gri_logger::singleton() -{ - omni_mutex_lock l(s_singleton_mutex); - gri_logger_sptr r; - - if (r = s_singleton.lock()) - return r; - - r = gri_logger_sptr(new gri_logger("gri_logger.log")); - s_singleton = r; - return r; -} - - -gri_logger::gri_logger(const char *filename) -{ - d_poster = new gri_log_poster(filename); -} - -gri_logger::~gri_logger() -{ - d_poster->kill(); - d_poster->join(NULL); -} - -void -gri_logger::write(const void *buf, size_t count) -{ - omni_mutex_lock l(d_write_mutex); - gr_buffer_sptr writer = d_poster->writer(); - - // either write it all, or drop it on the ground - if (count <= (size_t) writer->space_available()){ - memcpy(writer->write_pointer(), buf, count); - writer->update_write_pointer(count); - d_poster->post(); - } - else { - d_poster->note_writer_overrun(); - } -} - -void -gri_logger::printf(const char *format, ...) -{ - va_list ap; - char buf[4096]; - int n; - - va_start(ap, format); - n = vsnprintf(buf, sizeof(buf), format, ap); - va_end(ap); - if (n > -1 && n < (ssize_t) sizeof(buf)) - write(buf, n); -} - -#endif diff --git a/gnuradio-core/src/lib/io/gri_logger.h b/gnuradio-core/src/lib/io/gri_logger.h deleted file mode 100644 index 0a1414540..000000000 --- a/gnuradio-core/src/lib/io/gri_logger.h +++ /dev/null @@ -1,59 +0,0 @@ -/* -*- c++ -*- */ -/* - * Copyright 2006,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_GRI_LOGGER_H -#define INCLUDED_GRI_LOGGER_H - -#if 0 // This needs reimplementation with boost threads and synchronization - -#include <stddef.h> -#include <gnuradio/omnithread.h> -#include <gr_buffer.h> - -class gri_log_poster; -class gri_logger; -typedef boost::shared_ptr<gri_logger> gri_logger_sptr; - - -/*! - * \brief non-blocking logging to a file. - * - * In reality, this may block, but only for a bounded time. - * Trust me, it's safe to use from portaudio and JACK callbacks. - */ -class gri_logger -{ - gri_log_poster *d_poster; - omni_mutex d_write_mutex; - -public: - static gri_logger_sptr singleton(); - - gri_logger(const char *filename); - ~gri_logger(); - - void write(const void *buf, size_t count); - void printf(const char *format, ...); -}; - -#endif - -#endif /* INCLUDED_GRI_LOGGER_H */ diff --git a/gnuradio-core/src/lib/runtime/gr_top_block_impl.cc b/gnuradio-core/src/lib/runtime/gr_top_block_impl.cc index 50d480d00..9cad687fb 100644 --- a/gnuradio-core/src/lib/runtime/gr_top_block_impl.cc +++ b/gnuradio-core/src/lib/runtime/gr_top_block_impl.cc @@ -90,7 +90,7 @@ gr_top_block_impl::~gr_top_block_impl() void gr_top_block_impl::start() { - gr_lock_guard l(d_mutex); + gruel::scoped_lock l(d_mutex); if (d_state != IDLE) throw std::runtime_error("top_block::start: top block already running or wait() not called after previous stop()"); @@ -131,14 +131,14 @@ gr_top_block_impl::wait() void gr_top_block_impl::lock() { - gr_lock_guard lock(d_mutex); + gruel::scoped_lock lock(d_mutex); d_lock_count++; } void gr_top_block_impl::unlock() { - gr_lock_guard lock(d_mutex); + gruel::scoped_lock lock(d_mutex); if (d_lock_count <= 0){ d_lock_count = 0; // fix it, then complain diff --git a/gnuradio-core/src/lib/runtime/gr_top_block_impl.h b/gnuradio-core/src/lib/runtime/gr_top_block_impl.h index 35fb44ef9..ef28dd829 100644 --- a/gnuradio-core/src/lib/runtime/gr_top_block_impl.h +++ b/gnuradio-core/src/lib/runtime/gr_top_block_impl.h @@ -24,10 +24,7 @@ #define INCLUDED_GR_TOP_BLOCK_IMPL_H #include <gr_scheduler.h> -#include <boost/thread.hpp> - -typedef boost::mutex gr_mutex; // FIXME move these elsewhere -typedef boost::lock_guard<boost::mutex> gr_lock_guard; +#include <gruel/thread.h> /*! *\brief Abstract implementation details of gr_top_block @@ -69,7 +66,7 @@ protected: gr_flat_flowgraph_sptr d_ffg; gr_scheduler_sptr d_scheduler; - gr_mutex d_mutex; // protects d_state and d_lock_count + gruel::mutex d_mutex; // protects d_state and d_lock_count tb_state d_state; int d_lock_count; diff --git a/gnuradio-core/src/python/gnuradio/Makefile.am b/gnuradio-core/src/python/gnuradio/Makefile.am index dcc0017b3..f0516f2fd 100644 --- a/gnuradio-core/src/python/gnuradio/Makefile.am +++ b/gnuradio-core/src/python/gnuradio/Makefile.am @@ -30,6 +30,7 @@ grpython_PYTHON = \ eng_notation.py \ eng_option.py \ modulation_utils.py \ + modulation_utils2.py \ ofdm_packet_utils.py \ packet_utils.py \ gr_unittest.py \ diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am b/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am index 68d683623..7b24fb69d 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am +++ b/gnuradio-core/src/python/gnuradio/blks2impl/Makefile.am @@ -1,5 +1,5 @@ # -# Copyright 2005,2007,2009 Free Software Foundation, Inc. +# Copyright 2005,2007,2009,2010 Free Software Foundation, Inc. # # This file is part of GNU Radio # diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk.py b/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk.py index 860015c3f..55e4890f3 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk.py +++ b/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk.py @@ -233,10 +233,10 @@ class dbpsk_demod(gr.hier_block2): arity = pow(2,self.bits_per_symbol()) # Automatic gain control - scale = (1.0/16384.0) - self.pre_scaler = gr.multiply_const_cc(scale) # scale the signal from full-range to +-1 - #self.agc = gr.agc2_cc(0.6e-1, 1e-3, 1, 1, 100) - self.agc = gr.feedforward_agc_cc(16, 2.0) + #scale = (1.0/16384.0) + #self.pre_scaler = gr.multiply_const_cc(scale) # scale the signal from full-range to +-1 + self.agc = gr.agc2_cc(0.6e-1, 1e-3, 1, 1, 100) + #self.agc = gr.feedforward_agc_cc(16, 2.0) # RRC data filter ntaps = 11 * samples_per_symbol @@ -288,7 +288,7 @@ class dbpsk_demod(gr.hier_block2): self._setup_logging() # Connect and Initialize base class - self.connect(self, self.pre_scaler, self.agc, self.rrc_filter, self.receiver, + self.connect(self, self.agc, self.rrc_filter, self.receiver, self.diffdec, self.slicer, self.symbol_mapper, self.unpack, self) def samples_per_symbol(self): diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk2.py b/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk2.py index 4541b453b..d7bcf5390 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk2.py +++ b/gnuradio-core/src/python/gnuradio/blks2impl/dbpsk2.py @@ -1,5 +1,5 @@ # -# Copyright 2005,2006,2007 Free Software Foundation, Inc. +# Copyright 2009,2010 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -25,8 +25,8 @@ differential BPSK modulation and demodulation. """ -from gnuradio import gr, gru, modulation_utils -from math import pi, sqrt +from gnuradio import gr, gru, modulation_utils2 +from math import pi, sqrt, ceil import psk import cmath from pprint import pprint @@ -38,7 +38,8 @@ _def_gray_code = True _def_verbose = False _def_log = False -_def_costas_alpha = 0.1 +_def_freq_alpha = 0.010 +_def_phase_alpha = 0.1 _def_timing_alpha = 0.100 _def_timing_beta = 0.010 _def_timing_max_dev = 1.5 @@ -82,11 +83,9 @@ class dbpsk2_mod(gr.hier_block2): self._excess_bw = excess_bw self._gray_code = gray_code - if not isinstance(self._samples_per_symbol, int) or self._samples_per_symbol < 2: + if self._samples_per_symbol < 2: raise TypeError, ("sbp must be an integer >= 2, is %d" % self._samples_per_symbol) - ntaps = 11 * self._samples_per_symbol - arity = pow(2,self.bits_per_symbol()) # turn bytes into k-bit vectors @@ -103,15 +102,15 @@ class dbpsk2_mod(gr.hier_block2): self.chunks2symbols = gr.chunks_to_symbols_bc(psk.constellation[arity]) # pulse shaping filter - self.rrc_taps = gr.firdes.root_raised_cosine( - self._samples_per_symbol, # gain (samples_per_symbol since we're - # interpolating by samples_per_symbol) - self._samples_per_symbol, # sampling rate - 1.0, # symbol rate - self._excess_bw, # excess bandwidth (roll-off factor) + nfilts = 32 + ntaps = nfilts * 11 * int(self._samples_per_symbol) # make nfilts filters of ntaps each + self.rrc_taps = gr.firdes.root_raised_cosine( + nfilts, # gain + nfilts, # sampling rate based on 32 filters in resampler + 1.0, # symbol rate + self._excess_bw, # excess bandwidth (roll-off factor) ntaps) - self.rrc_filter = gr.interp_fir_filter_ccf(self._samples_per_symbol, - self.rrc_taps) + self.rrc_filter = gr.pfb_arb_resampler_ccf(self._samples_per_symbol, self.rrc_taps) # Connect self.connect(self, self.bytes2chunks, self.symbol_mapper, self.diffenc, @@ -146,7 +145,7 @@ class dbpsk2_mod(gr.hier_block2): """ Given command line options, create dictionary suitable for passing to __init__ """ - return modulation_utils.extract_kwargs_from_options(dbpsk2_mod.__init__, + return modulation_utils2.extract_kwargs_from_options(dbpsk2_mod.__init__, ('self',), options) extract_kwargs_from_options=staticmethod(extract_kwargs_from_options) @@ -182,12 +181,14 @@ class dbpsk2_demod(gr.hier_block2): def __init__(self, samples_per_symbol=_def_samples_per_symbol, excess_bw=_def_excess_bw, - costas_alpha=_def_costas_alpha, + freq_alpha=_def_freq_alpha, + phase_alpha=_def_phase_alpha, timing_alpha=_def_timing_alpha, timing_max_dev=_def_timing_max_dev, gray_code=_def_gray_code, verbose=_def_verbose, - log=_def_log): + log=_def_log, + sync_out=False): """ Hierarchical block for RRC-filtered differential BPSK demodulation @@ -198,9 +199,11 @@ class dbpsk2_demod(gr.hier_block2): @type samples_per_symbol: float @param excess_bw: Root-raised cosine filter excess bandwidth @type excess_bw: float - @param costas_alpha: loop filter gain - @type costas_alpha: float - @param timing_alpha: timing loop alpha gain + @param freq_alpha: loop filter gain for frequency recovery + @type freq_alpha: float + @param phase_alpha: loop filter gain for phase/fine frequency recovery + @type phase_alpha: float + @param timing_alpha: loop alpha gain for timing recovery @type timing_alpha: float @param timing_max: timing loop maximum rate deviations @type timing_max: float @@ -208,17 +211,23 @@ class dbpsk2_demod(gr.hier_block2): @type gray_code: bool @param verbose: Print information about modulator? @type verbose: bool - @param debug: Print modualtion data to files? - @type debug: bool + @param log: Print modualtion data to files? + @type log: bool + @param sync_out: Output a sync signal on :1? + @type sync_out: bool """ - - gr.hier_block2.__init__(self, "dbpsk2_demod", - gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature - gr.io_signature(1, 1, gr.sizeof_char)) # Output signature - + if sync_out: io_sig_out = gr.io_signaturev(2, 2, (gr.sizeof_char, gr.sizeof_gr_complex)) + else: io_sig_out = gr.io_signature(1, 1, gr.sizeof_char) + + gr.hier_block2.__init__(self, "dqpsk2_demod", + gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature + io_sig_out) # Output signature + self._samples_per_symbol = samples_per_symbol self._excess_bw = excess_bw - self._costas_alpha = costas_alpha + self._freq_alpha = freq_alpha + self._freq_beta = 0.10*self._freq_alpha + self._phase_alpha = phase_alpha self._timing_alpha = timing_alpha self._timing_beta = _def_timing_beta self._timing_max_dev=timing_max_dev @@ -233,24 +242,32 @@ class dbpsk2_demod(gr.hier_block2): self.agc = gr.agc2_cc(0.6e-1, 1e-3, 1, 1, 100) #self.agc = gr.feedforward_agc_cc(16, 1.0) - self._costas_beta = 0.25 * self._costas_alpha * self._costas_alpha - # Allow a frequency swing of +/- half of the sample rate - fmin = -0.5 - fmax = 0.5 - - self.clock_recov = gr.costas_loop_cc(self._costas_alpha, - self._costas_beta, - fmax, fmin, arity) + # Frequency correction + self.freq_recov = gr.fll_band_edge_cc(self._samples_per_symbol, self._excess_bw, + 11*int(self._samples_per_symbol), + self._freq_alpha, self._freq_beta) # symbol timing recovery with RRC data filter nfilts = 32 - ntaps = 11 * samples_per_symbol*nfilts - taps = gr.firdes.root_raised_cosine(nfilts, nfilts, 1.0/float(self._samples_per_symbol), self._excess_bw, ntaps) + ntaps = 11 * int(self._samples_per_symbol*nfilts) + taps = gr.firdes.root_raised_cosine(nfilts, nfilts, + 1.0/float(self._samples_per_symbol), + self._excess_bw, ntaps) self.time_recov = gr.pfb_clock_sync_ccf(self._samples_per_symbol, self._timing_alpha, taps, nfilts, nfilts/2, self._timing_max_dev) self.time_recov.set_beta(self._timing_beta) - + + # Perform phase / fine frequency correction + self._phase_beta = 0.25 * self._phase_alpha * self._phase_alpha + # Allow a frequency swing of +/- half of the sample rate + fmin = -0.5 + fmax = 0.5 + + self.phase_recov = gr.costas_loop_cc(self._phase_alpha, + self._phase_beta, + fmax, fmin, arity) + # Do differential decoding based on phase change of symbols self.diffdec = gr.diff_phasor_cc() @@ -275,9 +292,9 @@ class dbpsk2_demod(gr.hier_block2): # Connect self.connect(self, self.agc, - self.clock_recov, - self.time_recov, + self.freq_recov, self.time_recov, self.phase_recov, self.diffdec, self.slicer, self.symbol_mapper, self.unpack, self) + if sync_out: self.connect(self.time_recov, (self, 1)) def samples_per_symbol(self): return self._samples_per_symbol @@ -291,28 +308,27 @@ class dbpsk2_demod(gr.hier_block2): print "bits per symbol: %d" % self.bits_per_symbol() print "Gray code: %s" % self._gray_code print "RRC roll-off factor: %.2f" % self._excess_bw - print "Costas Loop alpha: %.2f" % self._costas_alpha - print "Costas Loop beta: %.2f" % self._costas_beta - print "Timing alpha gain: %.2f" % self._timing_alpha - print "Timing beta gain: %.2f" % self._timing_beta + print "FLL gain: %.2e" % self._freq_alpha + print "Timing alpha gain: %.2e" % self._timing_alpha + print "Timing beta gain: %.2e" % self._timing_beta print "Timing max dev: %.2f" % self._timing_max_dev + print "Phase track alpha: %.2e" % self._phase_alpha + print "Phase track beta: %.2e" % self._phase_beta def _setup_logging(self): print "Modulation logging turned on." - self.connect(self.pre_scaler, - gr.file_sink(gr.sizeof_gr_complex, "rx_prescaler.dat")) self.connect(self.agc, gr.file_sink(gr.sizeof_gr_complex, "rx_agc.dat")) - self.connect(self.rrc_filter, - gr.file_sink(gr.sizeof_gr_complex, "rx_rrc_filter.dat")) - self.connect(self.clock_recov, - gr.file_sink(gr.sizeof_gr_complex, "rx_clock_recov.dat")) + self.connect(self.freq_recov, + gr.file_sink(gr.sizeof_gr_complex, "rx_freq_recov.dat")) self.connect(self.time_recov, gr.file_sink(gr.sizeof_gr_complex, "rx_time_recov.dat")) + self.connect(self.phase_recov, + gr.file_sink(gr.sizeof_gr_complex, "rx_phase_recov.dat")) self.connect(self.diffdec, gr.file_sink(gr.sizeof_gr_complex, "rx_diffdec.dat")) self.connect(self.slicer, - gr.file_sink(gr.sizeof_char, "rx_slicer.dat")) + gr.file_sink(gr.sizeof_char, "rx_slicer.dat")) self.connect(self.symbol_mapper, gr.file_sink(gr.sizeof_char, "rx_symbol_mapper.dat")) self.connect(self.unpack, @@ -327,11 +343,13 @@ class dbpsk2_demod(gr.hier_block2): parser.add_option("", "--no-gray-code", dest="gray_code", action="store_false", default=_def_gray_code, help="disable gray coding on modulated bits (PSK)") - parser.add_option("", "--costas-alpha", type="float", default=None, - help="set Costas loop alpha value [default=%default] (PSK)") - parser.add_option("", "--gain-alpha", type="float", default=_def_timing_alpha, + parser.add_option("", "--freq-alpha", type="float", default=_def_freq_alpha, + help="set frequency lock loop alpha gain value [default=%default] (PSK)") + parser.add_option("", "--phase-alpha", type="float", default=_def_phase_alpha, + help="set phase tracking loop alpha value [default=%default] (PSK)") + parser.add_option("", "--timing-alpha", type="float", default=_def_timing_alpha, help="set timing symbol sync loop gain alpha value [default=%default] (GMSK/PSK)") - parser.add_option("", "--gain-beta", type="float", default=_def_timing_beta, + parser.add_option("", "--timing-beta", type="float", default=_def_timing_beta, help="set timing symbol sync loop gain beta value [default=%default] (GMSK/PSK)") parser.add_option("", "--timing-max-dev", type="float", default=_def_timing_max_dev, help="set timing symbol sync loop maximum deviation [default=%default] (GMSK/PSK)") @@ -341,11 +359,11 @@ class dbpsk2_demod(gr.hier_block2): """ Given command line options, create dictionary suitable for passing to __init__ """ - return modulation_utils.extract_kwargs_from_options( + return modulation_utils2.extract_kwargs_from_options( dbpsk2_demod.__init__, ('self',), options) extract_kwargs_from_options=staticmethod(extract_kwargs_from_options) # # Add these to the mod/demod registry # -modulation_utils.add_type_1_mod('dbpsk2', dbpsk2_mod) -modulation_utils.add_type_1_demod('dbpsk2', dbpsk2_demod) +modulation_utils2.add_type_1_mod('dbpsk2', dbpsk2_mod) +modulation_utils2.add_type_1_demod('dbpsk2', dbpsk2_demod) diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/dqpsk2.py b/gnuradio-core/src/python/gnuradio/blks2impl/dqpsk2.py index 9704ac98c..e1e627707 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/dqpsk2.py +++ b/gnuradio-core/src/python/gnuradio/blks2impl/dqpsk2.py @@ -1,5 +1,5 @@ # -# Copyright 2005,2006,2007,2009 Free Software Foundation, Inc. +# Copyright 2009,2010 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -25,7 +25,7 @@ differential QPSK modulation and demodulation. """ -from gnuradio import gr, gru, modulation_utils +from gnuradio import gr, gru, modulation_utils2 from math import pi, sqrt import psk import cmath @@ -38,7 +38,8 @@ _def_gray_code = True _def_verbose = False _def_log = False -_def_costas_alpha = 0.01 +_def_freq_alpha = 0.010 +_def_phase_alpha = 0.01 _def_timing_alpha = 0.100 _def_timing_beta = 0.010 _def_timing_max_dev = 1.5 @@ -82,8 +83,8 @@ class dqpsk2_mod(gr.hier_block2): self._excess_bw = excess_bw self._gray_code = gray_code - if not isinstance(samples_per_symbol, int) or samples_per_symbol < 2: - raise TypeError, ("sbp must be an integer >= 2, is %d" % samples_per_symbol) + if samples_per_symbol < 2: + raise TypeError, ("sbp must be >= 2, is %f" % samples_per_symbol) ntaps = 11 * samples_per_symbol @@ -105,14 +106,15 @@ class dqpsk2_mod(gr.hier_block2): self.chunks2symbols = gr.chunks_to_symbols_bc(rotated_const) # pulse shaping filter - self.rrc_taps = gr.firdes.root_raised_cosine( - self._samples_per_symbol, # gain (sps since we're interpolating by sps) - self._samples_per_symbol, # sampling rate - 1.0, # symbol rate - self._excess_bw, # excess bandwidth (roll-off factor) + nfilts = 32 + ntaps = 11 * int(nfilts * self._samples_per_symbol) # make nfilts filters of ntaps each + self.rrc_taps = gr.firdes.root_raised_cosine( + nfilts, # gain + nfilts, # sampling rate based on 32 filters in resampler + 1.0, # symbol rate + self._excess_bw, # excess bandwidth (roll-off factor) ntaps) - - self.rrc_filter = gr.interp_fir_filter_ccf(self._samples_per_symbol, self.rrc_taps) + self.rrc_filter = gr.pfb_arb_resampler_ccf(self._samples_per_symbol, self.rrc_taps) if verbose: self._print_verbage() @@ -166,7 +168,7 @@ class dqpsk2_mod(gr.hier_block2): """ Given command line options, create dictionary suitable for passing to __init__ """ - return modulation_utils.extract_kwargs_from_options(dqpsk2_mod.__init__, + return modulation_utils2.extract_kwargs_from_options(dqpsk2_mod.__init__, ('self',), options) extract_kwargs_from_options=staticmethod(extract_kwargs_from_options) @@ -182,12 +184,14 @@ class dqpsk2_demod(gr.hier_block2): def __init__(self, samples_per_symbol=_def_samples_per_symbol, excess_bw=_def_excess_bw, - costas_alpha=_def_costas_alpha, + freq_alpha=_def_freq_alpha, + phase_alpha=_def_phase_alpha, timing_alpha=_def_timing_alpha, timing_max_dev=_def_timing_max_dev, gray_code=_def_gray_code, verbose=_def_verbose, - log=_def_log): + log=_def_log, + sync_out=False): """ Hierarchical block for RRC-filtered DQPSK demodulation @@ -198,8 +202,10 @@ class dqpsk2_demod(gr.hier_block2): @type samples_per_symbol: float @param excess_bw: Root-raised cosine filter excess bandwidth @type excess_bw: float - @param costas_alpha: loop filter gain - @type costas_alphas: float + @param freq_alpha: loop filter gain for frequency recovery + @type freq_alpha: float + @param phase_alpha: loop filter gain + @type phase_alphas: float @param timing_alpha: timing loop alpha gain @type timing_alpha: float @param timing_max: timing loop maximum rate deviations @@ -208,17 +214,23 @@ class dqpsk2_demod(gr.hier_block2): @type gray_code: bool @param verbose: Print information about modulator? @type verbose: bool - @param debug: Print modualtion data to files? - @type debug: bool + @param log: Print modualtion data to files? + @type log: bool + @param sync_out: Output a sync signal on :1? + @type sync_out: bool """ + if sync_out: io_sig_out = gr.io_signaturev(2, 2, (gr.sizeof_char, gr.sizeof_gr_complex)) + else: io_sig_out = gr.io_signature(1, 1, gr.sizeof_char) gr.hier_block2.__init__(self, "dqpsk2_demod", gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature - gr.io_signature(1, 1, gr.sizeof_char)) # Output signature + io_sig_out) # Output signature self._samples_per_symbol = samples_per_symbol self._excess_bw = excess_bw - self._costas_alpha = costas_alpha + self._freq_alpha = freq_alpha + self._freq_beta = 0.25*self._freq_alpha**2 + self._phase_alpha = phase_alpha self._timing_alpha = timing_alpha self._timing_beta = _def_timing_beta self._timing_max_dev=timing_max_dev @@ -233,24 +245,35 @@ class dqpsk2_demod(gr.hier_block2): self.agc = gr.agc2_cc(0.6e-1, 1e-3, 1, 1, 100) #self.agc = gr.feedforward_agc_cc(16, 2.0) - self._costas_beta = 0.25 * self._costas_alpha * self._costas_alpha - # Allow a frequency swing of +/- half of the sample rate - fmin = -0.5 - fmax = 0.5 + # Frequency correction + self.freq_recov = gr.fll_band_edge_cc(self._samples_per_symbol, self._excess_bw, + 11*int(self._samples_per_symbol), + self._freq_alpha, self._freq_beta) - self.clock_recov = gr.costas_loop_cc(self._costas_alpha, - self._costas_beta, - fmax, fmin, arity) # symbol timing recovery with RRC data filter nfilts = 32 - ntaps = 11 * samples_per_symbol*nfilts - taps = gr.firdes.root_raised_cosine(nfilts, nfilts, 1.0/float(self._samples_per_symbol), self._excess_bw, ntaps) + ntaps = 11 * int(samples_per_symbol*nfilts) + taps = gr.firdes.root_raised_cosine(nfilts, nfilts, + 1.0/float(self._samples_per_symbol), + self._excess_bw, ntaps) self.time_recov = gr.pfb_clock_sync_ccf(self._samples_per_symbol, self._timing_alpha, taps, nfilts, nfilts/2, self._timing_max_dev) self.time_recov.set_beta(self._timing_beta) + + # Perform phase / fine frequency correction + self._phase_beta = 0.25 * self._phase_alpha * self._phase_alpha + # Allow a frequency swing of +/- half of the sample rate + fmin = -0.5 + fmax = 0.5 + + self.phase_recov = gr.costas_loop_cc(self._phase_alpha, + self._phase_beta, + fmax, fmin, arity) + + # Perform Differential decoding on the constellation self.diffdec = gr.diff_phasor_cc() @@ -275,9 +298,9 @@ class dqpsk2_demod(gr.hier_block2): # Connect self.connect(self, self.agc, - self.clock_recov, - self.time_recov, + self.freq_recov, self.time_recov, self.phase_recov, self.diffdec, self.slicer, self.symbol_mapper, self.unpack, self) + if sync_out: self.connect(self.time_recov, (self, 1)) def samples_per_symbol(self): return self._samples_per_symbol @@ -291,24 +314,23 @@ class dqpsk2_demod(gr.hier_block2): print "bits per symbol: %d" % self.bits_per_symbol() print "Gray code: %s" % self._gray_code print "RRC roll-off factor: %.2f" % self._excess_bw - print "Costas Loop alpha: %.2e" % self._costas_alpha - print "Costas Loop beta: %.2e" % self._costas_beta + print "FLL gain: %.2f" % self._freq_alpha print "Timing alpha gain: %.2f" % self._timing_alpha print "Timing beta gain: %.2f" % self._timing_beta print "Timing max dev: %.2f" % self._timing_max_dev + print "Phase track alpha: %.2e" % self._phase_alpha + print "Phase track beta: %.2e" % self._phase_beta def _setup_logging(self): print "Modulation logging turned on." - self.connect(self.pre_scaler, - gr.file_sink(gr.sizeof_gr_complex, "rx_prescaler.dat")) self.connect(self.agc, gr.file_sink(gr.sizeof_gr_complex, "rx_agc.dat")) - self.connect(self.rrc_filter, - gr.file_sink(gr.sizeof_gr_complex, "rx_rrc_filter.dat")) - self.connect(self.clock_recov, - gr.file_sink(gr.sizeof_gr_complex, "rx_clock_recov.dat")) + self.connect(self.freq_recov, + gr.file_sink(gr.sizeof_gr_complex, "rx_freq_recov.dat")) self.connect(self.time_recov, gr.file_sink(gr.sizeof_gr_complex, "rx_time_recov.dat")) + self.connect(self.phase_recov, + gr.file_sink(gr.sizeof_gr_complex, "rx_phase_recov.dat")) self.connect(self.diffdec, gr.file_sink(gr.sizeof_gr_complex, "rx_diffdec.dat")) self.connect(self.slicer, @@ -320,18 +342,20 @@ class dqpsk2_demod(gr.hier_block2): def add_options(parser): """ - Adds modulation-specific options to the standard parser + Adds DQPSK demodulation-specific options to the standard parser """ parser.add_option("", "--excess-bw", type="float", default=_def_excess_bw, help="set RRC excess bandwith factor [default=%default] (PSK)") parser.add_option("", "--no-gray-code", dest="gray_code", action="store_false", default=_def_gray_code, help="disable gray coding on modulated bits (PSK)") - parser.add_option("", "--costas-alpha", type="float", default=_def_costas_alpha, - help="set Costas loop alpha value [default=%default] (PSK)") - parser.add_option("", "--gain-alpha", type="float", default=_def_timing_alpha, + parser.add_option("", "--freq-alpha", type="float", default=_def_freq_alpha, + help="set frequency lock loop alpha gain value [default=%default] (PSK)") + parser.add_option("", "--phase-alpha", type="float", default=_def_phase_alpha, + help="set phase tracking loop alpha value [default=%default] (PSK)") + parser.add_option("", "--timing-alpha", type="float", default=_def_timing_alpha, help="set timing symbol sync loop gain alpha value [default=%default] (GMSK/PSK)") - parser.add_option("", "--gain-beta", type="float", default=_def_timing_beta, + parser.add_option("", "--timing-beta", type="float", default=_def_timing_beta, help="set timing symbol sync loop gain beta value [default=%default] (GMSK/PSK)") parser.add_option("", "--timing-max-dev", type="float", default=_def_timing_max_dev, help="set timing symbol sync loop maximum deviation [default=%default] (GMSK/PSK)") @@ -341,7 +365,7 @@ class dqpsk2_demod(gr.hier_block2): """ Given command line options, create dictionary suitable for passing to __init__ """ - return modulation_utils.extract_kwargs_from_options( + return modulation_utils2.extract_kwargs_from_options( dqpsk2_demod.__init__, ('self',), options) extract_kwargs_from_options=staticmethod(extract_kwargs_from_options) @@ -349,5 +373,5 @@ class dqpsk2_demod(gr.hier_block2): # # Add these to the mod/demod registry # -modulation_utils.add_type_1_mod('dqpsk2', dqpsk2_mod) -modulation_utils.add_type_1_demod('dqpsk2', dqpsk2_demod) +modulation_utils2.add_type_1_mod('dqpsk2', dqpsk2_mod) +modulation_utils2.add_type_1_demod('dqpsk2', dqpsk2_demod) diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/generic_usrp.py b/gnuradio-core/src/python/gnuradio/blks2impl/generic_usrp.py index 5abbaf9eb..82d1eca13 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/generic_usrp.py +++ b/gnuradio-core/src/python/gnuradio/blks2impl/generic_usrp.py @@ -121,12 +121,6 @@ class _generic_usrp_base(object): def set_auto_tr(self, enable): if self._type == USRP1_TYPE: return self._subdev.set_auto_tr(enable) - def __del__(self): - try: # Avoid weak reference error - del self._u - del self._subdev - except: pass - ######################################################################## # generic usrp source ######################################################################## diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/logpwrfft.py b/gnuradio-core/src/python/gnuradio/blks2impl/logpwrfft.py index 7ef40be40..200c4cfbe 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/logpwrfft.py +++ b/gnuradio-core/src/python/gnuradio/blks2impl/logpwrfft.py @@ -28,7 +28,7 @@ class _logpwrfft_base(gr.hier_block2): Create a log10(abs(fft)) stream chain, with real or complex input. """ - def __init__(self, sample_rate, fft_size, ref_scale, frame_rate, avg_alpha, average): + def __init__(self, sample_rate, fft_size, ref_scale, frame_rate, avg_alpha, average, win=None): """ Create an log10(abs(fft)) stream chain. Provide access to the setting the filter and sample rate. @@ -38,6 +38,7 @@ class _logpwrfft_base(gr.hier_block2): @param frame_rate Output frame rate @param avg_alpha FFT averaging (over time) constant [0.0-1.0] @param average Whether to average [True, False] + @param win the window taps generation function """ gr.hier_block2.__init__(self, self._name, gr.io_signature(1, 1, self._item_size), # Input signature @@ -48,7 +49,8 @@ class _logpwrfft_base(gr.hier_block2): vec_rate=frame_rate, vec_len=fft_size) - fft_window = window.blackmanharris(fft_size) + if win is None: win = window.blackmanharris + fft_window = win(fft_size) fft = self._fft_block[0](fft_size, True, fft_window) window_power = sum(map(lambda x: x*x, fft_window)) diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py index c45ae4d1a..a479ed48e 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py +++ b/gnuradio-core/src/python/gnuradio/blks2impl/pfb_channelizer.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright 2009 Free Software Foundation, Inc. +# Copyright 2009,2010 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -29,16 +29,18 @@ class pfb_channelizer_ccf(gr.hier_block2): 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): + def __init__(self, numchans, taps, oversample_rate=1): 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._oversample_rate = oversample_rate self.s2ss = gr.stream_to_streams(gr.sizeof_gr_complex, self._numchans) - self.pfb = gr.pfb_channelizer_ccf(self._numchans, self._taps) + self.pfb = gr.pfb_channelizer_ccf(self._numchans, self._taps, + self._oversample_rate) self.v2s = gr.vector_to_streams(gr.sizeof_gr_complex, self._numchans) self.connect(self, self.s2ss) diff --git a/gnuradio-core/src/python/gnuradio/blks2impl/stream_to_vector_decimator.py b/gnuradio-core/src/python/gnuradio/blks2impl/stream_to_vector_decimator.py index 39c8b5050..8f75729c9 100644 --- a/gnuradio-core/src/python/gnuradio/blks2impl/stream_to_vector_decimator.py +++ b/gnuradio-core/src/python/gnuradio/blks2impl/stream_to_vector_decimator.py @@ -72,8 +72,7 @@ class stream_to_vector_decimator(gr.hier_block2): self.one_in_n.set_n(self._decim) def _update_decimator(self): - self._decim = max(1, int(round(self._sample_rate/self._vec_len/self._vec_rate))) - self.one_in_n.set_n(self._decim) + self.set_decimation(self._sample_rate/self._vec_len/self._vec_rate) def decimation(self): """ diff --git a/gnuradio-core/src/python/gnuradio/gr/qa_scrambler.py b/gnuradio-core/src/python/gnuradio/gr/qa_scrambler.py index 76b0e62fa..aecf49293 100755 --- a/gnuradio-core/src/python/gnuradio/gr/qa_scrambler.py +++ b/gnuradio-core/src/python/gnuradio/gr/qa_scrambler.py @@ -40,5 +40,25 @@ class test_scrambler(gr_unittest.TestCase): self.tb.run() self.assertEqual(tuple(src_data[:-8]), dst.data()[8:]) # skip garbage during synchronization + def test_additive_scrambler(self): + src_data = (1,)*1000 + src = gr.vector_source_b(src_data, False) + scrambler = gr.additive_scrambler_bb(0x8a, 0x7f, 7) + descrambler = gr.additive_scrambler_bb(0x8a, 0x7f, 7) + dst = gr.vector_sink_b() + self.tb.connect(src, scrambler, descrambler, dst) + self.tb.run() + self.assertEqual(src_data, dst.data()) + + def test_additive_scrambler_reset(self): + src_data = (1,)*1000 + src = gr.vector_source_b(src_data, False) + scrambler = gr.additive_scrambler_bb(0x8a, 0x7f, 7, 100) + descrambler = gr.additive_scrambler_bb(0x8a, 0x7f, 7, 100) + dst = gr.vector_sink_b() + self.tb.connect(src, scrambler, descrambler, dst) + self.tb.run() + self.assertEqual(src_data, dst.data()) + if __name__ == '__main__': gr_unittest.main () diff --git a/gnuradio-core/src/python/gnuradio/gruimpl/hexint.py b/gnuradio-core/src/python/gnuradio/gruimpl/hexint.py index 8d46e8192..f2808c448 100644 --- a/gnuradio-core/src/python/gnuradio/gruimpl/hexint.py +++ b/gnuradio-core/src/python/gnuradio/gruimpl/hexint.py @@ -30,3 +30,15 @@ def hexint(mask): if mask >= 2**31: return int(mask-2**32) return mask + +def hexshort(mask): + """ + Convert unsigned masks into signed shorts. + + This allows us to use hex constants like 0x8000 when talking to + our hardware and not get screwed by them getting treated as python + longs. + """ + if mask >= 2**15: + return int(mask-2**16) + return mask diff --git a/gnuradio-core/src/python/gnuradio/modulation_utils2.py b/gnuradio-core/src/python/gnuradio/modulation_utils2.py new file mode 100644 index 000000000..c5dba3e79 --- /dev/null +++ b/gnuradio-core/src/python/gnuradio/modulation_utils2.py @@ -0,0 +1,81 @@ +# +# Copyright 2010 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 this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +# + +""" +Miscellaneous utilities for managing mods and demods, as well as other items +useful in dealing with generalized handling of different modulations and demods. +""" + +import inspect + + +# Type 1 modulators accept a stream of bytes on their input and produce complex baseband output +_type_1_modulators = {} + +def type_1_mods(): + return _type_1_modulators + +def add_type_1_mod(name, mod_class): + _type_1_modulators[name] = mod_class + + +# Type 1 demodulators accept complex baseband input and produce a stream of bits, packed +# 1 bit / byte as their output. Their output is completely unambiguous. There is no need +# to resolve phase or polarity ambiguities. +_type_1_demodulators = {} + +def type_1_demods(): + return _type_1_demodulators + +def add_type_1_demod(name, demod_class): + _type_1_demodulators[name] = demod_class + + +def extract_kwargs_from_options(function, excluded_args, options): + """ + Given a function, a list of excluded arguments and the result of + parsing command line options, create a dictionary of key word + arguments suitable for passing to the function. The dictionary + will be populated with key/value pairs where the keys are those + that are common to the function's argument list (minus the + excluded_args) and the attributes in options. The values are the + corresponding values from options unless that value is None. + In that case, the corresponding dictionary entry is not populated. + + (This allows different modulations that have the same parameter + names, but different default values to coexist. The downside is + that --help in the option parser will list the default as None, + but in that case the default provided in the __init__ argument + list will be used since there is no kwargs entry.) + + @param function: the function whose parameter list will be examined + @param excluded_args: function arguments that are NOT to be added to the dictionary + @type excluded_args: sequence of strings + @param options: result of command argument parsing + @type options: optparse.Values + """ + # Try this in C++ ;) + args, varargs, varkw, defaults = inspect.getargspec(function) + d = {} + for kw in [a for a in args if a not in excluded_args]: + if hasattr(options, kw): + if getattr(options, kw) is not None: + d[kw] = getattr(options, kw) + return d diff --git a/gnuradio-core/src/python/gnuradio/packet_utils.py b/gnuradio-core/src/python/gnuradio/packet_utils.py index 1417c17fa..e36b05413 100644 --- a/gnuradio-core/src/python/gnuradio/packet_utils.py +++ b/gnuradio-core/src/python/gnuradio/packet_utils.py @@ -143,7 +143,7 @@ def make_packet(payload, samples_per_symbol, bits_per_symbol, (payload_with_crc), '\x55')) if pad_for_usrp: - pkt = pkt + (_npadding_bytes(len(pkt), samples_per_symbol, bits_per_symbol) * '\x55') + pkt = pkt + (_npadding_bytes(len(pkt), int(samples_per_symbol), bits_per_symbol) * '\x55') #print "make_packet: len(pkt) =", len(pkt) return pkt diff --git a/gnuradio-core/src/python/gnuradio/window.py b/gnuradio-core/src/python/gnuradio/window.py index fb4a10675..e109a9892 100644 --- a/gnuradio-core/src/python/gnuradio/window.py +++ b/gnuradio-core/src/python/gnuradio/window.py @@ -1,5 +1,5 @@ # -# Copyright 2004,2005 Free Software Foundation, Inc. +# Copyright 2004,2005,2009 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -153,32 +153,6 @@ def riemann(fft_size): j -= 1 return window -def blackmanharris(fft_size): - a0 = 0.35875 - a1 = 0.48829 - a2 = 0.14128 - a3 = 0.01168 - window = [0 for i in range(fft_size)] - for index in xrange(fft_size): - window[index] = a0 - window[index] -= a1*math.cos(2.0*math.pi*(index+0.5)/(fft_size - 1)) - window[index] += a2*math.cos(4.0*math.pi*(index+0.5)/(fft_size - 1)) - window[index] -= a3*math.cos(6.0*math.pi*(index+0.5)/(fft_size - 1)) - return window - -def nuttall(fft_size): - a0 = 0.3635819 - a1 = 0.4891775 - a2 = 0.1365995 - a3 = 0.0106411 - window = [0 for i in range(fft_size)] - for index in xrange(fft_size): - window[index] = a0 - window[index] -= a1*math.cos(2.0*math.pi*(index+0.5)/(fft_size - 1)) - window[index] += a2*math.cos(4.0*math.pi*(index+0.5)/(fft_size - 1)) - window[index] -= a3*math.cos(6.0*math.pi*(index+0.5)/(fft_size - 1)) - return window - def kaiser(fft_size,beta): ibeta = 1.0/izero(beta) inm1 = 1.0/(fft_size) @@ -187,4 +161,20 @@ def kaiser(fft_size,beta): window[index] = izero(beta*math.sqrt(1.0 - (index * inm1)*(index * inm1))) * ibeta return window - +# Closure to generate functions to create cos windows + +def coswindow(coeffs): + def closure(fft_size): + window = [0] * fft_size + #print list(enumerate(coeffs)) + for w_index in range(fft_size): + for (c_index, coeff) in enumerate(coeffs): + window[w_index] += (-1)**c_index * coeff * math.cos(2.0*c_index*math.pi*(w_index+0.5)/(fft_size-1)) + return window + return closure + +blackmanharris = coswindow((0.35875,0.48829,0.14128,0.01168)) +nuttall = coswindow((0.3635819,0.4891775,0.1365995,0.0106411)) # Wikipedia calls this Blackman-Nuttall +nuttall_cfd = coswindow((0.355768,0.487396,0.144232,0.012604)) # Wikipedia calls this Nuttall, continuous first deriv +flattop = coswindow((1.0,1.93,1.29,0.388,0.032)) # Flat top window, coeffs from Wikipedia +rectangular = lambda fft_size: [1]*fft_size |