diff options
| -rw-r--r-- | gr-filter/examples/CMakeLists.txt | 2 | ||||
| -rwxr-xr-x | gr-filter/examples/channelize.py | 195 | ||||
| -rwxr-xr-x | gr-filter/examples/chirp_channelize.py | 204 | ||||
| -rw-r--r-- | gr-filter/grc/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | gr-filter/grc/filter_block_tree.xml | 1 | ||||
| -rw-r--r-- | gr-filter/grc/pfb_channelizer.xml | 62 | ||||
| -rw-r--r-- | gr-filter/include/filter/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | gr-filter/include/filter/pfb_channelizer_ccf.h | 205 | ||||
| -rw-r--r-- | gr-filter/lib/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | gr-filter/lib/hilbert_fc_impl.h | 1 | ||||
| -rw-r--r-- | gr-filter/lib/pfb_channelizer_ccf_impl.cc | 240 | ||||
| -rw-r--r-- | gr-filter/lib/pfb_channelizer_ccf_impl.h | 72 | ||||
| -rw-r--r-- | gr-filter/python/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | gr-filter/python/__init__.py | 2 | ||||
| -rw-r--r-- | gr-filter/python/pfb.py | 74 | ||||
| -rw-r--r-- | gr-filter/python/qa_pfb_channelizer.py | 104 | ||||
| -rw-r--r-- | gr-filter/swig/filter_swig.i | 3 |
17 files changed, 1167 insertions, 2 deletions
diff --git a/gr-filter/examples/CMakeLists.txt b/gr-filter/examples/CMakeLists.txt index 12a3c42e0..2eca31726 100644 --- a/gr-filter/examples/CMakeLists.txt +++ b/gr-filter/examples/CMakeLists.txt @@ -24,6 +24,8 @@ GR_PYTHON_INSTALL(PROGRAMS fir_filter_fff.py fir_filter_ccc.py fft_filter_ccc.py + channelize.py + chirp_channelize.py DESTINATION ${GR_PKG_FILTER_EXAMPLES_DIR} COMPONENT "filter_python" ) diff --git a/gr-filter/examples/channelize.py b/gr-filter/examples/channelize.py new file mode 100755 index 000000000..affce7b57 --- /dev/null +++ b/gr-filter/examples/channelize.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python +# +# Copyright 2009,2012 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +from gnuradio import filter +import sys, time + +try: + import scipy + from scipy import fftpack +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab + from pylab import mlab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 2000000 # number of samples to use + self._fs = 1000 # initial sampling rate + self._M = M = 9 # Number of channels to channelize + self._ifs = M*self._fs # initial sampling rate + + # Create a set of taps for the PFB channelizer + self._taps = gr.firdes.low_pass_2(1, self._ifs, 475.50, 50, + attenuation_dB=100, + window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._M)) + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + # Create a set of signals at different frequencies + # freqs lists the frequencies of the signals that get stored + # in the list "signals", which then get summed together + self.signals = list() + self.add = gr.add_cc() + freqs = [-70, -50, -30, -10, 10, 20, 40, 60, 80] + for i in xrange(len(freqs)): + f = freqs[i] + (M/2-M+i+1)*self._fs + self.signals.append(gr.sig_source_c(self._ifs, gr.GR_SIN_WAVE, f, 1)) + self.connect(self.signals[i], (self.add,i)) + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the channelizer filter + self.pfb = filter.pfb.channelizer_ccf(self._M, self._taps, 1) + + # Construct a vector sink for the input signal to the channelizer + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.add, self.head, self.pfb) + self.connect(self.add, self.snk_i) + + # Use this to play with the channel mapping + #self.pfb.set_channel_map([5,6,7,8,0,1,2,3,4]) + + # Create a vector sink for each of M output channels of the filter and connect it + self.snks = list() + for i in xrange(self._M): + self.snks.append(gr.vector_sink_c()) + self.connect((self.pfb, i), self.snks[i]) + + +def main(): + tstart = time.time() + + tb = pfb_top_block() + tb.run() + + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig_in = pylab.figure(1, figsize=(16,9), facecolor="w") + fig1 = pylab.figure(2, figsize=(16,9), facecolor="w") + fig2 = pylab.figure(3, figsize=(16,9), facecolor="w") + + Ns = 1000 + Ne = 10000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._ifs + + # Plot the input signal on its own figure + d = tb.snk_i.data()[Ns:Ne] + spin_f = fig_in.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(X)) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + pin_f = spin_f.plot(f_in, X_in, "b") + spin_f.set_xlim([min(f_in), max(f_in)+1]) + spin_f.set_ylim([-200.0, 50.0]) + + spin_f.set_title("Input Signal", weight="bold") + spin_f.set_xlabel("Frequency (Hz)") + spin_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + spin_t = fig_in.add_subplot(2, 1, 2) + pin_t = spin_t.plot(t_in, x_in.real, "b") + pin_t = spin_t.plot(t_in, x_in.imag, "r") + + spin_t.set_xlabel("Time (s)") + spin_t.set_ylabel("Amplitude") + + Ncols = int(scipy.floor(scipy.sqrt(tb._M))) + Nrows = int(scipy.floor(tb._M / Ncols)) + if(tb._M % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = tb._fs + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + for i in xrange(len(tb.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = tb.snks[i].data()[Ns:Ne] + + sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i) + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(X)) + f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size)) + p2_f = sp1_f.plot(f_o, X_o, "b") + sp1_f.set_xlim([min(f_o), max(f_o)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title(("Channel %d" % i), weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i) + p2_o = sp2_o.plot(t_o, x_o.real, "b") + p2_o = sp2_o.plot(t_o, x_o.imag, "r") + sp2_o.set_xlim([min(t_o), max(t_o)+1]) + sp2_o.set_ylim([-2, 2]) + + sp2_o.set_title(("Channel %d" % i), weight="bold") + sp2_o.set_xlabel("Time (s)") + sp2_o.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gr-filter/examples/chirp_channelize.py b/gr-filter/examples/chirp_channelize.py new file mode 100755 index 000000000..99bf3b6e2 --- /dev/null +++ b/gr-filter/examples/chirp_channelize.py @@ -0,0 +1,204 @@ +#!/usr/bin/env python +# +# Copyright 2009,2012 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +from gnuradio import filter +import sys, time + +try: + import scipy + from scipy import fftpack +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab + from pylab import mlab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 200000 # number of samples to use + self._fs = 9000 # initial sampling rate + self._M = 9 # Number of channels to channelize + + # Create a set of taps for the PFB channelizer + self._taps = gr.firdes.low_pass_2(1, self._fs, 500, 20, + attenuation_dB=10, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._M)) + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + repeated = True + if(repeated): + self.vco_input = gr.sig_source_f(self._fs, gr.GR_SIN_WAVE, 0.25, 110) + else: + amp = 100 + data = scipy.arange(0, amp, amp/float(self._N)) + self.vco_input = gr.vector_source_f(data, False) + + # Build a VCO controlled by either the sinusoid or single chirp tone + # Then convert this to a complex signal + self.vco = gr.vco_f(self._fs, 225, 1) + self.f2c = gr.float_to_complex() + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the channelizer filter + self.pfb = filter.pfb.channelizer_ccf(self._M, self._taps) + + # Construct a vector sink for the input signal to the channelizer + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.vco_input, self.vco, self.f2c) + self.connect(self.f2c, self.head, self.pfb) + self.connect(self.f2c, self.snk_i) + + # Create a vector sink for each of M output channels of the filter and connect it + self.snks = list() + for i in xrange(self._M): + self.snks.append(gr.vector_sink_c()) + self.connect((self.pfb, i), self.snks[i]) + + +def main(): + tstart = time.time() + + tb = pfb_top_block() + tb.run() + + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig_in = pylab.figure(1, figsize=(16,9), facecolor="w") + fig1 = pylab.figure(2, figsize=(16,9), facecolor="w") + fig2 = pylab.figure(3, figsize=(16,9), facecolor="w") + fig3 = pylab.figure(4, figsize=(16,9), facecolor="w") + + Ns = 650 + Ne = 20000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._fs + + # Plot the input signal on its own figure + d = tb.snk_i.data()[Ns:Ne] + spin_f = fig_in.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + pin_f = spin_f.plot(f_in, X_in, "b") + spin_f.set_xlim([min(f_in), max(f_in)+1]) + spin_f.set_ylim([-200.0, 50.0]) + + spin_f.set_title("Input Signal", weight="bold") + spin_f.set_xlabel("Frequency (Hz)") + spin_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + spin_t = fig_in.add_subplot(2, 1, 2) + pin_t = spin_t.plot(t_in, x_in.real, "b") + pin_t = spin_t.plot(t_in, x_in.imag, "r") + + spin_t.set_xlabel("Time (s)") + spin_t.set_ylabel("Amplitude") + + Ncols = int(scipy.floor(scipy.sqrt(tb._M))) + Nrows = int(scipy.floor(tb._M / Ncols)) + if(tb._M % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = tb._fs / tb._M + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + for i in xrange(len(tb.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = tb.snks[i].data()[Ns:Ne] + + sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i) + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(X)) + f_o = freq + p2_f = sp1_f.plot(f_o, X_o, "b") + sp1_f.set_xlim([min(f_o), max(f_o)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title(("Channel %d" % i), weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i) + p2_o = sp2_o.plot(t_o, x_o.real, "b") + p2_o = sp2_o.plot(t_o, x_o.imag, "r") + sp2_o.set_xlim([min(t_o), max(t_o)+1]) + sp2_o.set_ylim([-2, 2]) + + sp2_o.set_title(("Channel %d" % i), weight="bold") + sp2_o.set_xlabel("Time (s)") + sp2_o.set_ylabel("Amplitude") + + + sp3 = fig3.add_subplot(1,1,1) + p3 = sp3.plot(t_o, x_o.real) + sp3.set_xlim([min(t_o), max(t_o)+1]) + sp3.set_ylim([-2, 2]) + + sp3.set_title("All Channels") + sp3.set_xlabel("Time (s)") + sp3.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gr-filter/grc/CMakeLists.txt b/gr-filter/grc/CMakeLists.txt index 6ace8a8b4..1c69e47c3 100644 --- a/gr-filter/grc/CMakeLists.txt +++ b/gr-filter/grc/CMakeLists.txt @@ -23,6 +23,7 @@ install(FILES fft_filter_xxx.xml filter_delay_fc.xml hilbert_fc.xml + pfb_channelizer.xml DESTINATION ${GRC_BLOCKS_DIR} COMPONENT "filter_python" ) diff --git a/gr-filter/grc/filter_block_tree.xml b/gr-filter/grc/filter_block_tree.xml index fbc4e7cb5..0c685d9f6 100644 --- a/gr-filter/grc/filter_block_tree.xml +++ b/gr-filter/grc/filter_block_tree.xml @@ -34,5 +34,6 @@ <block>fft_filter_xxx</block> <block>filter_delay_fc</block> <block>hilbert_fc</block> + <block>pfb_channelizer_ccf</block> </cat> </cat> diff --git a/gr-filter/grc/pfb_channelizer.xml b/gr-filter/grc/pfb_channelizer.xml new file mode 100644 index 000000000..114abc0f0 --- /dev/null +++ b/gr-filter/grc/pfb_channelizer.xml @@ -0,0 +1,62 @@ +<?xml version="1.0"?> +<!-- +################################################### +##Polyphase Channelizer +################################################### + --> +<block> + <name>Polyphase Channelizer</name> + <key>pfb_channelizer_ccf</key> + <import>from gnuradio import filter</import> + <import>from gnuradio.filter import firdes</import> + <make>filter.pfb.channelizer_ccf( + $nchans, + $taps, + $osr, + $atten) +self.$(id).set_channel_map($ch_map) + </make> + <!-- Set taps not implemented yet + <callback>set_taps($taps)</callback> + --> + <callback>set_channel_map($ch_map)</callback> + + <param> + <name>Channels</name> + <key>nchans</key> + <type>int</type> + </param> + <param> + <name>Taps</name> + <key>taps</key> + <value>None</value> + <type>real_vector</type> + </param> + <param> + <name>Over Sample Ratio</name> + <key>osr</key> + <value>1.0</value> + <type>real</type> + </param> + <param> + <name>Attenuation</name> + <key>atten</key> + <value>100</value> + <type>real</type> + </param> + <param> + <name>Channel Map</name> + <key>ch_map</key> + <value>[]</value> + <type>int_vector</type> + </param> + <sink> + <name>in</name> + <type>complex</type> + </sink> + <source> + <name>out</name> + <type>complex</type> + <nports>$nchans</nports> + </source> +</block> diff --git a/gr-filter/include/filter/CMakeLists.txt b/gr-filter/include/filter/CMakeLists.txt index d5452a130..c64178696 100644 --- a/gr-filter/include/filter/CMakeLists.txt +++ b/gr-filter/include/filter/CMakeLists.txt @@ -87,6 +87,7 @@ install(FILES fft_filter_ccc.h fft_filter_fff.h hilbert_fc.h + pfb_channelizer_ccf.h DESTINATION ${GR_INCLUDE_DIR}/gnuradio/filter COMPONENT "fft_devel" ) diff --git a/gr-filter/include/filter/pfb_channelizer_ccf.h b/gr-filter/include/filter/pfb_channelizer_ccf.h new file mode 100644 index 000000000..ab79f7036 --- /dev/null +++ b/gr-filter/include/filter/pfb_channelizer_ccf.h @@ -0,0 +1,205 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009,2010,2012 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_PFB_CHANNELIZER_CCF_H +#define INCLUDED_FILTER_PFB_CHANNELIZER_CCF_H + +#include <filter/api.h> +#include <gr_block.h> +#include <gruel/thread.h> + +namespace gr { + namespace filter { + + /*! + * \class gr_pfb_channelizer_ccf + * + * \brief Polyphase filterbank channelizer with + * gr_complex input, gr_complex output and float taps + * + * \ingroup filter_blk + * \ingroup pfb_blk + * + * This block takes in complex inputs and channelizes it to <EM>M</EM> + * channels of equal bandwidth. Each of the resulting channels is + * decimated to the new rate that is the input sampling rate + * <EM>fs</EM> divided by the number of channels, <EM>M</EM>. + * + * The PFB channelizer code takes the taps generated above and builds + * a set of filters. The set contains <EM>M</EM> number of filters + * and each filter contains ceil(taps.size()/decim) number of taps. + * Each tap from the filter prototype is sequentially inserted into + * the next filter. When all of the input taps are used, the remaining + * filters in the filterbank are filled out with 0's to make sure each + * filter has the same number of taps. + * + * Each filter operates using the gr_fir filter classs of GNU Radio, + * which takes the input stream at <EM>i</EM> and performs the inner + * product calculation to <EM>i+(n-1)</EM> where <EM>n</EM> is the + * number of filter taps. To efficiently handle this in the GNU Radio + * structure, each filter input must come from its own input + * stream. So the channelizer must be provided with <EM>M</EM> streams + * where the input stream has been deinterleaved. This is most easily + * done using the gr_stream_to_streams block. + * + * The output is then produced as a vector, where index <EM>i</EM> in + * the vector is the next sample from the <EM>i</EM>th channel. This + * is most easily handled by sending the output to a + * gr_vector_to_streams block to handle the conversion and passing + * <EM>M</EM> streams out. + * + * The input and output formatting is done using a hier_block2 called + * pfb_channelizer_ccf. This can take in a single stream and outputs + * <EM>M</EM> streams based on the behavior described above. + * + * The filter's taps should be based on the input sampling rate. + * + * For example, using the GNU Radio's firdes utility to building + * filters, we build a low-pass filter with a sampling rate of + * <EM>fs</EM>, a 3-dB bandwidth of <EM>BW</EM> and a transition + * bandwidth of <EM>TB</EM>. We can also specify the out-of-band + * attenuation to use, <EM>ATT</EM>, and the filter window + * function (a Blackman-harris window in this case). The first input + * is the gain of the filter, which we specify here as unity. + * + * <B><EM>self._taps = gr.firdes.low_pass_2(1, fs, BW, TB, + * attenuation_dB=ATT, window=gr.firdes.WIN_BLACKMAN_hARRIS)</EM></B> + * + * The 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. + * + * <B><EM>f. harris, "Multirate Signal Processing for Communication + * Systems," Upper Saddle River, NJ: Prentice Hall, Inc. 2004.</EM></B> + * + */ + + class FILTER_API pfb_channelizer_ccf : virtual public gr_block + { + public: + // gr::filter::pfb_channelizer_ccf::sptr + typedef boost::shared_ptr<pfb_channelizer_ccf> sptr; + + /*! + * 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 + * rateis 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. + */ + static FILTER_API sptr make(unsigned int numchans, + const std::vector<float> &taps, + float oversample_rate); + + /*! + * Resets the filterbank's filter taps with the new prototype filter + * \param taps (vector/list of floats) The prototype filter to populate the filterbank. + */ + virtual void set_taps(const std::vector<float> &taps) = 0; + + /*! + * Print all of the filterbank taps to screen. + */ + virtual void print_taps() = 0; + + /*! + * Return a vector<vector<>> of the filterbank taps + */ + virtual std::vector<std::vector<float> > taps() const = 0; + + /*! + * Set the channel map. Channels are numbers as: + * + * N/2+1 | ... | N-1 | 0 | 1 | 2 | ... | N/2 + * <------------------- 0 --------------------> + * freq + * + * So output stream 0 comes from channel 0, etc. Setting a new + * channel map allows the user to specify which channel in frequency + * he/she wants to got to which output stream. + * + * The map should have the same number of elements as the number + * of output connections from the block. The minimum value of + * the map is 0 (for the 0th channel) and the maximum number is + * N-1 where N is the number of channels. + * + * We specify M as the number of output connections made where M + * <= N, so only M out of N channels are driven to an output + * stream. The number of items in the channel map should be at + * least M long. If there are more channels specified, any value + * in the map over M-1 will be ignored. If the size of the map + * is less than M the behavior is unknown (we don't wish to + * check every entry into the work function). + * + * This means that if the channelizer is splitting the signal up + * into N channels but only M channels are specified in the map + * (where M <= N), then M output streams must be connected and + * the map and the channel numbers used must be less than + * N-1. Output channel number can be reused, too. By default, + * 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. + */ + virtual std::vector<int> channel_map() const = 0; + }; + + } /* namespace filter */ +} /* namespace gr */ + +#endif /* INCLUDED_FILTER_PFB_CHANNELIZER_CCF_H */ diff --git a/gr-filter/lib/CMakeLists.txt b/gr-filter/lib/CMakeLists.txt index e1fae56a3..5930ec9bb 100644 --- a/gr-filter/lib/CMakeLists.txt +++ b/gr-filter/lib/CMakeLists.txt @@ -117,6 +117,7 @@ list(APPEND filter_sources fft_filter_ccc_impl.cc fft_filter_fff_impl.cc hilbert_fc_impl.cc + pfb_channelizer_ccf_impl.cc ) list(APPEND filter_libs diff --git a/gr-filter/lib/hilbert_fc_impl.h b/gr-filter/lib/hilbert_fc_impl.h index de4d754b4..d2b41b573 100644 --- a/gr-filter/lib/hilbert_fc_impl.h +++ b/gr-filter/lib/hilbert_fc_impl.h @@ -25,7 +25,6 @@ #include <filter/hilbert_fc.h> #include <filter/fir_filter.h> -#include <gr_io_signature.h> #include <gr_types.h> namespace gr { diff --git a/gr-filter/lib/pfb_channelizer_ccf_impl.cc b/gr-filter/lib/pfb_channelizer_ccf_impl.cc new file mode 100644 index 000000000..4006db789 --- /dev/null +++ b/gr-filter/lib/pfb_channelizer_ccf_impl.cc @@ -0,0 +1,240 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009,2010,2012 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 "pfb_channelizer_ccf_impl.h" +#include <gr_io_signature.h> +#include <cstdio> + +namespace gr { + namespace filter { + + pfb_channelizer_ccf::sptr pfb_channelizer_ccf::make(unsigned int numchans, + const std::vector<float> &taps, + float oversample_rate) + { + return gnuradio::get_initial_sptr(new pfb_channelizer_ccf_impl(numchans, taps, + oversample_rate)); + } + + pfb_channelizer_ccf_impl::pfb_channelizer_ccf_impl(unsigned int numchans, + 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, numchans, sizeof(gr_complex))), + d_updated(false), d_numchans(numchans), d_oversample_rate(oversample_rate) + { + // 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 fltp = modf(numchans / oversample_rate, &intp); + if(fltp != 0.0) + throw std::invalid_argument("pfb_channelizer: oversample rate must be N/i for i in [1, N]"); + + set_relative_rate(1.0/intp); + + d_filters = std::vector<kernel::fir_filter_ccf*>(d_numchans); + d_channel_map.resize(d_numchans); + + // Create an 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_ccf(1, vtaps); + d_channel_map[i] = i; + } + + // 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_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); + } + + pfb_channelizer_ccf_impl::~pfb_channelizer_ccf_impl() + { + delete [] d_idxlut; + delete d_fft; + + for(unsigned int i = 0; i < d_numchans; i++) { + delete d_filters[i]; + } + } + + void + pfb_channelizer_ccf_impl::set_taps (const std::vector<float> &taps) + { + gruel::scoped_lock guard(d_mutex); + unsigned int i,j; + + unsigned int ntaps = taps.size(); + d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_numchans); + + // Create d_numchan vectors to store each channel's taps + d_taps.resize(d_numchans); + + // Make a vector of the taps plus fill it out with 0's to fill + // each polyphase filter with exactly d_taps_per_filter + std::vector<float> tmp_taps; + tmp_taps = taps; + while((float)(tmp_taps.size()) < d_numchans*d_taps_per_filter) { + tmp_taps.push_back(0.0); + } + + // Partition the filter + for(i = 0; i < d_numchans; i++) { + // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out + 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]); + } + + // Set the history to ensure enough input items for each filter + set_history (d_taps_per_filter+1); + + d_updated = true; + } + + void + pfb_channelizer_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"); + } + } + + std::vector< std::vector<float> > + pfb_channelizer_ccf_impl::taps() const + { + return d_taps; + } + + void + pfb_channelizer_ccf_impl::set_channel_map(const std::vector<int> &map) + { + gruel::scoped_lock guard(d_mutex); + + if(map.size() > 0) { + 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_numchans) || (min < 0)) { + throw std::invalid_argument("pfb_channelizer_ccf_impl::set_channel_map: map range out of bounds.\n"); + } + d_channel_map = map; + } + } + + std::vector<int> + pfb_channelizer_ccf_impl::channel_map() const + { + return d_channel_map; + } + + int + pfb_channelizer_ccf_impl::general_work(int noutput_items, + gr_vector_int &ninput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) + { + gruel::scoped_lock guard(d_mutex); + + 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. + } + + size_t noutputs = output_items.size(); + + int n=1, i=-1, j=0, oo=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--; + } + + 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(); + + // Send to output channels + for(unsigned int nn = 0; nn < noutputs; nn++) { + out = (gr_complex*)output_items[nn]; + out[oo] = d_fft->get_outbuf()[d_channel_map[nn]]; + } + oo++; + } + + consume_each(toconsume); + return noutput_items; + } + + } /* namespace filter */ +} /* namespace gr */ diff --git a/gr-filter/lib/pfb_channelizer_ccf_impl.h b/gr-filter/lib/pfb_channelizer_ccf_impl.h new file mode 100644 index 000000000..ee9327663 --- /dev/null +++ b/gr-filter/lib/pfb_channelizer_ccf_impl.h @@ -0,0 +1,72 @@ +/* -*- c++ -*- */ +/* + * Copyright 2009,2010,2012 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_PFB_CHANNELIZER_CCF_IMPL_H +#define INCLUDED_FILTER_PFB_CHANNELIZER_CCF_IMPL_H + +#include <filter/pfb_channelizer_ccf.h> +#include <filter/fir_filter.h> +#include <fft/fft.h> + +namespace gr { + namespace filter { + + class FILTER_API pfb_channelizer_ccf_impl : public pfb_channelizer_ccf + { + private: + bool d_updated; + unsigned int d_numchans; + float d_oversample_rate; + std::vector<kernel::fir_filter_ccf*> d_filters; + std::vector< std::vector<float> > d_taps; + unsigned int d_taps_per_filter; + fft::fft_complex *d_fft; + int *d_idxlut; + int d_rate_ratio; + int d_output_multiple; + std::vector<int> d_channel_map; + gruel::mutex d_mutex; // mutex to protect set/work access + + public: + pfb_channelizer_ccf_impl(unsigned int numchans, + const std::vector<float> &taps, + float oversample_rate); + + ~pfb_channelizer_ccf_impl(); + + void set_taps(const std::vector<float> &taps); + void print_taps(); + std::vector<std::vector<float> > taps() const; + + void set_channel_map(const std::vector<int> &map); + std::vector<int> channel_map() const; + + int general_work(int noutput_items, + gr_vector_int &ninput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + }; + + } /* namespace filter */ +} /* namespace gr */ + +#endif diff --git a/gr-filter/python/CMakeLists.txt b/gr-filter/python/CMakeLists.txt index e52cac759..99a1aa3a9 100644 --- a/gr-filter/python/CMakeLists.txt +++ b/gr-filter/python/CMakeLists.txt @@ -23,6 +23,7 @@ include(GrPython) GR_PYTHON_INSTALL( FILES __init__.py + pfb.py DESTINATION ${GR_PYTHON_DIR}/gnuradio/filter COMPONENT "filter_python" ) diff --git a/gr-filter/python/__init__.py b/gr-filter/python/__init__.py index 56dd2dc5a..90b5ce0a4 100644 --- a/gr-filter/python/__init__.py +++ b/gr-filter/python/__init__.py @@ -25,4 +25,4 @@ processing blocks for FILTER and related functions. ''' from filter_swig import * - +import pfb diff --git a/gr-filter/python/pfb.py b/gr-filter/python/pfb.py new file mode 100644 index 000000000..9c7e18e31 --- /dev/null +++ b/gr-filter/python/pfb.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python +# +# Copyright 2009,2010,2012 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, optfir +import filter_swig as filter + +class channelizer_ccf(gr.hier_block2): + ''' + Make a Polyphase Filter channelizer (complex in, complex out, floating-point taps) + + This simplifies the interface by allowing a single input stream to connect to this block. + It will then output a stream for each channel. + ''' + def __init__(self, numchans, taps=None, oversample_rate=1, atten=100): + gr.hier_block2.__init__(self, "pfb_channelizer_ccf", + gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature + gr.io_signature(numchans, numchans, gr.sizeof_gr_complex)) # Output signature + + self._nchans = numchans + self._oversample_rate = oversample_rate + + if taps is not None: + self._taps = taps + else: + # Create a filter that covers the full bandwidth of the input signal + bw = 0.4 + tb = 0.2 + ripple = 0.1 + made = False + while not made: + try: + self._taps = optfir.low_pass(1, self._nchans, bw, bw+tb, ripple, atten) + made = True + except RuntimeError: + ripple += 0.01 + made = False + print("Warning: set ripple to %.4f dB. If this is a problem, adjust the attenuation or create your own filter taps." % (ripple)) + + # Build in an exit strategy; if we've come this far, it ain't working. + if(ripple >= 1.0): + raise RuntimeError("optfir could not generate an appropriate filter.") + + self.s2ss = gr.stream_to_streams(gr.sizeof_gr_complex, self._nchans) + self.pfb = filter.pfb_channelizer_ccf(self._nchans, self._taps, + self._oversample_rate) + self.connect(self, self.s2ss) + + for i in xrange(self._nchans): + self.connect((self.s2ss,i), (self.pfb,i)) + self.connect((self.pfb,i), (self,i)) + + def set_channel_map(self, newmap): + self.pfb.set_channel_map(newmap) + + diff --git a/gr-filter/python/qa_pfb_channelizer.py b/gr-filter/python/qa_pfb_channelizer.py new file mode 100644 index 000000000..b52c80e8b --- /dev/null +++ b/gr-filter/python/qa_pfb_channelizer.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python +# +# Copyright 2012 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, gr_unittest +import filter_swig as filter +import math + +class test_pfb_channelizer(gr_unittest.TestCase): + + def setUp(self): + self.tb = gr.top_block() + + def tearDown(self): + self.tb = None + + def test_000(self): + N = 1000 # number of samples to use + M = 5 # Number of channels to channelize + fs = 1000 # baseband sampling rate + ifs = M*fs # input samp rate to channelizer + + taps = gr.firdes.low_pass_2(1, ifs, 500, 50, + attenuation_dB=80, + window=gr.firdes.WIN_BLACKMAN_hARRIS) + + signals = list() + add = gr.add_cc() + freqs = [-200, -100, 0, 100, 200] + for i in xrange(len(freqs)): + f = freqs[i] + (M/2-M+i+1)*fs + signals.append(gr.sig_source_c(ifs, gr.GR_SIN_WAVE, f, 1)) + self.tb.connect(signals[i], (add,i)) + + head = gr.head(gr.sizeof_gr_complex, N) + s2ss = gr.stream_to_streams(gr.sizeof_gr_complex, M) + pfb = filter.pfb_channelizer_ccf(M, taps, 1) + + self.tb.connect(add, head, s2ss) + + snks = list() + for i in xrange(M): + snks.append(gr.vector_sink_c()) + self.tb.connect((s2ss,i), (pfb,i)) + self.tb.connect((pfb, i), snks[i]) + + self.tb.run() + + Ntest = 50 + L = len(snks[0].data()) + t = map(lambda x: float(x)/fs, xrange(L)) + + # Adjusted phase rotations for data + p0 = 0 + p1 = 1.6335486 + p2 = -3.01609 + p3 = 3.01609 + p4 = -1.6335486 + + # Create known data as complex sinusoids at the different baseband freqs + # the different channel numbering is due to channelizer output order. + expected0_data = map(lambda x: math.cos(2.*math.pi*freqs[2]*x+p0) + \ + 1j*math.sin(2.*math.pi*freqs[2]*x+p0), t) + expected1_data = map(lambda x: math.cos(2.*math.pi*freqs[3]*x+p1) + \ + 1j*math.sin(2.*math.pi*freqs[3]*x+p1), t) + expected2_data = map(lambda x: math.cos(2.*math.pi*freqs[4]*x+p2) + \ + 1j*math.sin(2.*math.pi*freqs[4]*x+p2), t) + expected3_data = map(lambda x: math.cos(2.*math.pi*freqs[0]*x+p3) + \ + 1j*math.sin(2.*math.pi*freqs[0]*x+p3), t) + expected4_data = map(lambda x: math.cos(2.*math.pi*freqs[1]*x+p4) + \ + 1j*math.sin(2.*math.pi*freqs[1]*x+p4), t) + + dst0_data = snks[0].data() + dst1_data = snks[1].data() + dst2_data = snks[2].data() + dst3_data = snks[3].data() + dst4_data = snks[4].data() + + self.assertComplexTuplesAlmostEqual(expected0_data[-Ntest:], dst0_data[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected1_data[-Ntest:], dst1_data[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected2_data[-Ntest:], dst2_data[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected3_data[-Ntest:], dst3_data[-Ntest:], 4) + self.assertComplexTuplesAlmostEqual(expected4_data[-Ntest:], dst4_data[-Ntest:], 4) + +if __name__ == '__main__': + gr_unittest.run(test_pfb_channelizer, "test_pfb_channelizer.xml") diff --git a/gr-filter/swig/filter_swig.i b/gr-filter/swig/filter_swig.i index 4f0caa59c..cc15b5722 100644 --- a/gr-filter/swig/filter_swig.i +++ b/gr-filter/swig/filter_swig.i @@ -39,6 +39,7 @@ #include "filter/fft_filter_ccc.h" #include "filter/fft_filter_fff.h" #include "filter/hilbert_fc.h" +#include "filter/pfb_channelizer_ccf.h" %} %include "filter/firdes.h" @@ -52,6 +53,7 @@ %include "filter/fft_filter_ccc.h" %include "filter/fft_filter_fff.h" %include "filter/hilbert_fc.h" +%include "filter/pfb_channelizer_ccf.h" GR_SWIG_BLOCK_MAGIC2(filter, dc_blocker_cc); GR_SWIG_BLOCK_MAGIC2(filter, dc_blocker_ff); @@ -62,3 +64,4 @@ GR_SWIG_BLOCK_MAGIC2(filter, fir_filter_ccc); GR_SWIG_BLOCK_MAGIC2(filter, fft_filter_ccc); GR_SWIG_BLOCK_MAGIC2(filter, fft_filter_fff); GR_SWIG_BLOCK_MAGIC2(filter, hilbert_fc); +GR_SWIG_BLOCK_MAGIC2(filter, pfb_channelizer_ccf); |