diff options
35 files changed, 1338 insertions, 308 deletions
diff --git a/gr-blocks/grc/blocks_null_sink.xml b/gr-blocks/grc/blocks_null_sink.xml index 2ae20e619a..f2907d20bd 100644 --- a/gr-blocks/grc/blocks_null_sink.xml +++ b/gr-blocks/grc/blocks_null_sink.xml @@ -45,10 +45,27 @@ <value>1</value> <type>int</type> </param> + <param> + <name>Num Inputs</name> + <key>num_inputs</key> + <value>1</value> + <type>int</type> + </param> + <param> + <name>Bus Connections</name> + <key>bus_conns</key> + <value>[[0,],]</value> + <type>raw</type> + <hide>part</hide> + </param> + <check>$num_inputs >= 1</check> <check>$vlen > 0</check> <sink> <name>in</name> <type>$type</type> <vlen>$vlen</vlen> + <nports>$num_inputs</nports> </sink> + <bus_structure_sink>$bus_conns</bus_structure_sink> </block> + diff --git a/gr-blocks/grc/blocks_null_source.xml b/gr-blocks/grc/blocks_null_source.xml index 01d3905cab..9c109e651f 100644 --- a/gr-blocks/grc/blocks_null_source.xml +++ b/gr-blocks/grc/blocks_null_source.xml @@ -45,10 +45,26 @@ <value>1</value> <type>int</type> </param> + <param> + <name>Num Outputs</name> + <key>num_outputs</key> + <value>1</value> + <type>int</type> + </param> + <param> + <name>Bus Connections</name> + <key>bus_conns</key> + <value>[[0,],]</value> + <type>raw</type> + <hide>part</hide> + </param> + <check>$num_outputs >= 1</check> <check>$vlen > 0</check> <source> <name>out</name> <type>$type</type> <vlen>$vlen</vlen> + <nports>$num_outputs</nports> </source> + <bus_structure_source>$bus_conns</bus_structure_source> </block> diff --git a/gr-blocks/grc/blocks_sub_xx.xml b/gr-blocks/grc/blocks_sub_xx.xml index ae01cf74aa..88b5ccbb57 100644 --- a/gr-blocks/grc/blocks_sub_xx.xml +++ b/gr-blocks/grc/blocks_sub_xx.xml @@ -48,7 +48,7 @@ <type>int</type> </param> <check>$vlen > 0</check> - <check>$num_inputs >= 2</check> + <check>$num_inputs >= 1</check> <sink> <name>in</name> <type>$type</type> diff --git a/gr-blocks/lib/null_sink_impl.cc b/gr-blocks/lib/null_sink_impl.cc index 41adeea0fd..14dd5eff24 100644 --- a/gr-blocks/lib/null_sink_impl.cc +++ b/gr-blocks/lib/null_sink_impl.cc @@ -39,8 +39,8 @@ namespace gr { null_sink_impl::null_sink_impl(size_t sizeof_stream_item) : sync_block("null_sink", - io_signature::make(1, 1, sizeof_stream_item), - io_signature::make(0, 0, 0)) + io_signature::make(1, -1, sizeof_stream_item), + io_signature::make(0, 0, 0)) { } diff --git a/gr-blocks/lib/null_source_impl.cc b/gr-blocks/lib/null_source_impl.cc index edf0104da1..9550dd6bce 100644 --- a/gr-blocks/lib/null_source_impl.cc +++ b/gr-blocks/lib/null_source_impl.cc @@ -41,7 +41,7 @@ namespace gr { null_source_impl::null_source_impl (size_t sizeof_stream_item) : sync_block("null_source", io_signature::make(0, 0, 0), - io_signature::make(1, 1, sizeof_stream_item)) + io_signature::make(1, -1, sizeof_stream_item)) { } @@ -54,8 +54,11 @@ namespace gr { gr_vector_const_void_star &input_items, gr_vector_void_star &output_items) { - void *optr = (void*)output_items[0]; - memset(optr, 0, noutput_items * output_signature()->sizeof_stream_item(0)); + void *optr; + for(size_t n = 0; n < input_items.size(); n++) { + optr = (void*)output_items[n]; + memset(optr, 0, noutput_items * output_signature()->sizeof_stream_item(n)); + } return noutput_items; } diff --git a/gr-filter/examples/reconstruction.py b/gr-filter/examples/reconstruction.py index fd8ba87a2f..0a83b5a4e0 100755 --- a/gr-filter/examples/reconstruction.py +++ b/gr-filter/examples/reconstruction.py @@ -72,7 +72,7 @@ def main(): src = blocks.vector_source_b(data.astype(scipy.uint8).tolist(), False) mod = digital.bpsk_mod(samples_per_symbol=2) - chan = filter.channel_model(npwr) + chan = channels.channel_model(npwr) rrc = filter.fft_filter_ccc(1, rrc_taps) # Split it up into pieces @@ -122,7 +122,7 @@ def main(): s12.set_title("Original Signal in Time") start = 1 - skip = 4 + skip = 2 s13 = f1.add_subplot(2,2,3) s13.plot(sin.real[start::skip], sin.imag[start::skip], "o") s13.set_title("Constellation") @@ -153,7 +153,7 @@ def main(): s32.plot(sout.imag[1000:1500], "o-r") s32.set_title("Reconstructed Signal in Time") - start = 2 + start = 0 skip = 4 s33 = f3.add_subplot(2,2,3) s33.plot(sout.real[start::skip], sout.imag[start::skip], "o") diff --git a/gr-filter/grc/filter_fft_filter_xxx.xml b/gr-filter/grc/filter_fft_filter_xxx.xml index 99b02cbab8..f843af2024 100644 --- a/gr-filter/grc/filter_fft_filter_xxx.xml +++ b/gr-filter/grc/filter_fft_filter_xxx.xml @@ -26,11 +26,18 @@ self.$(id).declare_sample_delay($samp_delay) <opt>taps:complex_vector</opt> </option> <option> + <name>Complex->Complex (Real Taps)</name> + <key>ccf</key> + <opt>input:complex</opt> + <opt>output:complex</opt> + <opt>taps:float_vector</opt> + </option> + <option> <name>Float->Float (Real Taps)</name> <key>fff</key> <opt>input:float</opt> <opt>output:float</opt> - <opt>taps:real_vector</opt> + <opt>taps:float_vector</opt> </option> </param> <param> diff --git a/gr-filter/grc/filter_pfb_channelizer.xml b/gr-filter/grc/filter_pfb_channelizer.xml index 26349cd10d..f8a51d2204 100644 --- a/gr-filter/grc/filter_pfb_channelizer.xml +++ b/gr-filter/grc/filter_pfb_channelizer.xml @@ -16,14 +16,13 @@ $atten) self.$(id).set_channel_map($ch_map) </make> - <!-- Set taps not implemented yet - <callback>set_taps($taps)</callback> - --> + <callback>set_taps($taps)</callback> <callback>set_channel_map($ch_map)</callback> <param> <name>Channels</name> <key>nchans</key> + <value>1</value> <type>int</type> </param> <param> @@ -50,6 +49,13 @@ self.$(id).set_channel_map($ch_map) <value>[]</value> <type>int_vector</type> </param> + <param> + <name>Bus Connections</name> + <key>bus_conns</key> + <value>[[0,],]</value> + <type>raw</type> + <hide>part</hide> + </param> <sink> <name>in</name> <type>complex</type> @@ -59,4 +65,5 @@ self.$(id).set_channel_map($ch_map) <type>complex</type> <nports>$nchans</nports> </source> + <bus_structure_source>$bus_conns</bus_structure_source> </block> diff --git a/gr-filter/grc/filter_pfb_decimator.xml b/gr-filter/grc/filter_pfb_decimator.xml index d57119636c..3a828762a0 100644 --- a/gr-filter/grc/filter_pfb_decimator.xml +++ b/gr-filter/grc/filter_pfb_decimator.xml @@ -13,7 +13,9 @@ $decim, $taps, $channel, - $atten) + $atten, + $fft_rot, + $fft_filts) </make> <callback>set_taps($taps)</callback> <callback>set_channel(int($channel))</callback> @@ -40,6 +42,36 @@ <value>100</value> <type>real</type> </param> + <param> + <name>Use FFT Rotator</name> + <key>fft_rot</key> + <value>True</value> + <type>raw</type> + <hide>part</hide> + <option> + <name>True</name> + <key>True</key> + </option> + <option> + <name>False</name> + <key>False</key> + </option> + </param> + <param> + <name>Use FFT Filters</name> + <key>fft_filts</key> + <value>True</value> + <type>raw</type> + <hide>part</hide> + <option> + <name>True</name> + <key>True</key> + </option> + <option> + <name>False</name> + <key>False</key> + </option> + </param> <sink> <name>in</name> <type>complex</type> diff --git a/gr-filter/grc/filter_pfb_synthesizer.xml b/gr-filter/grc/filter_pfb_synthesizer.xml index e84b25e62e..e7e1ae3951 100644 --- a/gr-filter/grc/filter_pfb_synthesizer.xml +++ b/gr-filter/grc/filter_pfb_synthesizer.xml @@ -45,6 +45,13 @@ self.$(id).set_channel_map($ch_map) <value>[]</value> <type>int_vector</type> </param> + <param> + <name>Bus Connections</name> + <key>bus_conns</key> + <value>[[0,],]</value> + <type>raw</type> + <hide>part</hide> + </param> <sink> <name>in</name> <type>complex</type> @@ -54,4 +61,5 @@ self.$(id).set_channel_map($ch_map) <name>out</name> <type>complex</type> </source> + <bus_structure_sink>$bus_conns</bus_structure_sink> </block> diff --git a/gr-filter/include/gnuradio/filter/CMakeLists.txt b/gr-filter/include/gnuradio/filter/CMakeLists.txt index d039af9b00..9dbc227887 100644 --- a/gr-filter/include/gnuradio/filter/CMakeLists.txt +++ b/gr-filter/include/gnuradio/filter/CMakeLists.txt @@ -94,6 +94,7 @@ install(FILES dc_blocker_ff.h filter_delay_fc.h fft_filter_ccc.h + fft_filter_ccf.h fft_filter_fff.h fractional_interpolator_cc.h fractional_interpolator_ff.h diff --git a/gr-filter/include/gnuradio/filter/fft_filter.h b/gr-filter/include/gnuradio/filter/fft_filter.h index 76271d9abf..24f6ce3134 100644 --- a/gr-filter/include/gnuradio/filter/fft_filter.h +++ b/gr-filter/include/gnuradio/filter/fft_filter.h @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2010,2012 Free Software Foundation, Inc. + * Copyright 2010,2012,2014 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -49,7 +49,7 @@ namespace gr { std::vector<float> d_tail; // state carried between blocks for overlap-add std::vector<float> d_taps; // stores time domain taps gr_complex *d_xformed_taps; // Fourier xformed taps - + void compute_sizes(int ntaps); int tailsize() const { return d_ntaps - 1; } @@ -77,7 +77,7 @@ namespace gr { * \param taps The filter taps (complex) */ int set_taps(const std::vector<float> &taps); - + /*! * \brief Set number of threads to use. */ @@ -92,12 +92,12 @@ namespace gr { * \brief Returns the number of taps in the filter. */ unsigned int ntaps() const; - + /*! * \brief Get number of threads being used. */ int nthreads() const; - + /*! l * \brief Perform the filter operation * @@ -108,7 +108,7 @@ l * \brief Perform the filter operation int filter(int nitems, const float *input, float *output); }; - + /*! * \brief Fast FFT filter with gr_complex input, gr_complex output and gr_complex taps * \ingroup filter_blk @@ -126,7 +126,7 @@ l * \brief Perform the filter operation std::vector<gr_complex> d_tail; // state carried between blocks for overlap-add std::vector<gr_complex> d_taps; // stores time domain taps gr_complex *d_xformed_taps; // Fourier xformed taps - + void compute_sizes(int ntaps); int tailsize() const { return d_ntaps - 1; } @@ -154,12 +154,12 @@ l * \brief Perform the filter operation * \param taps The filter taps (complex) */ int set_taps(const std::vector<gr_complex> &taps); - + /*! * \brief Set number of threads to use. */ void set_nthreads(int n); - + /*! * \brief Returns the taps. */ @@ -169,12 +169,100 @@ l * \brief Perform the filter operation * \brief Returns the number of taps in the filter. */ unsigned int ntaps() const; - + + /*! + * \brief Get number of threads being used. + */ + int nthreads() const; + + /*! + * \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); + }; + + + + /*! + * \brief Fast FFT filter with gr_complex input, gr_complex output and float taps + * \ingroup filter_blk + */ + class FILTER_API fft_filter_ccf + { + private: + int d_ntaps; + int d_nsamples; + int d_fftsize; // fftsize = ntaps + nsamples - 1 + int d_decimation; + fft::fft_complex *d_fwdfft; // forward "plan" + fft::fft_complex *d_invfft; // inverse "plan" + int d_nthreads; // number of FFTW threads to use + std::vector<gr_complex> d_tail; // state carried between blocks for overlap-add + std::vector<float> d_taps; // stores time domain taps + gr_complex *d_xformed_taps; // Fourier xformed 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 fft_filter_ccf). + * + * \param decimation The decimation rate of the filter (int) + * \param taps The filter taps (complex) + * \param nthreads The number of threads for the FFT to use (int) + */ + fft_filter_ccf(int decimation, + const std::vector<float> &taps, + int nthreads=1); + + ~fft_filter_ccf(); + + /*! + * \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<float> &taps); + + /*! + * \brief Set number of threads to use. + */ + void set_nthreads(int n); + + /*! + * \brief Returns the taps. + */ + std::vector<float> taps() const; + + /*! + * \brief Returns the number of taps in the filter. + */ + unsigned int ntaps() const; + + /*! + * \brief Returns the actual size of the filter. + * + * \details This value could be equal to ntaps, but we ofter + * build a longer filter to allow us to calculate a more + * efficient FFT. This value is the actual size of the filters + * used in the calculation of the overlap-and-save operation. + */ + unsigned int filtersize() const; + /*! * \brief Get number of threads being used. */ int nthreads() const; - + /*! * \brief Perform the filter operation * diff --git a/gr-filter/include/gnuradio/filter/fft_filter_ccf.h b/gr-filter/include/gnuradio/filter/fft_filter_ccf.h new file mode 100644 index 0000000000..b0230f8263 --- /dev/null +++ b/gr-filter/include/gnuradio/filter/fft_filter_ccf.h @@ -0,0 +1,89 @@ +/* -*- c++ -*- */ +/* + * Copyright 2014 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_FILTER_FFT_FILTER_CCF_H +#define INCLUDED_FILTER_FFT_FILTER_CCF_H + +#include <gnuradio/filter/api.h> +#include <gnuradio/sync_decimator.h> + +namespace gr { + namespace filter { + + /*! + * \brief Fast FFT filter with gr_complex input, gr_complex output and float taps + * \ingroup filter_blk + * + * \details + * This block implements a complex decimating filter using the + * fast convolution method via an FFT. The decimation factor is an + * interger that is greater than or equal to 1. + * + * The filter takes a set of complex (or real) taps to use in the + * filtering operation. These taps can be defined as anything that + * satisfies the user's filtering needs. For standard filters such + * as lowpass, highpass, bandpass, etc., the filter.firdes and + * filter.optfir classes provide convenient generating methods. + * + * This filter is implemented by using the FFTW package to perform + * the required FFTs. An optional argument, nthreads, may be + * passed to the constructor (or set using the set_nthreads member + * function) to split the FFT among N number of threads. This can + * improve performance on very large FFTs (that is, if the number + * of taps used is very large) if you have enough threads/cores to + * support it. + */ + class FILTER_API fft_filter_ccf : virtual public sync_decimator + { + public: + // gr::filter::fft_filter_ccf::sptr + typedef boost::shared_ptr<fft_filter_ccf> sptr; + + /*! + * Build an FFT filter blocks. + * + * \param decimation >= 1 + * \param taps complex filter taps + * \param nthreads number of threads for the FFT to use + */ + static sptr make(int decimation, + const std::vector<float> &taps, + int nthreads=1); + + virtual void set_taps(const std::vector<float> &taps) = 0; + virtual std::vector<float> taps() const = 0; + + /*! + * \brief Set number of threads to use. + */ + virtual void set_nthreads(int n) = 0; + + /*! + * \brief Get number of threads being used. + */ + virtual int nthreads() const = 0; + }; + + } /* namespace filter */ +} /* namespace gr */ + +#endif /* INCLUDED_FILTER_FFT_FILTER_CCF_H */ diff --git a/gr-filter/include/gnuradio/filter/pfb_channelizer_ccf.h b/gr-filter/include/gnuradio/filter/pfb_channelizer_ccf.h index f199f85b83..96ffd60daa 100644 --- a/gr-filter/include/gnuradio/filter/pfb_channelizer_ccf.h +++ b/gr-filter/include/gnuradio/filter/pfb_channelizer_ccf.h @@ -101,8 +101,14 @@ namespace gr { * <B><EM>f. harris, "Multirate Signal Processing for Communication * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004.</EM></B> * + * When dealing with oversampling, the above book is still a good + * reference along with this paper: + * + * <B><EM>E. Venosa, X. Chen, and fred harris, “Polyphase analysis + * filter bank down-converts unequal channel bandwidths with + * arbitrary center frequencies - design I,” in SDR’10-WinnComm, + * 2010.</EM></B> */ - class FILTER_API pfb_channelizer_ccf : virtual public block { public: @@ -151,7 +157,7 @@ namespace gr { * Print all of the filterbank taps to screen. */ virtual void print_taps() = 0; - + /*! * Return a vector<vector<>> of the filterbank taps */ @@ -189,7 +195,7 @@ namespace gr { * the map is [0...M-1] with M = N. */ virtual void set_channel_map(const std::vector<int> &map) = 0; - + /*! * Gets the current channel map. */ diff --git a/gr-filter/include/gnuradio/filter/pfb_decimator_ccf.h b/gr-filter/include/gnuradio/filter/pfb_decimator_ccf.h index da4eb2bd34..06f329a03e 100644 --- a/gr-filter/include/gnuradio/filter/pfb_decimator_ccf.h +++ b/gr-filter/include/gnuradio/filter/pfb_decimator_ccf.h @@ -99,10 +99,19 @@ namespace gr { * \param decim (unsigned integer) Specifies the decimation rate to use * \param taps (vector/list of floats) The prototype filter to populate the filterbank. * \param channel (unsigned integer) Selects the channel to return [default=0]. + * \param use_fft_rotator (bool) Rotate channels using FFT method instead of exp(phi). + * For larger values of \p channel, the FFT method will perform better. + * Generally, this value of \p channel is small (~5), but could be + * architecture-specific (Default: true). + * \param use_fft_filters (bool) Use FFT filters (fast convolution) instead of FIR filters. + * FFT filters perform better for larger numbers of taps but is + * architecture-specific (Default: true). */ static sptr make(unsigned int decim, - const std::vector<float> &taps, - unsigned int channel); + const std::vector<float> &taps, + unsigned int channel, + bool use_fft_rotator=true, + bool use_fft_filters=true); /*! * Resets the filterbank's filter taps with the new prototype filter diff --git a/gr-filter/include/gnuradio/filter/pfb_synthesizer_ccf.h b/gr-filter/include/gnuradio/filter/pfb_synthesizer_ccf.h index f1fbad7272..32465b61ac 100644 --- a/gr-filter/include/gnuradio/filter/pfb_synthesizer_ccf.h +++ b/gr-filter/include/gnuradio/filter/pfb_synthesizer_ccf.h @@ -34,6 +34,57 @@ namespace gr { * \brief Polyphase synthesis filterbank with * gr_complex input, gr_complex output and float taps * \ingroup channelizers_blk + * + * \details + * + * The PFB sythesis filterbank combines multiple baseband signals + * into a single channelized signal. Each input stream is, + * essentially, modulated onto an output channel according the the + * channel mapping (see set_channel_map for details). + * + * Setting this filterbank up means selecting the number of output + * channels, the prototype filter, and whether to handle channels + * at 2x the sample rate (this is generally used only for + * reconstruction filtering). + * + * The number of channels sets the maximum number of channels to + * use, but not all input streams must be connected. For M total + * channels, we can connect inputs 0 to N where N < M-1. Because + * of the way GNU Radio handles stream connections, we must + * connect the channels consecutively, and so we must use the + * set_channel_map if the desired output channels are not the same + * as the the default mapping. This features gives us the + * flexibility to output to any given channel. Generally, we try + * to not use the channels at the edge of the spectrum to avoid + * issues with filtering and roll-off of the transmitter or + * receiver. + * + * When using the 2x sample rate mode, we specify the number of + * channels that will be used. However, the actual output signal + * will be twice this number of channels. This is mainly important + * to know when setting the channel map. For M channels, the + * channel mapping can specy from 0 to 2M-1 channels to output + * onto. + * + * For more details about this and the concepts of reconstruction + * filtering, see: + * + * <B><EM>f. harris, "Multirate Signal Processing for Communication + * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004.</EM></B> + * + * <B><EM>X. Chen, E. Venosa, and fred harris, “Polyphase analysis + * filter bank up-converts unequal channel bandwidths with + * arbitrary center frequencies - design ii,” in SDR’10-WinnComm, + * 2010.</EM></B> + * + * <B><EM>E. Venosa, X. Chen, and fred harris, “Polyphase analysis + * filter bank down-converts unequal channel bandwidths with + * arbitrary center frequencies - design I,” in SDR’10-WinnComm, + * 2010.</EM></B> + * + * <B><EM>f. j. harris, C. Dick, X. Chen, and E. Venosa, “Wideband 160- + * channel polyphase filter bank cable TV channeliser,” in IET + * Signal Processing, 2010.</EM></B> */ class FILTER_API pfb_synthesizer_ccf : virtual public sync_interpolator { @@ -50,8 +101,8 @@ namespace gr { * \param twox (bool) use 2x oversampling or not (default is no) */ static sptr make(unsigned int numchans, - const std::vector<float> &taps, - bool twox=false); + const std::vector<float> &taps, + bool twox=false); /*! * Resets the filterbank's filter taps with the new prototype filter diff --git a/gr-filter/include/gnuradio/filter/polyphase_filterbank.h b/gr-filter/include/gnuradio/filter/polyphase_filterbank.h index fad04e3294..1745a470ab 100644 --- a/gr-filter/include/gnuradio/filter/polyphase_filterbank.h +++ b/gr-filter/include/gnuradio/filter/polyphase_filterbank.h @@ -26,6 +26,7 @@ #include <gnuradio/filter/api.h> #include <gnuradio/filter/fir_filter.h> +#include <gnuradio/filter/fft_filter.h> #include <gnuradio/fft/fft.h> namespace gr { @@ -101,11 +102,12 @@ namespace gr { { protected: unsigned int d_nfilts; - std::vector<kernel::fir_filter_ccf*> d_filters; + std::vector<kernel::fir_filter_ccf*> d_fir_filters; + std::vector<kernel::fft_filter_ccf*> d_fft_filters; std::vector< std::vector<float> > d_taps; unsigned int d_taps_per_filter; fft::fft_complex *d_fft; - + public: /*! * Build the polyphase filterbank decimator. @@ -113,9 +115,11 @@ namespace gr { * channels <EM>M</EM> * \param taps (vector/list of floats) The prototype filter to * populate the filterbank. + * \param fft_forward (bool) use a forward or inverse FFT (default=false). */ polyphase_filterbank(unsigned int nfilts, - const std::vector<float> &taps); + const std::vector<float> &taps, + bool fft_forward=false); ~polyphase_filterbank(); @@ -126,7 +130,7 @@ namespace gr { * \param taps (vector/list of floats) The prototype filter to * populate the filterbank. */ - void set_taps(const std::vector<float> &taps); + virtual void set_taps(const std::vector<float> &taps); /*! * Print all of the filterbank taps to screen. diff --git a/gr-filter/lib/CMakeLists.txt b/gr-filter/lib/CMakeLists.txt index 9f17092bb1..9889a4ddc1 100644 --- a/gr-filter/lib/CMakeLists.txt +++ b/gr-filter/lib/CMakeLists.txt @@ -128,6 +128,7 @@ list(APPEND filter_sources dc_blocker_ff_impl.cc filter_delay_fc_impl.cc fft_filter_ccc_impl.cc + fft_filter_ccf_impl.cc fft_filter_fff_impl.cc fractional_interpolator_cc_impl.cc fractional_interpolator_ff_impl.cc diff --git a/gr-filter/lib/fft_filter.cc b/gr-filter/lib/fft_filter.cc index e9c381fa2e..fcf5f64c63 100644 --- a/gr-filter/lib/fft_filter.cc +++ b/gr-filter/lib/fft_filter.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2010,2012 Free Software Foundation, Inc. + * Copyright 2010,2012,2014 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -345,6 +345,170 @@ namespace gr { return nitems; } + + /**************************************************************/ + + + fft_filter_ccf::fft_filter_ccf(int decimation, + const std::vector<float> &taps, + int nthreads) + : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), + d_invfft(0), d_nthreads(nthreads), d_xformed_taps(NULL) + { + set_taps(taps); + } + + fft_filter_ccf::~fft_filter_ccf() + { + delete d_fwdfft; + delete d_invfft; + if(d_xformed_taps != NULL) + volk_free(d_xformed_taps); + } + + /* + * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps + */ + int + fft_filter_ccf::set_taps(const std::vector<float> &taps) + { + int i = 0; + d_taps = taps; + 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] = gr_complex(taps[i] * scale, 0.0f); + + for(; i < d_fftsize; i++) + in[i] = gr_complex(0.0f, 0.0f); + + 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 + fft_filter_ccf::compute_sizes(int ntaps) + { + int old_fftsize = d_fftsize; + d_ntaps = ntaps; + d_fftsize = (int) (2 * pow(2.0, ceil(log(double(ntaps)) / log(2.0)))); + d_nsamples = d_fftsize - d_ntaps + 1; + + if(VERBOSE) { + std::cerr << "fft_filter_ccf: ntaps = " << d_ntaps + << " fftsize = " << d_fftsize + << " nsamples = " << d_nsamples << std::endl; + } + + // compute new plans + if(d_fftsize != old_fftsize) { + delete d_fwdfft; + delete d_invfft; + if(d_xformed_taps != NULL) + volk_free(d_xformed_taps); + d_fwdfft = new fft::fft_complex(d_fftsize, true, d_nthreads); + d_invfft = new fft::fft_complex(d_fftsize, false, d_nthreads); + d_xformed_taps = (gr_complex*)volk_malloc(sizeof(gr_complex)*d_fftsize, + volk_get_alignment()); + } + } + + void + fft_filter_ccf::set_nthreads(int n) + { + d_nthreads = n; + if(d_fwdfft) + d_fwdfft->set_nthreads(n); + if(d_invfft) + d_invfft->set_nthreads(n); + } + + std::vector<float> + fft_filter_ccf::taps() const + { + return d_taps; + } + + unsigned int + fft_filter_ccf::ntaps() const + { + return d_ntaps; + } + + unsigned int + fft_filter_ccf::filtersize() const + { + return d_fftsize; + } + + int + fft_filter_ccf::nthreads() const + { + return d_nthreads; + } + + int + fft_filter_ccf::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; + gr_complex *c = d_invfft->get_inbuf(); + + volk_32fc_x2_multiply_32fc_a(c, a, b, d_fftsize); + + 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)); + } + + return nitems; + } + + } /* namespace kernel */ } /* namespace filter */ } /* namespace gr */ diff --git a/gr-filter/lib/fft_filter_ccf_impl.cc b/gr-filter/lib/fft_filter_ccf_impl.cc new file mode 100644 index 0000000000..00e1842fff --- /dev/null +++ b/gr-filter/lib/fft_filter_ccf_impl.cc @@ -0,0 +1,120 @@ +/* -*- c++ -*- */ +/* + * Copyright 2014 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 "fft_filter_ccf_impl.h" +#include <gnuradio/io_signature.h> + +#include <math.h> +#include <assert.h> +#include <stdexcept> + +namespace gr { + namespace filter { + + fft_filter_ccf::sptr + fft_filter_ccf::make(int decimation, + const std::vector<float> &taps, + int nthreads) + { + return gnuradio::get_initial_sptr + (new fft_filter_ccf_impl + (decimation, taps, nthreads)); + } + + fft_filter_ccf_impl::fft_filter_ccf_impl(int decimation, + const std::vector<float> &taps, + int nthreads) + : sync_decimator("fft_filter_ccf", + io_signature::make(1, 1, sizeof(gr_complex)), + io_signature::make(1, 1, sizeof(gr_complex)), + decimation), + d_updated(false) + { + set_history(1); + + d_filter = new kernel::fft_filter_ccf(decimation, taps, nthreads); + + d_new_taps = taps; + d_nsamples = d_filter->set_taps(taps); + set_output_multiple(d_nsamples); + } + + fft_filter_ccf_impl::~fft_filter_ccf_impl() + { + delete d_filter; + } + + void + fft_filter_ccf_impl::set_taps(const std::vector<float> &taps) + { + d_new_taps = taps; + d_updated = true; + } + + std::vector<float> + fft_filter_ccf_impl::taps() const + { + return d_new_taps; + } + + void + fft_filter_ccf_impl::set_nthreads(int n) + { + if(d_filter) + d_filter->set_nthreads(n); + } + + int + fft_filter_ccf_impl::nthreads() const + { + if(d_filter) + return d_filter->nthreads(); + else + return 0; + } + + int + fft_filter_ccf_impl::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]; + + if (d_updated){ + d_nsamples = d_filter->set_taps(d_new_taps); + d_updated = false; + set_output_multiple(d_nsamples); + return 0; // output multiple may have changed + } + + d_filter->filter(noutput_items, in, out); + + return noutput_items; + } + + } /* namespace filter */ +} /* namespace gr */ diff --git a/gr-filter/lib/fft_filter_ccf_impl.h b/gr-filter/lib/fft_filter_ccf_impl.h new file mode 100644 index 0000000000..2e0ddd46fb --- /dev/null +++ b/gr-filter/lib/fft_filter_ccf_impl.h @@ -0,0 +1,62 @@ +/* -*- c++ -*- */ +/* + * Copyright 2014 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_FILTER_FFT_FILTER_CCF_IMPL_H +#define INCLUDED_FILTER_FFT_FILTER_CCF_IMPL_H + +#include <gnuradio/filter/api.h> +#include <gnuradio/filter/fft_filter.h> +#include <gnuradio/filter/fft_filter_ccf.h> + +namespace gr { + namespace filter { + + class FILTER_API fft_filter_ccf_impl : public fft_filter_ccf + { + private: + int d_nsamples; + bool d_updated; + kernel::fft_filter_ccf *d_filter; + std::vector<float> d_new_taps; + + public: + fft_filter_ccf_impl(int decimation, + const std::vector<float> &taps, + int nthreads=1); + + ~fft_filter_ccf_impl(); + + void set_taps(const std::vector<float> &taps); + std::vector<float> taps() const; + + void set_nthreads(int n); + int nthreads() const; + + int work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + }; + + } /* namespace filter */ +} /* namespace gr */ + +#endif /* INCLUDED_FILTER_FFT_FILTER_CCF_IMPL_H */ diff --git a/gr-filter/lib/pfb_channelizer_ccf_impl.cc b/gr-filter/lib/pfb_channelizer_ccf_impl.cc index c28434b6cb..62223e40da 100644 --- a/gr-filter/lib/pfb_channelizer_ccf_impl.cc +++ b/gr-filter/lib/pfb_channelizer_ccf_impl.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009,2010,2012 Free Software Foundation, Inc. + * Copyright 2009,2010,2012,2014 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -26,25 +26,28 @@ #include "pfb_channelizer_ccf_impl.h" #include <gnuradio/io_signature.h> +#include <stdio.h> namespace gr { namespace filter { - - pfb_channelizer_ccf::sptr pfb_channelizer_ccf::make(unsigned int nfilts, - const std::vector<float> &taps, - float oversample_rate) + + pfb_channelizer_ccf::sptr + pfb_channelizer_ccf::make(unsigned int nfilts, + const std::vector<float> &taps, + float oversample_rate) { - return gnuradio::get_initial_sptr(new pfb_channelizer_ccf_impl(nfilts, taps, - oversample_rate)); + return gnuradio::get_initial_sptr + (new pfb_channelizer_ccf_impl(nfilts, taps, + oversample_rate)); } pfb_channelizer_ccf_impl::pfb_channelizer_ccf_impl(unsigned int nfilts, const std::vector<float> &taps, float oversample_rate) : block("pfb_channelizer_ccf", - io_signature::make(nfilts, nfilts, sizeof(gr_complex)), - io_signature::make(1, nfilts, sizeof(gr_complex))), - polyphase_filterbank(nfilts, taps), + io_signature::make(nfilts, nfilts, sizeof(gr_complex)), + io_signature::make(1, nfilts, sizeof(gr_complex))), + polyphase_filterbank(nfilts, taps, false), d_updated(false), d_oversample_rate(oversample_rate) { // The over sampling rate must be rationally related to the number of channels @@ -59,15 +62,22 @@ namespace gr { set_relative_rate(1.0/intp); - // Default channel map + // Default channel map. The channel map specifies which input + // goes to which output channel; so out[0] comes from + // channel_map[0]. d_channel_map.resize(d_nfilts); for(unsigned int i = 0; i < d_nfilts; i++) { d_channel_map[i] = i; } - // 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. + // We use a look up table to set the index of the FFT input + // buffer, which equivalently performs the FFT shift operation + // on every other turn when the rate_ratio>1. Also, this + // performs the index 'flip' where the first input goes into the + // last filter. In the pfb_decimator_ccf, we directly index the + // input_items buffers starting with this last; here we start + // with the first and put it into the fft object properly for + // the same effect. d_rate_ratio = (int)rintf(d_nfilts / d_oversample_rate); d_idxlut = new int[d_nfilts]; for(unsigned int i = 0; i < d_nfilts; i++) { @@ -81,10 +91,8 @@ namespace gr { d_output_multiple++; set_output_multiple(d_output_multiple); - // History is the length of each filter arm plus 1. - // The +1 comes from the channel mapping in the work function - // where we start n=1 so that we can look at in[n-1] - set_history(d_taps_per_filter+1); + // Use set_taps to also set the history requirement + set_taps(taps); } pfb_channelizer_ccf_impl::~pfb_channelizer_ccf_impl() @@ -118,7 +126,7 @@ namespace gr { pfb_channelizer_ccf_impl::set_channel_map(const std::vector<int> &map) { gr::thread::scoped_lock guard(d_mutex); - + if(map.size() > 0) { unsigned int max = (unsigned int)*std::max_element(map.begin(), map.end()); if(max >= d_nfilts) { @@ -142,9 +150,9 @@ namespace gr { { gr::thread::scoped_lock guard(d_mutex); - gr_complex *in = (gr_complex *) input_items[0]; - gr_complex *out = (gr_complex *) output_items[0]; - + gr_complex *in = (gr_complex*)input_items[0]; + gr_complex *out = (gr_complex*)output_items[0]; + if(d_updated) { d_updated = false; return 0; // history requirements may have changed. @@ -152,6 +160,18 @@ namespace gr { size_t noutputs = output_items.size(); + // The following algorithm looks more complex in order to handle + // the cases where we want more that 1 sps for each + // channel. Otherwise, this would boil down into a single loop + // that operates from input_items[0] to [d_nfilts]. + + // When dealing with osps>1, we start not at the last filter, + // but nfilts/osps and then wrap around to the next symbol into + // the other set of filters. + // For details of this operation, see: + // fred harris, Multirate Signal Processing For Communication + // Systems. Upper Saddle River, NJ: Prentice Hall, 2004. + int n=1, i=-1, j=0, oo=0, last; int toconsume = (int)rintf(noutput_items/d_oversample_rate); while(n <= toconsume) { @@ -159,16 +179,16 @@ namespace gr { i = (i + d_rate_ratio) % d_nfilts; last = i; while(i >= 0) { - in = (gr_complex*)input_items[j]; - d_fft->get_inbuf()[d_idxlut[j]] = d_filters[i]->filter(&in[n]); + in = (gr_complex*)input_items[j]; + d_fft->get_inbuf()[d_idxlut[j]] = d_fir_filters[i]->filter(&in[n]); j++; i--; } i = d_nfilts-1; while(i > last) { - in = (gr_complex*)input_items[j]; - d_fft->get_inbuf()[d_idxlut[j]] = d_filters[i]->filter(&in[n-1]); + in = (gr_complex*)input_items[j]; + d_fft->get_inbuf()[d_idxlut[j]] = d_fir_filters[i]->filter(&in[n-1]); j++; i--; } diff --git a/gr-filter/lib/pfb_channelizer_ccf_impl.h b/gr-filter/lib/pfb_channelizer_ccf_impl.h index 6d727676b2..bfc53f8848 100644 --- a/gr-filter/lib/pfb_channelizer_ccf_impl.h +++ b/gr-filter/lib/pfb_channelizer_ccf_impl.h @@ -31,7 +31,7 @@ namespace gr { namespace filter { - + class FILTER_API pfb_channelizer_ccf_impl : public pfb_channelizer_ccf, kernel::polyphase_filterbank { private: diff --git a/gr-filter/lib/pfb_decimator_ccf_impl.cc b/gr-filter/lib/pfb_decimator_ccf_impl.cc index 19cf20f410..7b4bf73b39 100644 --- a/gr-filter/lib/pfb_decimator_ccf_impl.cc +++ b/gr-filter/lib/pfb_decimator_ccf_impl.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2009,2010,2012 Free Software Foundation, Inc. + * Copyright 2009,2010,2012,2014 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -26,34 +26,53 @@ #include "pfb_decimator_ccf_impl.h" #include <gnuradio/io_signature.h> +#include <gnuradio/expj.h> +#include <volk/volk.h> namespace gr { namespace filter { - + pfb_decimator_ccf::sptr pfb_decimator_ccf::make(unsigned int decim, const std::vector<float> &taps, - unsigned int channel) + unsigned int channel, + bool use_fft_rotator, + bool use_fft_filters) { return gnuradio::get_initial_sptr - (new pfb_decimator_ccf_impl(decim, taps, channel)); + (new pfb_decimator_ccf_impl(decim, taps, channel, + use_fft_rotator, + use_fft_filters)); } - pfb_decimator_ccf_impl::pfb_decimator_ccf_impl(unsigned int decim, const std::vector<float> &taps, - unsigned int channel) + unsigned int channel, + bool use_fft_rotator, + bool use_fft_filters) : sync_block("pfb_decimator_ccf", - io_signature::make(decim, decim, sizeof(gr_complex)), - io_signature::make(1, 1, sizeof(gr_complex))), + io_signature::make(decim, decim, sizeof(gr_complex)), + io_signature::make(1, 1, sizeof(gr_complex))), polyphase_filterbank(decim, taps), - d_updated(false), d_chan(channel) + d_updated(false), d_chan(channel), + d_use_fft_rotator(use_fft_rotator), + d_use_fft_filters(use_fft_filters) { d_rate = decim; d_rotator = new gr_complex[d_rate]; + for(unsigned int i = 0; i < d_rate; i++) { + d_rotator[i] = gr_expj(i*d_chan*2*M_PI/d_rate); + } set_relative_rate(1.0/(float)decim); - set_history(d_taps_per_filter); + + if(d_use_fft_filters) { + set_history(1); + set_output_multiple(d_fft_filters[0]->filtersize() - d_fft_filters[0]->ntaps() + 1); + } + else { + set_history(d_taps_per_filter); + } } pfb_decimator_ccf_impl::~pfb_decimator_ccf_impl() @@ -90,8 +109,6 @@ namespace gr { d_chan = chan; } -#define ROTATEFFT - int pfb_decimator_ccf_impl::work(int noutput_items, gr_vector_const_void_star &input_items, @@ -99,47 +116,145 @@ namespace gr { { gr::thread::scoped_lock guard(d_mutex); - gr_complex *in; - gr_complex *out = (gr_complex *)output_items[0]; - if(d_updated) { d_updated = false; return 0; // history requirements may have changed. } + if(d_use_fft_rotator) { + if(d_use_fft_filters) { + return work_fft_fft(noutput_items, input_items, output_items); + } + else { + return work_fir_fft(noutput_items, input_items, output_items); + } + } + else { + if(d_use_fft_filters) { + return work_fft_exp(noutput_items, input_items, output_items); + } + else { + return work_fir_exp(noutput_items, input_items, output_items); + } + } + } + + int + pfb_decimator_ccf_impl::work_fir_exp(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) + { + gr_complex *in; + gr_complex *out = (gr_complex *)output_items[0]; + int i; for(i = 0; i < noutput_items; i++) { // Move through filters from bottom to top out[i] = 0; for(int j = d_rate-1; j >= 0; j--) { - // Take in the items from the first input stream to d_rate + // Take items from M-1 to 0; filter and rotate in = (gr_complex*)input_items[d_rate - 1 - j]; + out[i] += d_fir_filters[j]->filter(&in[i])*d_rotator[j]; + } + } - // Filter current input stream from bottom filter to top - // The rotate them by expj(j*k*2pi/M) where M is the number of filters - // (the decimation rate) and k is the channel number to extract - - // This is the real math that goes on; we abuse the FFT to do this quickly - // for decimation rates > N where N is a small number (~5): - // out[i] += d_filters[j]->filter(&in[i])*gr_expj(j*d_chan*2*M_PI/d_rate); -#ifdef ROTATEFFT - d_fft->get_inbuf()[j] = d_filters[j]->filter(&in[i]); -#else - out[i] += d_filters[j]->filter(&in[i])*d_rotator[i]; -#endif + return noutput_items; + } + + int + pfb_decimator_ccf_impl::work_fir_fft(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) + { + gr_complex *in; + gr_complex *out = (gr_complex *)output_items[0]; + + int i; + for(i = 0; i < noutput_items; i++) { + // Move through filters from bottom to top + out[i] = 0; + for(unsigned int j = 0; j < d_rate; j++) { + // Take in the items from the first input stream to d_rate + in = (gr_complex*)input_items[d_rate-1-j]; + d_fft->get_inbuf()[j] = d_fir_filters[j]->filter(&in[i]); } -#ifdef ROTATEFFT - // Perform the FFT to do the complex multiply despinning for all channels - d_fft->execute(); + // Perform the FFT to do the complex multiply despinning for all channels + d_fft->execute(); - // Select only the desired channel out - out[i] = d_fft->get_outbuf()[d_chan]; -#endif + // Select only the desired channel out + out[i] = d_fft->get_outbuf()[d_chan]; + } + + return noutput_items; + } + + int + pfb_decimator_ccf_impl::work_fft_exp(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) + { + gr_complex *in; + gr_complex *out = (gr_complex *)output_items[0]; + + int i; + gr_complex *tmp = fft::malloc_complex(noutput_items*d_rate); + + // Filter each input stream by the FFT filters; do all + // noutput_items at once to avoid repeated calls to the FFT + // setup and operation. + for(unsigned int j = 0; j < d_rate; j++) { + in = (gr_complex*)input_items[d_rate-j-1]; + d_fft_filters[j]->filter(noutput_items, in, &(tmp[j*noutput_items])); + } + + // Rotate and add filter outputs (k=channel number; M=number of + // channels; and x[j] is the output of filter j. + // y[i] = \sum_{j=0}{M-1} (x[j][i]*exp(2j \pi j k / M)) + for(i = 0; i < noutput_items; i++) { + out[i] = 0; + for(unsigned int j = 0; j < d_rate; j++) { + out[i] += tmp[j*noutput_items+i]*d_rotator[j]; + } } return noutput_items; } + int + pfb_decimator_ccf_impl::work_fft_fft(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) + { + gr_complex *in; + gr_complex *out = (gr_complex *)output_items[0]; + + int i; + gr_complex *tmp = fft::malloc_complex(noutput_items*d_rate); + + for(unsigned int j = 0; j < d_rate; j++) { + in = (gr_complex*)input_items[d_rate-j-1]; + d_fft_filters[j]->filter(noutput_items, in, &tmp[j*noutput_items]); + } + + // Performs the rotate and add operations by implementing it as + // an FFT. + for(i = 0; i < noutput_items; i++) { + for(unsigned int j = 0; j < d_rate; j++) { + d_fft->get_inbuf()[j] = tmp[j*noutput_items + i]; + } + + // Perform the FFT to do the complex multiply despinning for all channels + d_fft->execute(); + + // Select only the desired channel out + out[i] = d_fft->get_outbuf()[d_chan]; + } + + fft::free(tmp); + return noutput_items; + } + + } /* namespace filter */ } /* namespace gr */ diff --git a/gr-filter/lib/pfb_decimator_ccf_impl.h b/gr-filter/lib/pfb_decimator_ccf_impl.h index 2df8a506f0..3397701cf9 100644 --- a/gr-filter/lib/pfb_decimator_ccf_impl.h +++ b/gr-filter/lib/pfb_decimator_ccf_impl.h @@ -26,26 +26,42 @@ #include <gnuradio/filter/pfb_decimator_ccf.h> #include <gnuradio/filter/polyphase_filterbank.h> -#include <gnuradio/filter/fir_filter.h> -#include <gnuradio/fft/fft.h> #include <gnuradio/thread/thread.h> namespace gr { namespace filter { - + class FILTER_API pfb_decimator_ccf_impl : public pfb_decimator_ccf, kernel::polyphase_filterbank { private: bool d_updated; unsigned int d_rate; unsigned int d_chan; + bool d_use_fft_rotator; + bool d_use_fft_filters; gr_complex *d_rotator; gr::thread::mutex d_mutex; // mutex to protect set/work access - + + inline int work_fir_exp(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + inline int work_fir_fft(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + inline int work_fft_exp(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + inline int work_fft_fft(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + + public: pfb_decimator_ccf_impl(unsigned int decim, const std::vector<float> &taps, - unsigned int channel); + unsigned int channel, + bool use_fft_rotator=true, + bool use_fft_filters=true); ~pfb_decimator_ccf_impl(); diff --git a/gr-filter/lib/pfb_interpolator_ccf_impl.cc b/gr-filter/lib/pfb_interpolator_ccf_impl.cc index a0ed73fd4e..cf79c103d3 100644 --- a/gr-filter/lib/pfb_interpolator_ccf_impl.cc +++ b/gr-filter/lib/pfb_interpolator_ccf_impl.cc @@ -29,7 +29,7 @@ namespace gr { namespace filter { - + pfb_interpolator_ccf::sptr pfb_interpolator_ccf::make(unsigned int interp, const std::vector<float> &taps) @@ -94,7 +94,7 @@ namespace gr { while(i < noutput_items) { for(unsigned int j = 0; j < d_rate; j++) { - out[i] = d_filters[j]->filter(&in[count]); + out[i] = d_fir_filters[j]->filter(&in[count]); i++; } count++; diff --git a/gr-filter/lib/pfb_synthesizer_ccf_impl.cc b/gr-filter/lib/pfb_synthesizer_ccf_impl.cc index 5f0e330ab4..bdfc158a5d 100644 --- a/gr-filter/lib/pfb_synthesizer_ccf_impl.cc +++ b/gr-filter/lib/pfb_synthesizer_ccf_impl.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2010,2012 Free Software Foundation, Inc. + * Copyright 2010,2012,2014 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -33,45 +33,43 @@ namespace gr { pfb_synthesizer_ccf::sptr pfb_synthesizer_ccf::make(unsigned int numchans, - const std::vector<float> &taps, - bool twox) + const std::vector<float> &taps, bool twox) { return gnuradio::get_initial_sptr - (new pfb_synthesizer_ccf_impl(numchans, taps, twox)); + (new pfb_synthesizer_ccf_impl(numchans, taps, twox)); } - - pfb_synthesizer_ccf_impl::pfb_synthesizer_ccf_impl(unsigned int numchans, - const std::vector<float> &taps, - bool twox) + pfb_synthesizer_ccf_impl::pfb_synthesizer_ccf_impl + (unsigned int numchans, const std::vector<float> &taps, bool twox) : sync_interpolator("pfb_synthesizer_ccf", - io_signature::make(1, numchans, sizeof(gr_complex)), - io_signature::make(1, 1, sizeof(gr_complex)), - (twox ? numchans/2 : numchans)), - d_updated(false), d_numchans(numchans), d_state(0) + io_signature::make(1, numchans, sizeof(gr_complex)), + io_signature::make(1, 1, sizeof(gr_complex)), + numchans), + d_updated(false), d_numchans(numchans), d_state(0) { // set up 2x multiplier; if twox==True, set to 2, otherwise to 1 d_twox = (twox ? 2 : 1); if(d_numchans % d_twox != 0) { - throw std::invalid_argument("pfb_synthesizer_ccf: number of channels must be even for 2x oversampling.\n"); + throw std::invalid_argument("pfb_synthesizer_ccf_impl: number of channels must be even for 2x oversampling.\n"); } - d_filters = std::vector<kernel::fir_filter_with_buffer_ccf*>(d_numchans); - d_channel_map.resize(d_numchans); + d_filters = std::vector<kernel::fir_filter_with_buffer_ccf*>(d_twox*d_numchans); + d_channel_map.resize(d_twox*d_numchans); - // Create a FIR filter for each channel and zero out the taps - std::vector<float> vtaps(0, d_numchans); - for(unsigned int i = 0; i < d_numchans; i++) { - d_filters[i] = new kernel::fir_filter_with_buffer_ccf(vtaps); - d_channel_map[i] = i; + // Create an FIR filter for each channel and zero out the taps + // and set the default channel map + std::vector<float> vtaps(0, d_twox*d_numchans); + for(unsigned int i = 0; i < d_twox*d_numchans; i++) { + d_filters[i] = new kernel::fir_filter_with_buffer_ccf(vtaps); + d_channel_map[i] = i; } // Now, actually set the filters' taps set_taps(taps); // Create the IFFT to handle the input channel rotations - d_fft = new fft::fft_complex(d_numchans, false); - memset(d_fft->get_inbuf(), 0, d_numchans*sizeof(gr_complex)); + d_fft = new fft::fft_complex(d_twox*d_numchans, false); + memset(d_fft->get_inbuf(), 0, d_twox*d_numchans*sizeof(gr_complex)); set_output_multiple(d_numchans); } @@ -79,8 +77,8 @@ namespace gr { pfb_synthesizer_ccf_impl::~pfb_synthesizer_ccf_impl() { delete d_fft; - for(unsigned int i = 0; i < d_numchans; i++) { - delete d_filters[i]; + for(unsigned int i = 0; i < d_twox*d_numchans; i++) { + delete d_filters[i]; } } @@ -89,20 +87,29 @@ namespace gr { { gr::thread::scoped_lock guard(d_mutex); + // The different modes, 1x or 2x the sampling rate, have + // different filtering partitions. if(d_twox == 1) - set_taps1(taps); + set_taps1(taps); else - set_taps2(taps); - - // Set the history to ensure enough input items for each filter - set_history(d_taps_per_filter+1); + set_taps2(taps); + // Because we keep our own buffers inside the filters, we don't + // need history. + set_history(1); d_updated = true; } void pfb_synthesizer_ccf_impl::set_taps1(const std::vector<float> &taps) { + // In this partitioning, we do a normal polyphase paritioning by + // deinterleaving the taps into each filter: + // + // Prototype filter: [t0, t1, t2, t3, t4, t5, t6] + // filter 0: [t0, t3, t6] + // filter 1: [t1, t4, 0 ] + // filter 2: [t2, t5, 0 ] unsigned int i,j; unsigned int ntaps = taps.size(); @@ -113,70 +120,77 @@ namespace gr { // 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; + std::vector<float> tmp_taps = taps; while((float)(tmp_taps.size()) < d_numchans*d_taps_per_filter) { - tmp_taps.push_back(0.0); + tmp_taps.push_back(0.0); } // Partition the filter for(i = 0; i < d_numchans; i++) { - // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out - d_taps[i] = std::vector<float>(d_taps_per_filter, 0); - for(j = 0; j < d_taps_per_filter; j++) { - d_taps[i][j] = tmp_taps[i + j*d_numchans]; // add taps to channels in reverse order - } - - // Build a filter for each channel and add it's taps to it - d_filters[i]->set_taps(d_taps[i]); + // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out + d_taps[i] = std::vector<float>(d_taps_per_filter, 0); + for(j = 0; j < d_taps_per_filter; j++) { + d_taps[i][j] = tmp_taps[i + j*d_numchans]; // add taps to channels in reverse order + } + + // Set the filter taps for each channel + d_filters[i]->set_taps(d_taps[i]); } } void - pfb_synthesizer_ccf_impl::set_taps2 (const std::vector<float> &taps) + pfb_synthesizer_ccf_impl::set_taps2(const std::vector<float> &taps) { + // In this partitioning, create 2M filters each in Z^2; the + // second half of the filters are also delayed by Z^1: + // + // Prototype filter: [t0, t1, t2, t3, t4, t5, t6] + // filter 0: [t0, 0, t6] + // filter 1: [t2, 0, 0] + // filter 2: [t4, 0, 0] + // filter 3: [ 0, t1, 0] + // filter 4: [ 0, t3, 0 ] + // filter 5: [ 0, t5, 0 ] + unsigned int i,j; int state = 0; - unsigned int ntaps = 2*taps.size(); + unsigned int ntaps = taps.size(); d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_numchans); // Create d_numchan vectors to store each channel's taps - d_taps.resize(d_numchans); + d_taps.resize(d_twox*d_numchans); // Make a vector of the taps plus fill it out with 0's to fill // each polyphase filter with exactly d_taps_per_filter - std::vector<float> tmp_taps; - tmp_taps = taps; + std::vector<float> tmp_taps = taps; while((float)(tmp_taps.size()) < d_numchans*d_taps_per_filter) { - tmp_taps.push_back(0.0); + tmp_taps.push_back(0.0); + //tmp_taps.insert(tmp_taps.begin(), 0.0); } // Partition the filter - unsigned int halfchans = d_numchans/2; - for(i = 0; i < halfchans; 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); - d_taps[halfchans+i] = std::vector<float>(d_taps_per_filter, 0); - state = 0; - for(j = 0; j < d_taps_per_filter; j++) { - // add taps to channels in reverse order - // Zero out every other tap - if(state == 0) { - d_taps[i][j] = tmp_taps[i + j*halfchans]; - d_taps[halfchans + i][j] = 0; - state = 1; - } - else { - d_taps[i][j] = 0; - d_taps[halfchans + i][j] = tmp_taps[i + j*halfchans]; - state = 0; - } - } - - // Build a filter for each channel and add it's taps to it - d_filters[i]->set_taps(d_taps[i]); - d_filters[halfchans + i]->set_taps(d_taps[halfchans + i]); + for(i = 0; i < d_numchans; i++) { + // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out + d_taps[i] = std::vector<float>(d_taps_per_filter, 0); + d_taps[d_numchans+i] = std::vector<float>(d_taps_per_filter, 0); + state = 0; + for(j = 0; j < d_taps_per_filter; j++) { + if(state == 0) { + d_taps[i][j] = 0; + d_taps[d_numchans + i][j] = tmp_taps[i + j*d_numchans]; + state = 1; + } + else { + d_taps[i][j] = tmp_taps[i + j*d_numchans]; + d_taps[d_numchans + i][j] = 0; + state = 0; + } + } + + // Set the filter taps for each channel + d_filters[i]->set_taps(d_taps[i]); + d_filters[d_numchans + i]->set_taps(d_taps[d_numchans + i]); } } @@ -184,15 +198,16 @@ namespace gr { pfb_synthesizer_ccf_impl::print_taps() { unsigned int i, j; - for(i = 0; i < d_numchans; i++) { - printf("filter[%d]: [", i); - for(j = 0; j < d_taps_per_filter; j++) { - printf(" %.4e", d_taps[i][j]); - } - printf("]\n\n"); + for(i = 0; i < d_twox*d_numchans; i++) { + printf("filter[%d]: [", i); + for(j = 0; j < d_taps_per_filter; j++) { + printf(" %.4e", d_taps[i][j]); + } + printf("]\n\n"); } } + std::vector< std::vector<float> > pfb_synthesizer_ccf_impl::taps() const { @@ -205,14 +220,15 @@ namespace gr { gr::thread::scoped_lock guard(d_mutex); if(map.size() > 0) { - unsigned int max = (unsigned int)*std::max_element(map.begin(), map.end()); - if(max >= d_numchans) { - throw std::invalid_argument("gr_pfb_synthesizer_ccf::set_channel_map: map range out of bounds.\n"); - } - d_channel_map = map; - - // Zero out fft buffer so that unused channels are always 0 - memset(d_fft->get_inbuf(), 0, d_numchans*sizeof(gr_complex)); + unsigned int max = (unsigned int)*std::max_element(map.begin(), map.end()); + unsigned int min = (unsigned int)*std::min_element(map.begin(), map.end()); + if((max >= d_twox*d_numchans) || (min < 0)) { + throw std::invalid_argument("pfb_synthesizer_ccf_impl::set_channel_map: map range out of bounds.\n"); + } + d_channel_map = map; + + // Zero out fft buffer so that unused channels are always 0 + memset(d_fft->get_inbuf(), 0,d_twox*d_numchans*sizeof(gr_complex)); } } @@ -224,63 +240,63 @@ namespace gr { int pfb_synthesizer_ccf_impl::work(int noutput_items, - gr_vector_const_void_star &input_items, - gr_vector_void_star &output_items) + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) { gr::thread::scoped_lock guard(d_mutex); - gr_complex *in = (gr_complex*)input_items[0]; - gr_complex *out = (gr_complex*)output_items[0]; + gr_complex *in = (gr_complex*) input_items[0]; + gr_complex *out = (gr_complex *) output_items[0]; if(d_updated) { - d_updated = false; - return 0; // history requirements may have changed. + d_updated = false; + return 0; // history requirements may have changed. } unsigned int n, i; size_t ninputs = input_items.size(); - // Algoritm for critically sampled channels + // Algorithm for critically sampled channels if(d_twox == 1) { - for(n = 0; n < noutput_items/d_numchans; n++) { - for(i = 0; i < ninputs; i++) { - in = (gr_complex*)input_items[i]; - d_fft->get_inbuf()[d_channel_map[i]] = in[n]; - } - - // spin through IFFT - d_fft->execute(); - - for(i = 0; i < d_numchans; i++) { - out[i] = d_filters[i]->filter(d_fft->get_outbuf()[i]); - } - out += d_numchans; - } + for(n = 0; n < noutput_items/d_numchans; n++) { + for(i = 0; i < ninputs; i++) { + in = (gr_complex*)input_items[i]; + d_fft->get_inbuf()[d_channel_map[i]] = in[n]; + } + + // spin through IFFT + d_fft->execute(); + + for(i = 0; i < d_numchans; i++) { + out[i] = d_filters[i]->filter(d_fft->get_outbuf()[i]); + } + out += d_numchans; + } } // Algorithm for oversampling by 2x else { - unsigned int halfchans = d_numchans/2; - for(n = 0; n < noutput_items/halfchans; n++) { - for(i = 0; i < ninputs; i++) { - in = (gr_complex*)input_items[i]; - d_fft->get_inbuf()[d_channel_map[i]] = in[n]; - } - - // spin through IFFT - d_fft->execute(); - - // Output is sum of two filters, but the input buffer to the filters must be circularly - // shifted by numchans every time through, done by using d_state to determine which IFFT - // buffer position to pull from. - for(i = 0; i < halfchans; i++) { - out[i] = d_filters[i]->filter(d_fft->get_outbuf()[d_state*halfchans+i]); - out[i] += d_filters[halfchans+i]->filter(d_fft->get_outbuf()[(d_state^1)*halfchans+i]); - } - d_state ^= 1; - - out += halfchans; - } + for(n = 0; n < noutput_items/d_numchans; n++) { + for(i = 0; i < ninputs; i++) { + //in = (gr_complex*)input_items[ninputs-i-1]; + in = (gr_complex*)input_items[i]; + d_fft->get_inbuf()[d_channel_map[i]] = in[n]; + } + + // spin through IFFT + d_fft->execute(); + + // Output is sum of two filters, but the input buffer to the filters must be circularly + // shifted by numchans every time through, done by using d_state to determine which IFFT + // buffer position to pull from. + for(i = 0; i < d_numchans; i++) { + out[i] = d_filters[i]->filter(d_fft->get_outbuf()[d_state*d_numchans+i]); + out[i] += d_filters[d_numchans+i]->filter(d_fft->get_outbuf()[(d_state^1)*d_numchans+i]); + } + d_state ^= 1; + + out += d_numchans; + } } return noutput_items; diff --git a/gr-filter/lib/polyphase_filterbank.cc b/gr-filter/lib/polyphase_filterbank.cc index ab978e8aaf..5fffc02dce 100644 --- a/gr-filter/lib/polyphase_filterbank.cc +++ b/gr-filter/lib/polyphase_filterbank.cc @@ -1,6 +1,6 @@ /* -*- c++ -*- */ /* - * Copyright 2012 Free Software Foundation, Inc. + * Copyright 2012,2014 Free Software Foundation, Inc. * * This file is part of GNU Radio * @@ -32,29 +32,33 @@ namespace gr { namespace kernel { polyphase_filterbank::polyphase_filterbank(unsigned int nfilts, - const std::vector<float> &taps) + const std::vector<float> &taps, + bool fft_forward) : d_nfilts(nfilts) { - d_filters = std::vector<kernel::fir_filter_ccf*>(d_nfilts); + d_fir_filters = std::vector<kernel::fir_filter_ccf*>(d_nfilts); + d_fft_filters = std::vector<kernel::fft_filter_ccf*>(d_nfilts); // Create an FIR filter for each channel and zero out the taps - std::vector<float> vtaps(0, d_nfilts); + std::vector<float> vtaps(1, 0.0f); for(unsigned int i = 0; i < d_nfilts; i++) { - d_filters[i] = new kernel::fir_filter_ccf(1, vtaps); + d_fir_filters[i] = new kernel::fir_filter_ccf(1, vtaps); + d_fft_filters[i] = new kernel::fft_filter_ccf(1, vtaps); } // Now, actually set the filters' taps set_taps(taps); // Create the FFT to handle the output de-spinning of the channels - d_fft = new fft::fft_complex(d_nfilts, false); + d_fft = new fft::fft_complex(d_nfilts, fft_forward); } polyphase_filterbank::~polyphase_filterbank() { delete d_fft; for(unsigned int i = 0; i < d_nfilts; i++) { - delete d_filters[i]; + delete d_fir_filters[i]; + delete d_fft_filters[i]; } } @@ -62,7 +66,6 @@ namespace gr { polyphase_filterbank::set_taps(const std::vector<float> &taps) { unsigned int i,j; - unsigned int ntaps = taps.size(); d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_nfilts); @@ -71,16 +74,11 @@ namespace gr { // 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; + std::vector<float> tmp_taps = taps; while((float)(tmp_taps.size()) < d_nfilts*d_taps_per_filter) { - tmp_taps.push_back(0.0); + tmp_taps.push_back(0.0); } - // Reverse taps here; set_taps in filter will then re-reverse, - // but order them propely, anyways. - std::reverse(tmp_taps.begin(), tmp_taps.end()); - // Partition the filter for(i = 0; i < d_nfilts; i++) { // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out @@ -89,8 +87,9 @@ namespace gr { d_taps[i][j] = tmp_taps[i + j*d_nfilts]; } - // Build a filter for each channel and add it's taps to it - d_filters[i]->set_taps(d_taps[i]); + // Set the filter taps for each channel + d_fir_filters[i]->set_taps(d_taps[i]); + d_fft_filters[i]->set_taps(d_taps[i]); } } diff --git a/gr-filter/python/filter/pfb.py b/gr-filter/python/filter/pfb.py index 7f7375cb7f..732812ef2a 100644 --- a/gr-filter/python/filter/pfb.py +++ b/gr-filter/python/filter/pfb.py @@ -77,7 +77,11 @@ class channelizer_ccf(gr.hier_block2): def set_channel_map(self, newmap): self.pfb.set_channel_map(newmap) + def set_taps(self, taps): + self.pfb.set_taps(taps) + def taps(self): + return self.pfb.taps() class interpolator_ccf(gr.hier_block2): ''' @@ -122,6 +126,8 @@ class interpolator_ccf(gr.hier_block2): self.connect(self, self.pfb) self.connect(self.pfb, self) + def set_taps(self, taps): + self.pfb.set_taps(taps) class decimator_ccf(gr.hier_block2): ''' @@ -130,10 +136,11 @@ class decimator_ccf(gr.hier_block2): This simplifies the interface by allowing a single input stream to connect to this block. It will then output a stream that is the decimated output stream. ''' - def __init__(self, decim, taps=None, channel=0, atten=100): - gr.hier_block2.__init__(self, "pfb_decimator_ccf", - gr.io_signature(1, 1, gr.sizeof_gr_complex), - gr.io_signature(1, 1, gr.sizeof_gr_complex)) + def __init__(self, decim, taps=None, channel=0, atten=100, + use_fft_rotators=True, use_fft_filters=True): + gr.hier_block2.__init__(self, "pfb_decimator_ccf", + gr.io_signature(1, 1, gr.sizeof_gr_complex), + gr.io_signature(1, 1, gr.sizeof_gr_complex)) self._decim = decim self._channel = channel @@ -160,7 +167,8 @@ class decimator_ccf(gr.hier_block2): raise RuntimeError("optfir could not generate an appropriate filter.") self.s2ss = blocks.stream_to_streams(gr.sizeof_gr_complex, self._decim) - self.pfb = filter.pfb_decimator_ccf(self._decim, self._taps, self._channel) + self.pfb = filter.pfb_decimator_ccf(self._decim, self._taps, self._channel, + use_fft_rotators, use_fft_filters) self.connect(self, self.s2ss) diff --git a/gr-filter/python/filter/qa_fft_filter.py b/gr-filter/python/filter/qa_fft_filter.py index fafd2330ae..172b945441 100755 --- a/gr-filter/python/filter/qa_fft_filter.py +++ b/gr-filter/python/filter/qa_fft_filter.py @@ -65,6 +65,18 @@ def reference_filter_fff(dec, taps, input): tb.run() return dst.data() +def reference_filter_ccf(dec, taps, input): + """ + compute result using conventional fir filter + """ + tb = gr.top_block() + src = blocks.vector_source_c(input) + op = filter.fir_filter_ccf(dec, taps) + dst = blocks.vector_sink_c() + tb.connect(src, op, dst) + tb.run() + return dst.data() + def print_complex(x): for i in x: @@ -207,6 +219,125 @@ class test_fft_filter(gr_unittest.TestCase): self.assert_fft_ok2(expected_result, result_data) # ---------------------------------------------------------------- + # test _ccf version + # ---------------------------------------------------------------- + + def test_ccf_001(self): + tb = gr.top_block() + src_data = (0,1,2,3,4,5,6,7) + taps = (1,) + expected_result = tuple([complex(x) for x in (0,1,2,3,4,5,6,7)]) + src = blocks.vector_source_c(src_data) + op = filter.fft_filter_ccf(1, taps) + dst = blocks.vector_sink_c() + tb.connect(src, op, dst) + tb.run() + result_data = dst.data() + #print 'expected:', expected_result + #print 'results: ', result_data + self.assertComplexTuplesAlmostEqual (expected_result, result_data, 5) + + + def test_ccf_002(self): + # Test nthreads + tb = gr.top_block() + src_data = (0,1,2,3,4,5,6,7) + taps = (2,) + nthreads = 2 + expected_result = tuple([2 * complex(x) for x in (0,1,2,3,4,5,6,7)]) + src = blocks.vector_source_c(src_data) + op = filter.fft_filter_ccf(1, taps, nthreads) + dst = blocks.vector_sink_c() + tb.connect(src, op, dst) + tb.run() + result_data = dst.data() + #print 'expected:', expected_result + #print 'results: ', result_data + self.assertComplexTuplesAlmostEqual (expected_result, result_data, 5) + + def test_ccf_003(self): + tb = gr.top_block() + src_data = (0,1,2,3,4,5,6,7) + taps = (2,) + expected_result = tuple([2 * complex(x) for x in (0,1,2,3,4,5,6,7)]) + src = blocks.vector_source_c(src_data) + op = filter.fft_filter_ccf(1, taps) + dst = blocks.vector_sink_c() + tb.connect(src, op, dst) + tb.run() + result_data = dst.data() + #print 'expected:', expected_result + #print 'results: ', result_data + self.assertComplexTuplesAlmostEqual (expected_result, result_data, 5) + + + def test_ccf_004(self): + random.seed(0) + for i in xrange(25): + # sys.stderr.write("\n>>> Loop = %d\n" % (i,)) + src_len = 4*1024 + src_data = make_random_complex_tuple(src_len) + ntaps = int(random.uniform(2, 1000)) + taps = make_random_float_tuple(ntaps) + expected_result = reference_filter_ccf(1, taps, src_data) + + src = blocks.vector_source_c(src_data) + op = filter.fft_filter_ccf(1, taps) + dst = blocks.vector_sink_c() + tb = gr.top_block() + tb.connect(src, op, dst) + tb.run() + result_data = dst.data() + del tb + self.assert_fft_ok2(expected_result, result_data) + + def test_ccf_005(self): + random.seed(0) + for i in xrange(25): + # sys.stderr.write("\n>>> Loop = %d\n" % (i,)) + dec = i + 1 + src_len = 4*1024 + src_data = make_random_complex_tuple(src_len) + ntaps = int(random.uniform(2, 100)) + taps = make_random_float_tuple(ntaps) + expected_result = reference_filter_ccf(dec, taps, src_data) + + src = blocks.vector_source_c(src_data) + op = filter.fft_filter_ccf(dec, taps) + dst = blocks.vector_sink_c() + tb = gr.top_block() + tb.connect(src, op, dst) + tb.run() + del tb + result_data = dst.data() + + self.assert_fft_ok2(expected_result, result_data) + + def test_ccf_006(self): + # Test decimating with nthreads=2 + random.seed(0) + nthreads = 2 + for i in xrange(25): + # sys.stderr.write("\n>>> Loop = %d\n" % (i,)) + dec = i + 1 + src_len = 4*1024 + src_data = make_random_complex_tuple(src_len) + ntaps = int(random.uniform(2, 100)) + taps = make_random_float_tuple(ntaps) + expected_result = reference_filter_ccf(dec, taps, src_data) + + src = blocks.vector_source_c(src_data) + op = filter.fft_filter_ccc(dec, taps, nthreads) + dst = blocks.vector_sink_c() + tb = gr.top_block() + tb.connect(src, op, dst) + tb.run() + del tb + result_data = dst.data() + + self.assert_fft_ok2(expected_result, result_data) + + # ---------------------------------------------------------------- # test _fff version # ---------------------------------------------------------------- diff --git a/gr-filter/python/filter/qa_pfb_channelizer.py b/gr-filter/python/filter/qa_pfb_channelizer.py index 4c108e8443..46c6e7b5ae 100755 --- a/gr-filter/python/filter/qa_pfb_channelizer.py +++ b/gr-filter/python/filter/qa_pfb_channelizer.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright 2012,2013 Free Software Foundation, Inc. +# Copyright 2012-2014 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -67,18 +67,18 @@ class test_pfb_channelizer(gr_unittest.TestCase): self.tb.connect((s2ss,i), (pfb,i)) self.tb.connect((pfb, i), snks[i]) - self.tb.run() + self.tb.run() Ntest = 50 L = len(snks[0].data()) # Adjusted phase rotations for data - p0 = -2*math.pi * 0 / M - p1 = -2*math.pi * 1 / M - p2 = -2*math.pi * 2 / M - p3 = -2*math.pi * 3 / M - p4 = -2*math.pi * 4 / M - + p0 = 0.11058379158914133 + p1 = 4.5108246571401693 + p2 = 3.9739891674564594 + p3 = 2.2820531095511924 + p4 = 1.3782797467397869 + # Filter delay is the normal delay of each arm tpf = math.ceil(len(taps) / float(M)) delay = -(tpf - 1.0) / 2.0 diff --git a/gr-filter/python/filter/qa_pfb_decimator.py b/gr-filter/python/filter/qa_pfb_decimator.py index ea7a15570c..4366e85eec 100755 --- a/gr-filter/python/filter/qa_pfb_decimator.py +++ b/gr-filter/python/filter/qa_pfb_decimator.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright 2012,2013 Free Software Foundation, Inc. +# Copyright 2012-2014 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -29,7 +29,7 @@ def sig_source_c(samp_rate, freq, amp, N): 1j*math.sin(2.*math.pi*freq*x), t) return y -def run_test(tb, channel): +def run_test(tb, channel, fft_rotate, fft_filter): N = 1000 # number of samples to use M = 5 # Number of channels fs = 5000.0 # baseband sampling rate @@ -50,19 +50,24 @@ def run_test(tb, channel): tb.connect(signals[i], (add,i)) s2ss = blocks.stream_to_streams(gr.sizeof_gr_complex, M) - pfb = filter.pfb_decimator_ccf(M, taps, channel) + pfb = filter.pfb_decimator_ccf(M, taps, channel, fft_rotate, fft_filter) snk = blocks.vector_sink_c() tb.connect(add, s2ss) for i in xrange(M): tb.connect((s2ss,i), (pfb,i)) tb.connect(pfb, snk) - tb.run() + tb.run() L = len(snk.data()) - # Each channel is rotated by 2pi/M - phase = -2*math.pi * channel / M + # Adjusted phase rotations for data + phase = [ 0.11058476216852586, + 4.5108246571401693, + 3.9739891674564594, + 2.2820531095511924, + 1.3782797467397869] + phase = phase[channel] # Filter delay is the normal delay of each arm tpf = math.ceil(len(taps) / float(M)) @@ -90,33 +95,63 @@ class test_pfb_decimator(gr_unittest.TestCase): def test_000(self): Ntest = 50 - dst_data, expected_data = run_test(self.tb, 0) - - self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 4) + dst_data0, expected_data0 = run_test(self.tb, 0, False, False) + dst_data1, expected_data1 = run_test(self.tb, 0, False, True) + dst_data2, expected_data2 = run_test(self.tb, 0, True, False) + dst_data3, expected_data3 = run_test(self.tb, 0, True, True) + + self.assertComplexTuplesAlmostEqual(expected_data0[-Ntest:], dst_data0[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data1[-Ntest:], dst_data1[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data2[-Ntest:], dst_data2[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data3[-Ntest:], dst_data3[-Ntest:], 4) def test_001(self): Ntest = 50 - dst_data, expected_data = run_test(self.tb, 1) - - self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 4) + dst_data0, expected_data0 = run_test(self.tb, 1, False, False) + dst_data1, expected_data1 = run_test(self.tb, 1, False, True) + dst_data2, expected_data2 = run_test(self.tb, 1, True, False) + dst_data3, expected_data3 = run_test(self.tb, 1, True, True) + + self.assertComplexTuplesAlmostEqual(expected_data0[-Ntest:], dst_data0[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data1[-Ntest:], dst_data1[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data2[-Ntest:], dst_data2[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data3[-Ntest:], dst_data3[-Ntest:], 4) def test_002(self): Ntest = 50 - dst_data, expected_data = run_test(self.tb, 2) - - self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 4) + dst_data0, expected_data0 = run_test(self.tb, 2, False, False) + dst_data1, expected_data1 = run_test(self.tb, 2, False, True) + dst_data2, expected_data2 = run_test(self.tb, 2, True, False) + dst_data3, expected_data3 = run_test(self.tb, 2, True, True) + + self.assertComplexTuplesAlmostEqual(expected_data0[-Ntest:], dst_data0[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data1[-Ntest:], dst_data1[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data2[-Ntest:], dst_data2[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data3[-Ntest:], dst_data3[-Ntest:], 4) def test_003(self): Ntest = 50 - dst_data, expected_data = run_test(self.tb, 3) - - self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 4) + dst_data0, expected_data0 = run_test(self.tb, 3, False, False) + dst_data1, expected_data1 = run_test(self.tb, 3, False, True) + dst_data2, expected_data2 = run_test(self.tb, 3, True, False) + dst_data3, expected_data3 = run_test(self.tb, 3, True, True) + + self.assertComplexTuplesAlmostEqual(expected_data0[-Ntest:], dst_data0[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data1[-Ntest:], dst_data1[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data2[-Ntest:], dst_data2[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data3[-Ntest:], dst_data3[-Ntest:], 4) def test_004(self): Ntest = 50 - dst_data, expected_data = run_test(self.tb, 4) - - self.assertComplexTuplesAlmostEqual(expected_data[-Ntest:], dst_data[-Ntest:], 4) + dst_data0, expected_data0 = run_test(self.tb, 4, False, False) + dst_data1, expected_data1 = run_test(self.tb, 4, False, True) + dst_data2, expected_data2 = run_test(self.tb, 4, True, False) + dst_data3, expected_data3 = run_test(self.tb, 4, True, True) + + self.assertComplexTuplesAlmostEqual(expected_data0[-Ntest:], dst_data0[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data1[-Ntest:], dst_data1[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data2[-Ntest:], dst_data2[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected_data3[-Ntest:], dst_data3[-Ntest:], 4) if __name__ == '__main__': gr_unittest.run(test_pfb_decimator, "test_pfb_decimator.xml") diff --git a/gr-filter/python/filter/qa_pfb_interpolator.py b/gr-filter/python/filter/qa_pfb_interpolator.py index 33c22675af..b7ed4feef6 100755 --- a/gr-filter/python/filter/qa_pfb_interpolator.py +++ b/gr-filter/python/filter/qa_pfb_interpolator.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright 2012,2013 Free Software Foundation, Inc. +# Copyright 2012-2014 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -42,9 +42,9 @@ class test_pfb_interpolator(gr_unittest.TestCase): N = 1000 # number of samples to use M = 5 # Number of channels fs = 1000 # baseband sampling rate - ifs = M*fs # input samp rate to decimator + ofs = M*fs # output samp rate of interpolator - taps = filter.firdes.low_pass_2(M, ifs, fs/2, fs/10, + taps = filter.firdes.low_pass_2(M, ofs, fs/4, fs/10, attenuation_dB=80, window=filter.firdes.WIN_BLACKMAN_hARRIS) @@ -57,20 +57,16 @@ class test_pfb_interpolator(gr_unittest.TestCase): self.tb.connect(signal, pfb) self.tb.connect(pfb, snk) - self.tb.run() + self.tb.run() Ntest = 50 L = len(snk.data()) - # Can only get channel 0 out; no phase rotation - phase = 0 + # Phase rotation through the filters + phase = 4.8870112969978994 - # Calculate the filter delay - delay = -(len(taps) - 1) / 2.0 - (M-1) - delay = int(delay) - - # Create a time scale that's delayed to match the filter delay - t = map(lambda x: float(x)/ifs, xrange(delay, L+delay)) + # Create a time scale + t = map(lambda x: float(x)/ofs, xrange(0, L)) # Create known data as complex sinusoids for the baseband freq # of the extracted channel is due to decimator output order. diff --git a/gr-filter/python/filter/qa_pfb_synthesizer.py b/gr-filter/python/filter/qa_pfb_synthesizer.py index baba904088..0b3f8b27a2 100755 --- a/gr-filter/python/filter/qa_pfb_synthesizer.py +++ b/gr-filter/python/filter/qa_pfb_synthesizer.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright 2012,2013 Free Software Foundation, Inc. +# Copyright 2012-2014 Free Software Foundation, Inc. # # This file is part of GNU Radio # @@ -62,7 +62,14 @@ class test_pfb_synthesizer(gr_unittest.TestCase): self.tb.connect(pfb, snk) - self.tb.run() + self.tb.run() + + # Adjusted phase rotations for data + p0 = 0 + p1 = 2*math.pi*1.0/M + p2 = 2*math.pi*2.0/M + p3 = 2*math.pi*3.0/M + p4 = 2*math.pi*4.0/M Ntest = 1000 L = len(snk.data()) @@ -73,20 +80,19 @@ class test_pfb_synthesizer(gr_unittest.TestCase): freqs = [-2200, -1100, 0, 1100, 2200] expected_data = len(t)*[0,] for i in xrange(len(t)): - expected_data[i] = math.cos(2.*math.pi*freqs[0]*t[i]) + \ - 1j*math.sin(2.*math.pi*freqs[0]*t[i]) + \ - math.cos(2.*math.pi*freqs[1]*t[i]) + \ - 1j*math.sin(2.*math.pi*freqs[1]*t[i]) + \ - math.cos(2.*math.pi*freqs[2]*t[i]) + \ - 1j*math.sin(2.*math.pi*freqs[2]*t[i]) + \ - math.cos(2.*math.pi*freqs[3]*t[i]) + \ - 1j*math.sin(2.*math.pi*freqs[3]*t[i]) + \ - math.cos(2.*math.pi*freqs[4]*t[i]) + \ - 1j*math.sin(2.*math.pi*freqs[4]*t[i]) + expected_data[i] = math.cos(2.*math.pi*freqs[0]*t[i] + p3) + \ + 1j*math.sin(2.*math.pi*freqs[0]*t[i] + p3) + \ + math.cos(2.*math.pi*freqs[1]*t[i] + p4) + \ + 1j*math.sin(2.*math.pi*freqs[1]*t[i] + p4) + \ + math.cos(2.*math.pi*freqs[2]*t[i] + p0) + \ + 1j*math.sin(2.*math.pi*freqs[2]*t[i] + p0) + \ + math.cos(2.*math.pi*freqs[3]*t[i] + p1) + \ + 1j*math.sin(2.*math.pi*freqs[3]*t[i] + p1) + \ + math.cos(2.*math.pi*freqs[4]*t[i] + p2) + \ + 1j*math.sin(2.*math.pi*freqs[4]*t[i] + p2) dst_data = snk.data() - offset = 25 - self.assertComplexTuplesAlmostEqual(expected_data[2000-offset:2000-offset+Ntest], + self.assertComplexTuplesAlmostEqual(expected_data[2000:2000+Ntest], dst_data[2000:2000+Ntest], 4) if __name__ == '__main__': diff --git a/gr-filter/swig/filter_swig.i b/gr-filter/swig/filter_swig.i index bc9f204850..3eb3f3c634 100644 --- a/gr-filter/swig/filter_swig.i +++ b/gr-filter/swig/filter_swig.i @@ -40,6 +40,7 @@ #include "gnuradio/filter/fir_filter_fsf.h" #include "gnuradio/filter/fir_filter_scc.h" #include "gnuradio/filter/fft_filter_ccc.h" +#include "gnuradio/filter/fft_filter_ccf.h" #include "gnuradio/filter/fft_filter_fff.h" #include "gnuradio/filter/fractional_interpolator_cc.h" #include "gnuradio/filter/fractional_interpolator_ff.h" @@ -88,6 +89,7 @@ %include "gnuradio/filter/fir_filter_fsf.h" %include "gnuradio/filter/fir_filter_scc.h" %include "gnuradio/filter/fft_filter_ccc.h" +%include "gnuradio/filter/fft_filter_ccf.h" %include "gnuradio/filter/fft_filter_fff.h" %include "gnuradio/filter/fractional_interpolator_cc.h" %include "gnuradio/filter/fractional_interpolator_ff.h" @@ -133,6 +135,7 @@ GR_SWIG_BLOCK_MAGIC2(filter, fir_filter_fff); GR_SWIG_BLOCK_MAGIC2(filter, fir_filter_fsf); GR_SWIG_BLOCK_MAGIC2(filter, fir_filter_scc); GR_SWIG_BLOCK_MAGIC2(filter, fft_filter_ccc); +GR_SWIG_BLOCK_MAGIC2(filter, fft_filter_ccf); GR_SWIG_BLOCK_MAGIC2(filter, fft_filter_fff); GR_SWIG_BLOCK_MAGIC2(filter, fractional_interpolator_cc); GR_SWIG_BLOCK_MAGIC2(filter, fractional_interpolator_ff); |