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 23from __future__ import print_function 24from __future__ import division 25from __future__ import unicode_literals 26from gnuradio import gr 27from gnuradio import blocks 28from gnuradio import filter 29import sys, time 30import numpy 31 32try: 33 from gnuradio import analog 34except ImportError: 35 sys.stderr.write("Error: Program requires gr-analog.\n") 36 sys.exit(1) 37 38 39try: 40 from matplotlib import pyplot 41 from matplotlib import pyplot as mlab 42except ImportError: 43 sys.stderr.write("Error: Program requires matplotlib (see: matplotlib.sourceforge.net).\n") 44 sys.exit(1) 45 46class pfb_top_block(gr.top_block): 47 def __init__(self): 48 gr.top_block.__init__(self) 49 50 self._N = 10000000 # number of samples to use 51 self._fs = 10000 # initial sampling rate 52 self._decim = 20 # Decimation rate 53 54 # Generate the prototype filter taps for the decimators with a 200 Hz bandwidth 55 self._taps = filter.firdes.low_pass_2(1, self._fs, 56 200, 150, 57 attenuation_dB=120, 58 window=filter.firdes.WIN_BLACKMAN_hARRIS) 59 60 # Calculate the number of taps per channel for our own information 61 tpc = numpy.ceil(float(len(self._taps)) / float(self._decim)) 62 print("Number of taps: ", len(self._taps)) 63 print("Number of filters: ", self._decim) 64 print("Taps per channel: ", tpc) 65 66 # Build the input signal source 67 # We create a list of freqs, and a sine wave is generated and added to the source 68 # for each one of these frequencies. 69 self.signals = list() 70 self.add = blocks.add_cc() 71 freqs = [10, 20, 2040] 72 for i in range(len(freqs)): 73 self.signals.append(analog.sig_source_c(self._fs, analog.GR_SIN_WAVE, freqs[i], 1)) 74 self.connect(self.signals[i], (self.add,i)) 75 76 self.head = blocks.head(gr.sizeof_gr_complex, self._N) 77 78 # Construct a PFB decimator filter 79 self.pfb = filter.pfb.decimator_ccf(self._decim, self._taps, 0) 80 81 # Construct a standard FIR decimating filter 82 self.dec = filter.fir_filter_ccf(self._decim, self._taps) 83 84 self.snk_i = blocks.vector_sink_c() 85 86 # Connect the blocks 87 self.connect(self.add, self.head, self.pfb) 88 self.connect(self.add, self.snk_i) 89 90 # Create the sink for the decimated siganl 91 self.snk = blocks.vector_sink_c() 92 self.connect(self.pfb, self.snk) 93 94 95def main(): 96 tb = pfb_top_block() 97 98 tstart = time.time() 99 tb.run() 100 tend = time.time() 101 print("Run time: %f" % (tend - tstart)) 102 103 if 1: 104 fig1 = pyplot.figure(1, figsize=(16,9)) 105 fig2 = pyplot.figure(2, figsize=(16,9)) 106 107 Ns = 10000 108 Ne = 10000 109 110 fftlen = 8192 111 winfunc = numpy.blackman 112 fs = tb._fs 113 114 # Plot the input to the decimator 115 116 d = tb.snk_i.data()[Ns:Ns+Ne] 117 sp1_f = fig1.add_subplot(2, 1, 1) 118 119 X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen / 4, Fs=fs, 120 window = lambda d: d*winfunc(fftlen), 121 scale_by_freq=True) 122 X_in = 10.0*numpy.log10(abs(numpy.fft.fftshift(X))) 123 f_in = numpy.arange(-fs / 2.0, fs / 2.0, fs / float(X_in.size)) 124 p1_f = sp1_f.plot(f_in, X_in, "b") 125 sp1_f.set_xlim([min(f_in), max(f_in)+1]) 126 sp1_f.set_ylim([-200.0, 50.0]) 127 128 sp1_f.set_title("Input Signal", weight="bold") 129 sp1_f.set_xlabel("Frequency (Hz)") 130 sp1_f.set_ylabel("Power (dBW)") 131 132 Ts = 1.0 / fs 133 Tmax = len(d)*Ts 134 135 t_in = numpy.arange(0, Tmax, Ts) 136 x_in = numpy.array(d) 137 sp1_t = fig1.add_subplot(2, 1, 2) 138 p1_t = sp1_t.plot(t_in, x_in.real, "b") 139 p1_t = sp1_t.plot(t_in, x_in.imag, "r") 140 sp1_t.set_ylim([-tb._decim*1.1, tb._decim*1.1]) 141 142 sp1_t.set_xlabel("Time (s)") 143 sp1_t.set_ylabel("Amplitude") 144 145 146 # Plot the output of the decimator 147 fs_o = tb._fs / tb._decim 148 149 sp2_f = fig2.add_subplot(2, 1, 1) 150 d = tb.snk.data()[Ns:Ns+Ne] 151 X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen / 4, Fs=fs_o, 152 window = lambda d: d*winfunc(fftlen), 153 scale_by_freq=True) 154 X_o = 10.0*numpy.log10(abs(numpy.fft.fftshift(X))) 155 f_o = numpy.arange(-fs_o / 2.0, fs_o / 2.0, fs_o / float(X_o.size)) 156 p2_f = sp2_f.plot(f_o, X_o, "b") 157 sp2_f.set_xlim([min(f_o), max(f_o)+1]) 158 sp2_f.set_ylim([-200.0, 50.0]) 159 160 sp2_f.set_title("PFB Decimated Signal", weight="bold") 161 sp2_f.set_xlabel("Frequency (Hz)") 162 sp2_f.set_ylabel("Power (dBW)") 163 164 165 Ts_o = 1.0 / fs_o 166 Tmax_o = len(d)*Ts_o 167 168 x_o = numpy.array(d) 169 t_o = numpy.arange(0, Tmax_o, Ts_o) 170 sp2_t = fig2.add_subplot(2, 1, 2) 171 p2_t = sp2_t.plot(t_o, x_o.real, "b-o") 172 p2_t = sp2_t.plot(t_o, x_o.imag, "r-o") 173 sp2_t.set_ylim([-2.5, 2.5]) 174 175 sp2_t.set_xlabel("Time (s)") 176 sp2_t.set_ylabel("Amplitude") 177 178 pyplot.show() 179 180 181if __name__ == "__main__": 182 try: 183 main() 184 except KeyboardInterrupt: 185 pass 186 187