GNU Radio 3.6.4.1 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. These blocks and their documentation can be found in Polyphase Filterbank.

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 gnuradio-examples/python/pfb.

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 gnuradio-examples/python/pfb/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 = gr.firdes.low_pass_2(1, self._fs, 475.50, 50,
attenuation_dB=100,
window=gr.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 gnuradio-core/src/python/gnuradio/blks2impl/pfb_arb_resampler.py, which is exposed in Python as blks2.pfb_arb_resampler_ccf and blks2.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 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, blks2
24 import sys, time
25 
26 try:
27  import scipy
28  from scipy import fftpack
29 except ImportError:
30  print "Error: Program requires scipy (see: www.scipy.org)."
31  sys.exit(1)
32 
33 try:
34  import pylab
35  from pylab import mlab
36 except ImportError:
37  print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)."
38  sys.exit(1)
39 
40 class pfb_top_block(gr.top_block):
41  def __init__(self):
42  gr.top_block.__init__(self)
43 
44  self._N = 2000000 # number of samples to use
45  self._fs = 9000 # initial sampling rate
46  self._M = 9 # Number of channels to channelize
47 
48  # Create a set of taps for the PFB channelizer
49  self._taps = gr.firdes.low_pass_2(1, self._fs, 475.50, 50,
50  attenuation_dB=100, window=gr.firdes.WIN_BLACKMAN_hARRIS)
51 
52  # Calculate the number of taps per channel for our own information
53  tpc = scipy.ceil(float(len(self._taps)) / float(self._M))
54  print "Number of taps: ", len(self._taps)
55  print "Number of channels: ", self._M
56  print "Taps per channel: ", tpc
57 
58  # Create a set of signals at different frequencies
59  # freqs lists the frequencies of the signals that get stored
60  # in the list "signals", which then get summed together
61  self.signals = list()
62  self.add = gr.add_cc()
63  freqs = [-4070, -3050, -2030, -1010, 10, 1020, 2040, 3060, 4080]
64  for i in xrange(len(freqs)):
65  self.signals.append(gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freqs[i], 1))
66  self.connect(self.signals[i], (self.add,i))
67 
68  self.head = gr.head(gr.sizeof_gr_complex, self._N)
69 
70  # Construct the channelizer filter
71  self.pfb = blks2.pfb_channelizer_ccf(self._M, self._taps, 1)
72 
73  # Construct a vector sink for the input signal to the channelizer
74  self.snk_i = gr.vector_sink_c()
75 
76  # Connect the blocks
77  self.connect(self.add, self.head, self.pfb)
78  self.connect(self.add, self.snk_i)
79 
80  # Use this to play with the channel mapping
81  #self.pfb.set_channel_map([5,6,7,8,0,1,2,3,4])
82 
83  # Create a vector sink for each of M output channels of the filter and connect it
84  self.snks = list()
85  for i in xrange(self._M):
86  self.snks.append(gr.vector_sink_c())
87  self.connect((self.pfb, i), self.snks[i])
88 
89 
90 def main():
91  tstart = time.time()
92 
93  tb = pfb_top_block()
94  tb.run()
95 
96  tend = time.time()
97  print "Run time: %f" % (tend - tstart)
98 
99  if 1:
100  fig_in = pylab.figure(1, figsize=(16,9), facecolor="w")
101  fig1 = pylab.figure(2, figsize=(16,9), facecolor="w")
102  fig2 = pylab.figure(3, figsize=(16,9), facecolor="w")
103 
104  Ns = 1000
105  Ne = 10000
106 
107  fftlen = 8192
108  winfunc = scipy.blackman
109  fs = tb._fs
110 
111  # Plot the input signal on its own figure
112  d = tb.snk_i.data()[Ns:Ne]
113  spin_f = fig_in.add_subplot(2, 1, 1)
114 
115  X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs,
116  window = lambda d: d*winfunc(fftlen),
117  scale_by_freq=True)
118  X_in = 10.0*scipy.log10(abs(X))
119  f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size))
120  pin_f = spin_f.plot(f_in, X_in, "b")
121  spin_f.set_xlim([min(f_in), max(f_in)+1])
122  spin_f.set_ylim([-200.0, 50.0])
123 
124  spin_f.set_title("Input Signal", weight="bold")
125  spin_f.set_xlabel("Frequency (Hz)")
126  spin_f.set_ylabel("Power (dBW)")
127 
128 
129  Ts = 1.0/fs
130  Tmax = len(d)*Ts
131 
132  t_in = scipy.arange(0, Tmax, Ts)
133  x_in = scipy.array(d)
134  spin_t = fig_in.add_subplot(2, 1, 2)
135  pin_t = spin_t.plot(t_in, x_in.real, "b")
136  pin_t = spin_t.plot(t_in, x_in.imag, "r")
137 
138  spin_t.set_xlabel("Time (s)")
139  spin_t.set_ylabel("Amplitude")
140 
141  Ncols = int(scipy.floor(scipy.sqrt(tb._M)))
142  Nrows = int(scipy.floor(tb._M / Ncols))
143  if(tb._M % Ncols != 0):
144  Nrows += 1
145 
146  # Plot each of the channels outputs. Frequencies on Figure 2 and
147  # time signals on Figure 3
148  fs_o = tb._fs / tb._M
149  Ts_o = 1.0/fs_o
150  Tmax_o = len(d)*Ts_o
151  for i in xrange(len(tb.snks)):
152  # remove issues with the transients at the beginning
153  # also remove some corruption at the end of the stream
154  # this is a bug, probably due to the corner cases
155  d = tb.snks[i].data()[Ns:Ne]
156 
157  sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i)
158  X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o,
159  window = lambda d: d*winfunc(fftlen),
160  scale_by_freq=True)
161  X_o = 10.0*scipy.log10(abs(X))
162  f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size))
163  p2_f = sp1_f.plot(f_o, X_o, "b")
164  sp1_f.set_xlim([min(f_o), max(f_o)+1])
165  sp1_f.set_ylim([-200.0, 50.0])
166 
167  sp1_f.set_title(("Channel %d" % i), weight="bold")
168  sp1_f.set_xlabel("Frequency (Hz)")
169  sp1_f.set_ylabel("Power (dBW)")
170 
171  x_o = scipy.array(d)
172  t_o = scipy.arange(0, Tmax_o, Ts_o)
173  sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i)
174  p2_o = sp2_o.plot(t_o, x_o.real, "b")
175  p2_o = sp2_o.plot(t_o, x_o.imag, "r")
176  sp2_o.set_xlim([min(t_o), max(t_o)+1])
177  sp2_o.set_ylim([-2, 2])
178 
179  sp2_o.set_title(("Channel %d" % i), weight="bold")
180  sp2_o.set_xlabel("Time (s)")
181  sp2_o.set_ylabel("Amplitude")
182 
183  pylab.show()
184 
185 
186 if __name__ == "__main__":
187  try:
188  main()
189  except KeyboardInterrupt:
190  pass
191