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