git.gag.com Git - debian/gnuradio/blob - gnuradio-examples/python/pfb/decimate.py

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