GNU Radio 3.7.2 C++ API
Polyphase Filterbanks

Introduction

Polyphase filterbanks (PFB) are a very powerful set of filtering tools that can efficiently perform many multi-rate signal processing tasks. GNU Radio has a set of polyphase filterbank blocks to be used in all sorts of applications.

Usage

See the documentation for the individual blocks for details about what they can do and how they should be used. Furthermore, there are examples for these blocks in gr-filter/examples.

The main issue when using the PFB blocks is defining the prototype filter, which is passed to all of the blocks as a vector of taps. The taps from the prototype filter which get partitioned among the N channels of the channelizer.

An example of creating a set of filter taps for a PFB channelizer is found on line 49 of gr-filter/examples/channelizer.py and reproduced below. Notice that the sample rate is the sample rate at the input to the channelizer while the bandwidth and transition width are defined for the channel bandwidths. This makes a fairly long filter that is then split up between the N channels of the PFB.

self._fs = 9000 # input sample rate
self._M = 9 # Number of channels to channelize
self._taps = filter.firdes.low_pass_2(1, self._fs, 475.50, 50,
attenuation_dB=100,
window=filter.firdes.WIN_BLACKMAN_hARRIS)

In this example, the signal into the channelizer is sampled at 9 ksps (complex, so 9 kHz of bandwidth). The filter uses 9 channels, so each output channel will have a bandwidth and sample rate of 1 kHz. We want to pass most of the channel, so we define the channel bandwidth to be a low pass filter with a bandwidth of 475.5 Hz and a transition bandwidth of 50 Hz, but we have defined this using a sample rate of the original 9 kHz. The prototype filter has 819 taps to be divided up between the 9 channels, so each channel uses 91 taps. This is probably over-kill for a channelizer, and we could reduce the amount of taps per channel to a couple of dozen with no ill effects.

The basic rule when defining a set of taps for a PFB block is to think about the filter running at the highest rate it will see while the bandwidth is defined for the size of the channels. In the channelizer case, the highest rate is defined as the rate of the incoming signal, but in other PFB blocks, this is not so obvious.

Two very useful blocks to use are the arbitrary resampler and the clock synchronizer (for PAM signals). These PFBs are defined with a set number of filters based on the fidelity required from them, not the rate changes. By default, the filter_size is set to 32 for these blocks, which is a reasonable default for most tasks. Because the PFB uses this number of filters in the filterbank, the maximum rate of the bank is defined from this (see the theory of a polyphase interpolator for a justification of this). So the prototype filter is defined to use a sample rate of filter_size times the signal's sampling rate.

A helpful wrapper for the arbitrary resampler is found in gr-filter/python/pfb.py, which is exposed in Python as filter.pfb.arb_resampler_ccf and filter.pfb.arb_resampler_fff. This block is set up so that the user only needs to pass it the real number rate as the resampling rate. With just this information, this hierarchical block automatically creates a filter that fully passes the signal bandwidth being resampled but does not pass any out-of-band noise. See the code for this block for details of how the filter is constructed.

Of course, a user can create his or her own taps and use them in the arbitrary resampler for more specific requirements. Some of the UHD examples (gr-uhd/examples) use this ability to create a received matched filter or channel filter that also resamples the signal.

Examples

The following is an example of the using the channelizer. It creates the appropriate filter to channelizer 9 channels out of an original signal that is 9000 Hz wide, so each output channel is now 1000 Hz. The code then plots the PSD of the original signal to see the signals in the origina spectrum and then makes 9 plots for each of the channels.

NOTE: you need the Scipy and Matplotlib Python modules installed to run this example.

1 #!/usr/bin/env python
2 #
3 # Copyright 2009,2012,2013 Free Software Foundation, Inc.
4 #
5 # This file is part of GNU Radio
6 #
7 # GNU Radio is free software; you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation; either version 3, or (at your option)
10 # any later version.
11 #
12 # GNU Radio is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with GNU Radio; see the file COPYING. If not, write to
19 # the Free Software Foundation, Inc., 51 Franklin Street,
20 # Boston, MA 02110-1301, USA.
21 #
22 
23 from gnuradio import gr
24 from gnuradio import blocks
25 from gnuradio import filter
26 import sys, time
27 
28 try:
29  from gnuradio import analog
30 except ImportError:
31  sys.stderr.write("Error: Program requires gr-analog.\n")
32  sys.exit(1)
33 
34 try:
35  import scipy
36  from scipy import fftpack
37 except ImportError:
38  sys.stderr.write("Error: Program requires scipy (see: www.scipy.org).\n")
39  sys.exit(1)
40 
41 try:
42  import pylab
43  from pylab import mlab
44 except ImportError:
45  sys.stderr.write("Error: Program requires matplotlib (see: matplotlib.sourceforge.net).\n")
46  sys.exit(1)
47 
48 class pfb_top_block(gr.top_block):
49  def __init__(self):
50  gr.top_block.__init__(self)
51 
52  self._N = 2000000 # number of samples to use
53  self._fs = 1000 # initial sampling rate
54  self._M = M = 9 # Number of channels to channelize
55  self._ifs = M*self._fs # initial sampling rate
56 
57  # Create a set of taps for the PFB channelizer
58  self._taps = filter.firdes.low_pass_2(1, self._ifs, 475.50, 50,
59  attenuation_dB=100,
60  window=filter.firdes.WIN_BLACKMAN_hARRIS)
61 
62  # Calculate the number of taps per channel for our own information
63  tpc = scipy.ceil(float(len(self._taps)) / float(self._M))
64  print "Number of taps: ", len(self._taps)
65  print "Number of channels: ", self._M
66  print "Taps per channel: ", tpc
67 
68  # Create a set of signals at different frequencies
69  # freqs lists the frequencies of the signals that get stored
70  # in the list "signals", which then get summed together
71  self.signals = list()
72  self.add = blocks.add_cc()
73  freqs = [-70, -50, -30, -10, 10, 20, 40, 60, 80]
74  for i in xrange(len(freqs)):
75  f = freqs[i] + (M/2-M+i+1)*self._fs
76  self.signals.append(analog.sig_source_c(self._ifs, analog.GR_SIN_WAVE, f, 1))
77  self.connect(self.signals[i], (self.add,i))
78 
79  self.head = blocks.head(gr.sizeof_gr_complex, self._N)
80 
81  # Construct the channelizer filter
82  self.pfb = filter.pfb.channelizer_ccf(self._M, self._taps, 1)
83 
84  # Construct a vector sink for the input signal to the channelizer
85  self.snk_i = blocks.vector_sink_c()
86 
87  # Connect the blocks
88  self.connect(self.add, self.head, self.pfb)
89  self.connect(self.add, self.snk_i)
90 
91  # Use this to play with the channel mapping
92  #self.pfb.set_channel_map([5,6,7,8,0,1,2,3,4])
93 
94  # Create a vector sink for each of M output channels of the filter and connect it
95  self.snks = list()
96  for i in xrange(self._M):
97  self.snks.append(blocks.vector_sink_c())
98  self.connect((self.pfb, i), self.snks[i])
99 
100 
101 def main():
102  tstart = time.time()
103 
104  tb = pfb_top_block()
105  tb.run()
106 
107  tend = time.time()
108  print "Run time: %f" % (tend - tstart)
109 
110  if 1:
111  fig_in = pylab.figure(1, figsize=(16,9), facecolor="w")
112  fig1 = pylab.figure(2, figsize=(16,9), facecolor="w")
113  fig2 = pylab.figure(3, figsize=(16,9), facecolor="w")
114 
115  Ns = 1000
116  Ne = 10000
117 
118  fftlen = 8192
119  winfunc = scipy.blackman
120  fs = tb._ifs
121 
122  # Plot the input signal on its own figure
123  d = tb.snk_i.data()[Ns:Ne]
124  spin_f = fig_in.add_subplot(2, 1, 1)
125 
126  X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs,
127  window = lambda d: d*winfunc(fftlen),
128  scale_by_freq=True)
129  X_in = 10.0*scipy.log10(abs(X))
130  f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size))
131  pin_f = spin_f.plot(f_in, X_in, "b")
132  spin_f.set_xlim([min(f_in), max(f_in)+1])
133  spin_f.set_ylim([-200.0, 50.0])
134 
135  spin_f.set_title("Input Signal", weight="bold")
136  spin_f.set_xlabel("Frequency (Hz)")
137  spin_f.set_ylabel("Power (dBW)")
138 
139 
140  Ts = 1.0/fs
141  Tmax = len(d)*Ts
142 
143  t_in = scipy.arange(0, Tmax, Ts)
144  x_in = scipy.array(d)
145  spin_t = fig_in.add_subplot(2, 1, 2)
146  pin_t = spin_t.plot(t_in, x_in.real, "b")
147  pin_t = spin_t.plot(t_in, x_in.imag, "r")
148 
149  spin_t.set_xlabel("Time (s)")
150  spin_t.set_ylabel("Amplitude")
151 
152  Ncols = int(scipy.floor(scipy.sqrt(tb._M)))
153  Nrows = int(scipy.floor(tb._M / Ncols))
154  if(tb._M % Ncols != 0):
155  Nrows += 1
156 
157  # Plot each of the channels outputs. Frequencies on Figure 2 and
158  # time signals on Figure 3
159  fs_o = tb._fs
160  Ts_o = 1.0/fs_o
161  Tmax_o = len(d)*Ts_o
162  for i in xrange(len(tb.snks)):
163  # remove issues with the transients at the beginning
164  # also remove some corruption at the end of the stream
165  # this is a bug, probably due to the corner cases
166  d = tb.snks[i].data()[Ns:Ne]
167 
168  sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i)
169  X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o,
170  window = lambda d: d*winfunc(fftlen),
171  scale_by_freq=True)
172  X_o = 10.0*scipy.log10(abs(X))
173  f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size))
174  p2_f = sp1_f.plot(f_o, X_o, "b")
175  sp1_f.set_xlim([min(f_o), max(f_o)+1])
176  sp1_f.set_ylim([-200.0, 50.0])
177 
178  sp1_f.set_title(("Channel %d" % i), weight="bold")
179  sp1_f.set_xlabel("Frequency (Hz)")
180  sp1_f.set_ylabel("Power (dBW)")
181 
182  x_o = scipy.array(d)
183  t_o = scipy.arange(0, Tmax_o, Ts_o)
184  sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i)
185  p2_o = sp2_o.plot(t_o, x_o.real, "b")
186  p2_o = sp2_o.plot(t_o, x_o.imag, "r")
187  sp2_o.set_xlim([min(t_o), max(t_o)+1])
188  sp2_o.set_ylim([-2, 2])
189 
190  sp2_o.set_title(("Channel %d" % i), weight="bold")
191  sp2_o.set_xlabel("Time (s)")
192  sp2_o.set_ylabel("Amplitude")
193 
194  pylab.show()
195 
196 
197 if __name__ == "__main__":
198  try:
199  main()
200  except KeyboardInterrupt:
201  pass
202 

The PFB Arbitrary Resampler Kernel

GNU Radio has a PFB arbitrary resampler block that can be used to resample a signal to any arbitrary and real resampling rate. The resampling feature is one that could easily be useful to other blocks, and so we have extracted the kernel of the resampler into its own class that can be used as such.

The PFB arbitrary resampler is defined in pfb_arb_resampler.h and has the following constructor:

namespace gr {
namespace filter {
namespace kernel {
pfb_arb_resampler_XXX(float rate,
const std::vector<float> &taps,
unsigned int filter_size);
} /* namespace kernel */
} /* namespace filter */
} /* namespace gr */

Currently, only a 'ccf' and 'fff' version are defined. This kernel, like the block itself, takes in the resampling rate as a floating point number. The taps are passed as the baseband prototype filter, and the quantization error of the filter is determined by the filter_size parameter.

The prototype taps are generated like all other PFB filter taps. Specifically, we construct them generally as a lowpass filter at the maximum rate of the filter. In the case of these resamplers, the maximum rate is actually the number of filters.

A simple example follows. We construct a filter that will pass the entire passband of the original signal to be resampled. To make it easy, we work in normalized sample rates for this. The gain of the filter is set to filter_size to compensate for the upsampling, the sampling rate itself is also set to filter_size, which is assuming that the incoming signal is at a sampling rate of 1.0. We defined the passband to be 0.5 to pass the entire width of the original signal and set a transition band to 0.1. Note that this causes a bit of roll-off outside of the original passband and could lead to introducing some aliasing. More care should be taken to construct the passband and transition width of the filter for the given signal while keeping the total number of taps small. A stopband attenuation of 60 dB was used here, and again, this is a parameter we can adjust to alter the performance and size of the filter.

firdes.low_pass_2(filter_size, filter_size, 0.5, 0.1, 60)

As is typical with the PFB filters, a filter size of 32 is generally an appropriate trade-off of accuracy, performance, and memory. This should provide an error roughly equivalent to the quanization error of using 16-bit fixed point representation. Generally, increasing over 32 provides some accuracy benefits without a huge increase in computational demands.