diff options
Diffstat (limited to 'gnuradio-core/src/lib')
43 files changed, 2546 insertions, 644 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; |