From faab807cf5e8b4a4c950d1fd8ae6672296fb1ad9 Mon Sep 17 00:00:00 2001 From: Tom Rondeau Date: Tue, 3 Apr 2012 18:21:52 -0400 Subject: Rework example directories. Gets rid of gnuradio-examples, moves these to more appropriate components. gnuradio-core and grc now have their own examples directories for files directly related to them. --- gnuradio-core/src/examples/pfb/CMakeLists.txt | 36 ++ gnuradio-core/src/examples/pfb/channelize.py | 191 +++++++ gnuradio-core/src/examples/pfb/chirp_channelize.py | 203 +++++++ gnuradio-core/src/examples/pfb/decimate.py | 178 ++++++ gnuradio-core/src/examples/pfb/fmtest.py | 225 ++++++++ gnuradio-core/src/examples/pfb/interpolate.py | 233 ++++++++ gnuradio-core/src/examples/pfb/reconstruction.py | 131 +++++ gnuradio-core/src/examples/pfb/resampler.py | 127 +++++ gnuradio-core/src/examples/pfb/resampler_demo.grc | 598 +++++++++++++++++++++ gnuradio-core/src/examples/pfb/synth_filter.py | 83 +++ gnuradio-core/src/examples/pfb/synth_to_chan.py | 117 ++++ 11 files changed, 2122 insertions(+) create mode 100644 gnuradio-core/src/examples/pfb/CMakeLists.txt create mode 100755 gnuradio-core/src/examples/pfb/channelize.py create mode 100755 gnuradio-core/src/examples/pfb/chirp_channelize.py create mode 100755 gnuradio-core/src/examples/pfb/decimate.py create mode 100755 gnuradio-core/src/examples/pfb/fmtest.py create mode 100755 gnuradio-core/src/examples/pfb/interpolate.py create mode 100755 gnuradio-core/src/examples/pfb/reconstruction.py create mode 100755 gnuradio-core/src/examples/pfb/resampler.py create mode 100644 gnuradio-core/src/examples/pfb/resampler_demo.grc create mode 100755 gnuradio-core/src/examples/pfb/synth_filter.py create mode 100755 gnuradio-core/src/examples/pfb/synth_to_chan.py (limited to 'gnuradio-core/src/examples/pfb') diff --git a/gnuradio-core/src/examples/pfb/CMakeLists.txt b/gnuradio-core/src/examples/pfb/CMakeLists.txt new file mode 100644 index 000000000..899c47b90 --- /dev/null +++ b/gnuradio-core/src/examples/pfb/CMakeLists.txt @@ -0,0 +1,36 @@ +# Copyright 2011 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. + +include(GrPython) + +GR_PYTHON_INSTALL(PROGRAMS + channelize.py + chirp_channelize.py + decimate.py + fmtest.py + interpolate.py + resampler_demo.grc + resampler.py + synth_filter.py + synth_to_chan.py + reconstruction.py + DESTINATION ${GR_PKG_DATA_DIR}/examples/pfb + COMPONENT "core_python" +) + diff --git a/gnuradio-core/src/examples/pfb/channelize.py b/gnuradio-core/src/examples/pfb/channelize.py new file mode 100755 index 000000000..2fcb14a36 --- /dev/null +++ b/gnuradio-core/src/examples/pfb/channelize.py @@ -0,0 +1,191 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys, time + +try: + import scipy + from scipy import fftpack +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab + from pylab import mlab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 2000000 # number of samples to use + self._fs = 9000 # initial sampling rate + self._M = 9 # Number of channels to channelize + + # Create a set of taps for the PFB channelizer + self._taps = gr.firdes.low_pass_2(1, self._fs, 475.50, 50, + attenuation_dB=100, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._M)) + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + # Create a set of signals at different frequencies + # freqs lists the frequencies of the signals that get stored + # in the list "signals", which then get summed together + self.signals = list() + self.add = gr.add_cc() + freqs = [-4070, -3050, -2030, -1010, 10, 1020, 2040, 3060, 4080] + for i in xrange(len(freqs)): + self.signals.append(gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freqs[i], 1)) + self.connect(self.signals[i], (self.add,i)) + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the channelizer filter + self.pfb = blks2.pfb_channelizer_ccf(self._M, self._taps, 1) + + # Construct a vector sink for the input signal to the channelizer + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.add, self.head, self.pfb) + self.connect(self.add, self.snk_i) + + # Use this to play with the channel mapping + #self.pfb.set_channel_map([5,6,7,8,0,1,2,3,4]) + + # Create a vector sink for each of M output channels of the filter and connect it + self.snks = list() + for i in xrange(self._M): + self.snks.append(gr.vector_sink_c()) + self.connect((self.pfb, i), self.snks[i]) + + +def main(): + tstart = time.time() + + tb = pfb_top_block() + tb.run() + + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig_in = pylab.figure(1, figsize=(16,9), facecolor="w") + fig1 = pylab.figure(2, figsize=(16,9), facecolor="w") + fig2 = pylab.figure(3, figsize=(16,9), facecolor="w") + + Ns = 1000 + Ne = 10000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._fs + + # Plot the input signal on its own figure + d = tb.snk_i.data()[Ns:Ne] + spin_f = fig_in.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(X)) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + pin_f = spin_f.plot(f_in, X_in, "b") + spin_f.set_xlim([min(f_in), max(f_in)+1]) + spin_f.set_ylim([-200.0, 50.0]) + + spin_f.set_title("Input Signal", weight="bold") + spin_f.set_xlabel("Frequency (Hz)") + spin_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + spin_t = fig_in.add_subplot(2, 1, 2) + pin_t = spin_t.plot(t_in, x_in.real, "b") + pin_t = spin_t.plot(t_in, x_in.imag, "r") + + spin_t.set_xlabel("Time (s)") + spin_t.set_ylabel("Amplitude") + + Ncols = int(scipy.floor(scipy.sqrt(tb._M))) + Nrows = int(scipy.floor(tb._M / Ncols)) + if(tb._M % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = tb._fs / tb._M + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + for i in xrange(len(tb.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = tb.snks[i].data()[Ns:Ne] + + sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i) + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(X)) + f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size)) + p2_f = sp1_f.plot(f_o, X_o, "b") + sp1_f.set_xlim([min(f_o), max(f_o)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title(("Channel %d" % i), weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i) + p2_o = sp2_o.plot(t_o, x_o.real, "b") + p2_o = sp2_o.plot(t_o, x_o.imag, "r") + sp2_o.set_xlim([min(t_o), max(t_o)+1]) + sp2_o.set_ylim([-2, 2]) + + sp2_o.set_title(("Channel %d" % i), weight="bold") + sp2_o.set_xlabel("Time (s)") + sp2_o.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-core/src/examples/pfb/chirp_channelize.py b/gnuradio-core/src/examples/pfb/chirp_channelize.py new file mode 100755 index 000000000..951255d3b --- /dev/null +++ b/gnuradio-core/src/examples/pfb/chirp_channelize.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys, time + +try: + import scipy + from scipy import fftpack +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab + from pylab import mlab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 200000 # number of samples to use + self._fs = 9000 # initial sampling rate + self._M = 9 # Number of channels to channelize + + # Create a set of taps for the PFB channelizer + self._taps = gr.firdes.low_pass_2(1, self._fs, 500, 20, + attenuation_dB=10, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._M)) + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + repeated = True + if(repeated): + self.vco_input = gr.sig_source_f(self._fs, gr.GR_SIN_WAVE, 0.25, 110) + else: + amp = 100 + data = scipy.arange(0, amp, amp/float(self._N)) + self.vco_input = gr.vector_source_f(data, False) + + # Build a VCO controlled by either the sinusoid or single chirp tone + # Then convert this to a complex signal + self.vco = gr.vco_f(self._fs, 225, 1) + self.f2c = gr.float_to_complex() + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the channelizer filter + self.pfb = blks2.pfb_channelizer_ccf(self._M, self._taps) + + # Construct a vector sink for the input signal to the channelizer + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.vco_input, self.vco, self.f2c) + self.connect(self.f2c, self.head, self.pfb) + self.connect(self.f2c, self.snk_i) + + # Create a vector sink for each of M output channels of the filter and connect it + self.snks = list() + for i in xrange(self._M): + self.snks.append(gr.vector_sink_c()) + self.connect((self.pfb, i), self.snks[i]) + + +def main(): + tstart = time.time() + + tb = pfb_top_block() + tb.run() + + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig_in = pylab.figure(1, figsize=(16,9), facecolor="w") + fig1 = pylab.figure(2, figsize=(16,9), facecolor="w") + fig2 = pylab.figure(3, figsize=(16,9), facecolor="w") + fig3 = pylab.figure(4, figsize=(16,9), facecolor="w") + + Ns = 650 + Ne = 20000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._fs + + # Plot the input signal on its own figure + d = tb.snk_i.data()[Ns:Ne] + spin_f = fig_in.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + pin_f = spin_f.plot(f_in, X_in, "b") + spin_f.set_xlim([min(f_in), max(f_in)+1]) + spin_f.set_ylim([-200.0, 50.0]) + + spin_f.set_title("Input Signal", weight="bold") + spin_f.set_xlabel("Frequency (Hz)") + spin_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + spin_t = fig_in.add_subplot(2, 1, 2) + pin_t = spin_t.plot(t_in, x_in.real, "b") + pin_t = spin_t.plot(t_in, x_in.imag, "r") + + spin_t.set_xlabel("Time (s)") + spin_t.set_ylabel("Amplitude") + + Ncols = int(scipy.floor(scipy.sqrt(tb._M))) + Nrows = int(scipy.floor(tb._M / Ncols)) + if(tb._M % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = tb._fs / tb._M + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + for i in xrange(len(tb.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = tb.snks[i].data()[Ns:Ne] + + sp1_f = fig1.add_subplot(Nrows, Ncols, 1+i) + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(X)) + f_o = freq + p2_f = sp1_f.plot(f_o, X_o, "b") + sp1_f.set_xlim([min(f_o), max(f_o)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title(("Channel %d" % i), weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_o = fig2.add_subplot(Nrows, Ncols, 1+i) + p2_o = sp2_o.plot(t_o, x_o.real, "b") + p2_o = sp2_o.plot(t_o, x_o.imag, "r") + sp2_o.set_xlim([min(t_o), max(t_o)+1]) + sp2_o.set_ylim([-2, 2]) + + sp2_o.set_title(("Channel %d" % i), weight="bold") + sp2_o.set_xlabel("Time (s)") + sp2_o.set_ylabel("Amplitude") + + + sp3 = fig3.add_subplot(1,1,1) + p3 = sp3.plot(t_o, x_o.real) + sp3.set_xlim([min(t_o), max(t_o)+1]) + sp3.set_ylim([-2, 2]) + + sp3.set_title("All Channels") + sp3.set_xlabel("Time (s)") + sp3.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-core/src/examples/pfb/decimate.py b/gnuradio-core/src/examples/pfb/decimate.py new file mode 100755 index 000000000..643a2c241 --- /dev/null +++ b/gnuradio-core/src/examples/pfb/decimate.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys, time + +try: + import scipy + from scipy import fftpack +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab + from pylab import mlab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 10000000 # number of samples to use + self._fs = 10000 # initial sampling rate + self._decim = 20 # Decimation rate + + # Generate the prototype filter taps for the decimators with a 200 Hz bandwidth + self._taps = gr.firdes.low_pass_2(1, self._fs, 200, 150, + attenuation_dB=120, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._decim)) + print "Number of taps: ", len(self._taps) + print "Number of filters: ", self._decim + print "Taps per channel: ", tpc + + # Build the input signal source + # We create a list of freqs, and a sine wave is generated and added to the source + # for each one of these frequencies. + self.signals = list() + self.add = gr.add_cc() + freqs = [10, 20, 2040] + for i in xrange(len(freqs)): + self.signals.append(gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freqs[i], 1)) + self.connect(self.signals[i], (self.add,i)) + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct a PFB decimator filter + self.pfb = blks2.pfb_decimator_ccf(self._decim, self._taps, 0) + + # Construct a standard FIR decimating filter + self.dec = gr.fir_filter_ccf(self._decim, self._taps) + + self.snk_i = gr.vector_sink_c() + + # Connect the blocks + self.connect(self.add, self.head, self.pfb) + self.connect(self.add, self.snk_i) + + # Create the sink for the decimated siganl + self.snk = gr.vector_sink_c() + self.connect(self.pfb, self.snk) + + +def main(): + tb = pfb_top_block() + + tstart = time.time() + tb.run() + tend = time.time() + print "Run time: %f" % (tend - tstart) + + if 1: + fig1 = pylab.figure(1, figsize=(16,9)) + fig2 = pylab.figure(2, figsize=(16,9)) + + Ns = 10000 + Ne = 10000 + + fftlen = 8192 + winfunc = scipy.blackman + fs = tb._fs + + # Plot the input to the decimator + + d = tb.snk_i.data()[Ns:Ns+Ne] + sp1_f = fig1.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + p1_f = sp1_f.plot(f_in, X_in, "b") + sp1_f.set_xlim([min(f_in), max(f_in)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + sp1_f.set_title("Input Signal", weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + sp1_t = fig1.add_subplot(2, 1, 2) + p1_t = sp1_t.plot(t_in, x_in.real, "b") + p1_t = sp1_t.plot(t_in, x_in.imag, "r") + sp1_t.set_ylim([-tb._decim*1.1, tb._decim*1.1]) + + sp1_t.set_xlabel("Time (s)") + sp1_t.set_ylabel("Amplitude") + + + # Plot the output of the decimator + fs_o = tb._fs / tb._decim + + sp2_f = fig2.add_subplot(2, 1, 1) + d = tb.snk.data()[Ns:Ns+Ne] + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size)) + p2_f = sp2_f.plot(f_o, X_o, "b") + sp2_f.set_xlim([min(f_o), max(f_o)+1]) + sp2_f.set_ylim([-200.0, 50.0]) + + sp2_f.set_title("PFB Decimated Signal", weight="bold") + sp2_f.set_xlabel("Frequency (Hz)") + sp2_f.set_ylabel("Power (dBW)") + + + Ts_o = 1.0/fs_o + Tmax_o = len(d)*Ts_o + + x_o = scipy.array(d) + t_o = scipy.arange(0, Tmax_o, Ts_o) + sp2_t = fig2.add_subplot(2, 1, 2) + p2_t = sp2_t.plot(t_o, x_o.real, "b-o") + p2_t = sp2_t.plot(t_o, x_o.imag, "r-o") + sp2_t.set_ylim([-2.5, 2.5]) + + sp2_t.set_xlabel("Time (s)") + sp2_t.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-core/src/examples/pfb/fmtest.py b/gnuradio-core/src/examples/pfb/fmtest.py new file mode 100755 index 000000000..635ee4e9e --- /dev/null +++ b/gnuradio-core/src/examples/pfb/fmtest.py @@ -0,0 +1,225 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys, math, time + +try: + import scipy + from scipy import fftpack +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + + +class fmtx(gr.hier_block2): + def __init__(self, lo_freq, audio_rate, if_rate): + + gr.hier_block2.__init__(self, "build_fm", + gr.io_signature(1, 1, gr.sizeof_float), # Input signature + gr.io_signature(1, 1, gr.sizeof_gr_complex)) # Output signature + + fmtx = blks2.nbfm_tx (audio_rate, if_rate, max_dev=5e3, tau=75e-6) + + # Local oscillator + lo = gr.sig_source_c (if_rate, # sample rate + gr.GR_SIN_WAVE, # waveform type + lo_freq, #frequency + 1.0, # amplitude + 0) # DC Offset + mixer = gr.multiply_cc () + + self.connect (self, fmtx, (mixer, 0)) + self.connect (lo, (mixer, 1)) + self.connect (mixer, self) + +class fmtest(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._nsamples = 1000000 + self._audio_rate = 8000 + + # Set up N channels with their own baseband and IF frequencies + self._N = 5 + chspacing = 16000 + freq = [10, 20, 30, 40, 50] + f_lo = [0, 1*chspacing, -1*chspacing, 2*chspacing, -2*chspacing] + + self._if_rate = 4*self._N*self._audio_rate + + # Create a signal source and frequency modulate it + self.sum = gr.add_cc () + for n in xrange(self._N): + sig = gr.sig_source_f(self._audio_rate, gr.GR_SIN_WAVE, freq[n], 0.5) + fm = fmtx(f_lo[n], self._audio_rate, self._if_rate) + self.connect(sig, fm) + self.connect(fm, (self.sum, n)) + + self.head = gr.head(gr.sizeof_gr_complex, self._nsamples) + self.snk_tx = gr.vector_sink_c() + self.channel = blks2.channel_model(0.1) + + self.connect(self.sum, self.head, self.channel, self.snk_tx) + + + # Design the channlizer + self._M = 10 + bw = chspacing/2.0 + t_bw = chspacing/10.0 + self._chan_rate = self._if_rate / self._M + self._taps = gr.firdes.low_pass_2(1, self._if_rate, bw, t_bw, + attenuation_dB=100, + window=gr.firdes.WIN_BLACKMAN_hARRIS) + tpc = math.ceil(float(len(self._taps)) / float(self._M)) + + print "Number of taps: ", len(self._taps) + print "Number of channels: ", self._M + print "Taps per channel: ", tpc + + self.pfb = blks2.pfb_channelizer_ccf(self._M, self._taps) + + self.connect(self.channel, self.pfb) + + # Create a file sink for each of M output channels of the filter and connect it + self.fmdet = list() + self.squelch = list() + self.snks = list() + for i in xrange(self._M): + self.fmdet.append(blks2.nbfm_rx(self._audio_rate, self._chan_rate)) + self.squelch.append(blks2.standard_squelch(self._audio_rate*10)) + self.snks.append(gr.vector_sink_f()) + self.connect((self.pfb, i), self.fmdet[i], self.squelch[i], self.snks[i]) + + def num_tx_channels(self): + return self._N + + def num_rx_channels(self): + return self._M + +def main(): + + fm = fmtest() + + tstart = time.time() + fm.run() + tend = time.time() + + if 1: + fig1 = pylab.figure(1, figsize=(12,10), facecolor="w") + fig2 = pylab.figure(2, figsize=(12,10), facecolor="w") + fig3 = pylab.figure(3, figsize=(12,10), facecolor="w") + + Ns = 10000 + Ne = 100000 + + fftlen = 8192 + winfunc = scipy.blackman + + # Plot transmitted signal + fs = fm._if_rate + + d = fm.snk_tx.data()[Ns:Ns+Ne] + sp1_f = fig1.add_subplot(2, 1, 1) + + X,freq = sp1_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + visible=False) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + p1_f = sp1_f.plot(f_in, X_in, "b") + sp1_f.set_xlim([min(f_in), max(f_in)+1]) + sp1_f.set_ylim([-120.0, 20.0]) + + sp1_f.set_title("Input Signal", weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + sp1_t = fig1.add_subplot(2, 1, 2) + p1_t = sp1_t.plot(t_in, x_in.real, "b-o") + #p1_t = sp1_t.plot(t_in, x_in.imag, "r-o") + sp1_t.set_ylim([-5, 5]) + + # Set up the number of rows and columns for plotting the subfigures + Ncols = int(scipy.floor(scipy.sqrt(fm.num_rx_channels()))) + Nrows = int(scipy.floor(fm.num_rx_channels() / Ncols)) + if(fm.num_rx_channels() % Ncols != 0): + Nrows += 1 + + # Plot each of the channels outputs. Frequencies on Figure 2 and + # time signals on Figure 3 + fs_o = fm._audio_rate + for i in xrange(len(fm.snks)): + # remove issues with the transients at the beginning + # also remove some corruption at the end of the stream + # this is a bug, probably due to the corner cases + d = fm.snks[i].data()[Ns:Ne] + + sp2_f = fig2.add_subplot(Nrows, Ncols, 1+i) + X,freq = sp2_f.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs_o, + window = lambda d: d*winfunc(fftlen), + visible=False) + #X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + X_o = 10.0*scipy.log10(abs(X)) + #f_o = scipy.arange(-fs_o/2.0, fs_o/2.0, fs_o/float(X_o.size)) + f_o = scipy.arange(0, fs_o/2.0, fs_o/2.0/float(X_o.size)) + p2_f = sp2_f.plot(f_o, X_o, "b") + sp2_f.set_xlim([min(f_o), max(f_o)+0.1]) + sp2_f.set_ylim([-120.0, 20.0]) + sp2_f.grid(True) + + sp2_f.set_title(("Channel %d" % i), weight="bold") + sp2_f.set_xlabel("Frequency (kHz)") + sp2_f.set_ylabel("Power (dBW)") + + + Ts = 1.0/fs_o + Tmax = len(d)*Ts + t_o = scipy.arange(0, Tmax, Ts) + + x_t = scipy.array(d) + sp2_t = fig3.add_subplot(Nrows, Ncols, 1+i) + p2_t = sp2_t.plot(t_o, x_t.real, "b") + p2_t = sp2_t.plot(t_o, x_t.imag, "r") + sp2_t.set_xlim([min(t_o), max(t_o)+1]) + sp2_t.set_ylim([-1, 1]) + + sp2_t.set_xlabel("Time (s)") + sp2_t.set_ylabel("Amplitude") + + + pylab.show() + + +if __name__ == "__main__": + main() diff --git a/gnuradio-core/src/examples/pfb/interpolate.py b/gnuradio-core/src/examples/pfb/interpolate.py new file mode 100755 index 000000000..370cf26a7 --- /dev/null +++ b/gnuradio-core/src/examples/pfb/interpolate.py @@ -0,0 +1,233 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys, time + +try: + import scipy + from scipy import fftpack +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab + from pylab import mlab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +class pfb_top_block(gr.top_block): + def __init__(self): + gr.top_block.__init__(self) + + self._N = 100000 # number of samples to use + self._fs = 2000 # initial sampling rate + self._interp = 5 # Interpolation rate for PFB interpolator + self._ainterp = 5.5 # Resampling rate for the PFB arbitrary resampler + + # Frequencies of the signals we construct + freq1 = 100 + freq2 = 200 + + # Create a set of taps for the PFB interpolator + # This is based on the post-interpolation sample rate + self._taps = gr.firdes.low_pass_2(self._interp, self._interp*self._fs, freq2+50, 50, + attenuation_dB=120, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Create a set of taps for the PFB arbitrary resampler + # The filter size is the number of filters in the filterbank; 32 will give very low side-lobes, + # and larger numbers will reduce these even farther + # The taps in this filter are based on a sampling rate of the filter size since it acts + # internally as an interpolator. + flt_size = 32 + self._taps2 = gr.firdes.low_pass_2(flt_size, flt_size*self._fs, freq2+50, 150, + attenuation_dB=120, window=gr.firdes.WIN_BLACKMAN_hARRIS) + + # Calculate the number of taps per channel for our own information + tpc = scipy.ceil(float(len(self._taps)) / float(self._interp)) + print "Number of taps: ", len(self._taps) + print "Number of filters: ", self._interp + print "Taps per channel: ", tpc + + # Create a couple of signals at different frequencies + self.signal1 = gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freq1, 0.5) + self.signal2 = gr.sig_source_c(self._fs, gr.GR_SIN_WAVE, freq2, 0.5) + self.signal = gr.add_cc() + + self.head = gr.head(gr.sizeof_gr_complex, self._N) + + # Construct the PFB interpolator filter + self.pfb = blks2.pfb_interpolator_ccf(self._interp, self._taps) + + # Construct the PFB arbitrary resampler filter + self.pfb_ar = blks2.pfb_arb_resampler_ccf(self._ainterp, self._taps2, flt_size) + self.snk_i = gr.vector_sink_c() + + #self.pfb_ar.pfb.print_taps() + #self.pfb.pfb.print_taps() + + # Connect the blocks + self.connect(self.signal1, self.head, (self.signal,0)) + self.connect(self.signal2, (self.signal,1)) + self.connect(self.signal, self.pfb) + self.connect(self.signal, self.pfb_ar) + self.connect(self.signal, self.snk_i) + + # Create the sink for the interpolated signals + self.snk1 = gr.vector_sink_c() + self.snk2 = gr.vector_sink_c() + self.connect(self.pfb, self.snk1) + self.connect(self.pfb_ar, self.snk2) + + +def main(): + tb = pfb_top_block() + + tstart = time.time() + tb.run() + tend = time.time() + print "Run time: %f" % (tend - tstart) + + + if 1: + fig1 = pylab.figure(1, figsize=(12,10), facecolor="w") + fig2 = pylab.figure(2, figsize=(12,10), facecolor="w") + fig3 = pylab.figure(3, figsize=(12,10), facecolor="w") + + Ns = 10000 + Ne = 10000 + + fftlen = 8192 + winfunc = scipy.blackman + + # Plot input signal + fs = tb._fs + + d = tb.snk_i.data()[Ns:Ns+Ne] + sp1_f = fig1.add_subplot(2, 1, 1) + + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_in = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_in = scipy.arange(-fs/2.0, fs/2.0, fs/float(X_in.size)) + p1_f = sp1_f.plot(f_in, X_in, "b") + sp1_f.set_xlim([min(f_in), max(f_in)+1]) + sp1_f.set_ylim([-200.0, 50.0]) + + + sp1_f.set_title("Input Signal", weight="bold") + sp1_f.set_xlabel("Frequency (Hz)") + sp1_f.set_ylabel("Power (dBW)") + + Ts = 1.0/fs + Tmax = len(d)*Ts + + t_in = scipy.arange(0, Tmax, Ts) + x_in = scipy.array(d) + sp1_t = fig1.add_subplot(2, 1, 2) + p1_t = sp1_t.plot(t_in, x_in.real, "b-o") + #p1_t = sp1_t.plot(t_in, x_in.imag, "r-o") + sp1_t.set_ylim([-2.5, 2.5]) + + sp1_t.set_title("Input Signal", weight="bold") + sp1_t.set_xlabel("Time (s)") + sp1_t.set_ylabel("Amplitude") + + + # Plot output of PFB interpolator + fs_int = tb._fs*tb._interp + + sp2_f = fig2.add_subplot(2, 1, 1) + d = tb.snk1.data()[Ns:Ns+(tb._interp*Ne)] + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_o = scipy.arange(-fs_int/2.0, fs_int/2.0, fs_int/float(X_o.size)) + p2_f = sp2_f.plot(f_o, X_o, "b") + sp2_f.set_xlim([min(f_o), max(f_o)+1]) + sp2_f.set_ylim([-200.0, 50.0]) + + sp2_f.set_title("Output Signal from PFB Interpolator", weight="bold") + sp2_f.set_xlabel("Frequency (Hz)") + sp2_f.set_ylabel("Power (dBW)") + + Ts_int = 1.0/fs_int + Tmax = len(d)*Ts_int + + t_o = scipy.arange(0, Tmax, Ts_int) + x_o1 = scipy.array(d) + sp2_t = fig2.add_subplot(2, 1, 2) + p2_t = sp2_t.plot(t_o, x_o1.real, "b-o") + #p2_t = sp2_t.plot(t_o, x_o.imag, "r-o") + sp2_t.set_ylim([-2.5, 2.5]) + + sp2_t.set_title("Output Signal from PFB Interpolator", weight="bold") + sp2_t.set_xlabel("Time (s)") + sp2_t.set_ylabel("Amplitude") + + + # Plot output of PFB arbitrary resampler + fs_aint = tb._fs * tb._ainterp + + sp3_f = fig3.add_subplot(2, 1, 1) + d = tb.snk2.data()[Ns:Ns+(tb._interp*Ne)] + X,freq = mlab.psd(d, NFFT=fftlen, noverlap=fftlen/4, Fs=fs, + window = lambda d: d*winfunc(fftlen), + scale_by_freq=True) + X_o = 10.0*scipy.log10(abs(fftpack.fftshift(X))) + f_o = scipy.arange(-fs_aint/2.0, fs_aint/2.0, fs_aint/float(X_o.size)) + p3_f = sp3_f.plot(f_o, X_o, "b") + sp3_f.set_xlim([min(f_o), max(f_o)+1]) + sp3_f.set_ylim([-200.0, 50.0]) + + sp3_f.set_title("Output Signal from PFB Arbitrary Resampler", weight="bold") + sp3_f.set_xlabel("Frequency (Hz)") + sp3_f.set_ylabel("Power (dBW)") + + Ts_aint = 1.0/fs_aint + Tmax = len(d)*Ts_aint + + t_o = scipy.arange(0, Tmax, Ts_aint) + x_o2 = scipy.array(d) + sp3_f = fig3.add_subplot(2, 1, 2) + p3_f = sp3_f.plot(t_o, x_o2.real, "b-o") + p3_f = sp3_f.plot(t_o, x_o1.real, "m-o") + #p3_f = sp3_f.plot(t_o, x_o2.imag, "r-o") + sp3_f.set_ylim([-2.5, 2.5]) + + sp3_f.set_title("Output Signal from PFB Arbitrary Resampler", weight="bold") + sp3_f.set_xlabel("Time (s)") + sp3_f.set_ylabel("Amplitude") + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-core/src/examples/pfb/reconstruction.py b/gnuradio-core/src/examples/pfb/reconstruction.py new file mode 100755 index 000000000..c7909f7a5 --- /dev/null +++ b/gnuradio-core/src/examples/pfb/reconstruction.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python + +import scipy, math, pylab +from scipy import fftpack +from gnuradio import gr, digital, blks2 + +fftlen = 8192 + +def main(): + N = 10000 + fs = 2000.0 + Ts = 1.0/fs + t = scipy.arange(0, N*Ts, Ts) + + # When playing with the number of channels, be careful about the filter + # specs and the channel map of the synthesizer set below. + nchans = 10 + + # Build the filter(s) + bw = 1000 + tb = 400 + proto_taps = gr.firdes.low_pass_2(1, nchans*fs, bw, tb, 80, + gr.firdes.WIN_BLACKMAN_hARRIS) + print "Filter length: ", len(proto_taps) + + + # Create a modulated signal + npwr = 0.01 + data = scipy.random.randint(0, 256, N) + rrc_taps = gr.firdes.root_raised_cosine(1, 2, 1, 0.35, 41) + + src = gr.vector_source_b(data.astype(scipy.uint8).tolist(), False) + mod = digital.bpsk_mod(samples_per_symbol=2) + chan = gr.channel_model(npwr) + rrc = gr.fft_filter_ccc(1, rrc_taps) + + # Split it up into pieces + channelizer = blks2.pfb_channelizer_ccf(nchans, proto_taps, 2) + + # Put the pieces back together again + syn_taps = [nchans*t for t in proto_taps] + synthesizer = gr.pfb_synthesizer_ccf(nchans, syn_taps, True) + src_snk = gr.vector_sink_c() + snk = gr.vector_sink_c() + + # Remap the location of the channels + # Can be done in synth or channelizer (watch out for rotattions in + # the channelizer) + synthesizer.set_channel_map([ 0, 1, 2, 3, 4, + 15, 16, 17, 18, 19]) + + tb = gr.top_block() + tb.connect(src, mod, chan, rrc, channelizer) + tb.connect(rrc, src_snk) + + vsnk = [] + for i in xrange(nchans): + tb.connect((channelizer,i), (synthesizer, i)) + + vsnk.append(gr.vector_sink_c()) + tb.connect((channelizer,i), vsnk[i]) + + tb.connect(synthesizer, snk) + tb.run() + + sin = scipy.array(src_snk.data()[1000:]) + sout = scipy.array(snk.data()[1000:]) + + + # Plot original signal + fs_in = nchans*fs + f1 = pylab.figure(1, figsize=(16,12), facecolor='w') + s11 = f1.add_subplot(2,2,1) + s11.psd(sin, NFFT=fftlen, Fs=fs_in) + s11.set_title("PSD of Original Signal") + s11.set_ylim([-200, -20]) + + s12 = f1.add_subplot(2,2,2) + s12.plot(sin.real[1000:1500], "o-b") + s12.plot(sin.imag[1000:1500], "o-r") + s12.set_title("Original Signal in Time") + + start = 1 + skip = 4 + s13 = f1.add_subplot(2,2,3) + s13.plot(sin.real[start::skip], sin.imag[start::skip], "o") + s13.set_title("Constellation") + s13.set_xlim([-2, 2]) + s13.set_ylim([-2, 2]) + + # Plot channels + nrows = int(scipy.sqrt(nchans)) + ncols = int(scipy.ceil(float(nchans)/float(nrows))) + + f2 = pylab.figure(2, figsize=(16,12), facecolor='w') + for n in xrange(nchans): + s = f2.add_subplot(nrows, ncols, n+1) + s.psd(vsnk[n].data(), NFFT=fftlen, Fs=fs_in) + s.set_title("Channel {0}".format(n)) + s.set_ylim([-200, -20]) + + # Plot reconstructed signal + fs_out = 2*nchans*fs + f3 = pylab.figure(3, figsize=(16,12), facecolor='w') + s31 = f3.add_subplot(2,2,1) + s31.psd(sout, NFFT=fftlen, Fs=fs_out) + s31.set_title("PSD of Reconstructed Signal") + s31.set_ylim([-200, -20]) + + s32 = f3.add_subplot(2,2,2) + s32.plot(sout.real[1000:1500], "o-b") + s32.plot(sout.imag[1000:1500], "o-r") + s32.set_title("Reconstructed Signal in Time") + + start = 2 + skip = 4 + s33 = f3.add_subplot(2,2,3) + s33.plot(sout.real[start::skip], sout.imag[start::skip], "o") + s33.set_title("Constellation") + s33.set_xlim([-2, 2]) + s33.set_ylim([-2, 2]) + + pylab.show() + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + pass + diff --git a/gnuradio-core/src/examples/pfb/resampler.py b/gnuradio-core/src/examples/pfb/resampler.py new file mode 100755 index 000000000..7b296ca71 --- /dev/null +++ b/gnuradio-core/src/examples/pfb/resampler.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys + +try: + import scipy +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +class mytb(gr.top_block): + def __init__(self, fs_in, fs_out, fc, N=10000): + gr.top_block.__init__(self) + + rerate = float(fs_out) / float(fs_in) + print "Resampling from %f to %f by %f " %(fs_in, fs_out, rerate) + + # Creating our own taps + taps = gr.firdes.low_pass_2(32, 32, 0.25, 0.1, 80) + + self.src = gr.sig_source_c(fs_in, gr.GR_SIN_WAVE, fc, 1) + #self.src = gr.noise_source_c(gr.GR_GAUSSIAN, 1) + self.head = gr.head(gr.sizeof_gr_complex, N) + + # A resampler with our taps + self.resamp_0 = blks2.pfb_arb_resampler_ccf(rerate, taps, + flt_size=32) + + # A resampler that just needs a resampling rate. + # Filter is created for us and designed to cover + # entire bandwidth of the input signal. + # An optional atten=XX rate can be used here to + # specify the out-of-band rejection (default=80). + self.resamp_1 = blks2.pfb_arb_resampler_ccf(rerate) + + self.snk_in = gr.vector_sink_c() + self.snk_0 = gr.vector_sink_c() + self.snk_1 = gr.vector_sink_c() + + self.connect(self.src, self.head, self.snk_in) + self.connect(self.head, self.resamp_0, self.snk_0) + self.connect(self.head, self.resamp_1, self.snk_1) + +def main(): + fs_in = 8000 + fs_out = 20000 + fc = 1000 + N = 10000 + + tb = mytb(fs_in, fs_out, fc, N) + tb.run() + + + # Plot PSD of signals + nfftsize = 2048 + fig1 = pylab.figure(1, figsize=(10,10), facecolor="w") + sp1 = fig1.add_subplot(2,1,1) + sp1.psd(tb.snk_in.data(), NFFT=nfftsize, + noverlap=nfftsize/4, Fs = fs_in) + sp1.set_title(("Input Signal at f_s=%.2f kHz" % (fs_in/1000.0))) + sp1.set_xlim([-fs_in/2, fs_in/2]) + + sp2 = fig1.add_subplot(2,1,2) + sp2.psd(tb.snk_0.data(), NFFT=nfftsize, + noverlap=nfftsize/4, Fs = fs_out, + label="With our filter") + sp2.psd(tb.snk_1.data(), NFFT=nfftsize, + noverlap=nfftsize/4, Fs = fs_out, + label="With auto-generated filter") + sp2.set_title(("Output Signals at f_s=%.2f kHz" % (fs_out/1000.0))) + sp2.set_xlim([-fs_out/2, fs_out/2]) + sp2.legend() + + # Plot signals in time + Ts_in = 1.0/fs_in + Ts_out = 1.0/fs_out + t_in = scipy.arange(0, len(tb.snk_in.data())*Ts_in, Ts_in) + t_out = scipy.arange(0, len(tb.snk_0.data())*Ts_out, Ts_out) + + fig2 = pylab.figure(2, figsize=(10,10), facecolor="w") + sp21 = fig2.add_subplot(2,1,1) + sp21.plot(t_in, tb.snk_in.data()) + sp21.set_title(("Input Signal at f_s=%.2f kHz" % (fs_in/1000.0))) + sp21.set_xlim([t_in[100], t_in[200]]) + + sp22 = fig2.add_subplot(2,1,2) + sp22.plot(t_out, tb.snk_0.data(), + label="With our filter") + sp22.plot(t_out, tb.snk_1.data(), + label="With auto-generated filter") + sp22.set_title(("Output Signals at f_s=%.2f kHz" % (fs_out/1000.0))) + r = float(fs_out)/float(fs_in) + sp22.set_xlim([t_out[r * 100], t_out[r * 200]]) + sp22.legend() + + pylab.show() + +if __name__ == "__main__": + main() + diff --git a/gnuradio-core/src/examples/pfb/resampler_demo.grc b/gnuradio-core/src/examples/pfb/resampler_demo.grc new file mode 100644 index 000000000..468636a5c --- /dev/null +++ b/gnuradio-core/src/examples/pfb/resampler_demo.grc @@ -0,0 +1,598 @@ + + + Sun Aug 23 11:39:47 2009 + + options + + id + resampler_demo + + + _enabled + True + + + title + + + + author + + + + description + + + + window_size + 1280, 1024 + + + generate_options + wx_gui + + + category + Custom + + + run + True + + + realtime_scheduling + + + + _coordinate + (10, 10) + + + _rotation + 0 + + + + import + + id + import_0 + + + _enabled + True + + + import + import math + + + _coordinate + (11, 59) + + + _rotation + 0 + + + + variable + + id + rs_taps + + + _enabled + True + + + value + firdes.low_pass(nphases, nphases, frac_bw, 0.5-frac_bw) + + + _coordinate + (273, 154) + + + _rotation + 0 + + + + gr_add_const_vxx + + id + adder + + + _enabled + True + + + type + float + + + const + -1.0 + + + vlen + 1 + + + _coordinate + (227, 303) + + + _rotation + 0 + + + + gr_throttle + + id + throttle + + + _enabled + True + + + type + float + + + samples_per_second + samp_rate + + + vlen + 1 + + + _coordinate + (227, 493) + + + _rotation + 0 + + + + wxgui_fftsink2 + + id + orig_fft + + + _enabled + True + + + type + complex + + + title + Original Spectrum + + + samp_rate + samp_rate + + + baseband_freq + 0 + + + y_per_div + 10 + + + y_divs + 10 + + + ref_level + 30 + + + fft_size + 1024 + + + fft_rate + 30 + + + peak_hold + False + + + average + False + + + avg_alpha + 0 + + + grid_pos + 1, 0, 1, 3 + + + notebook + + + + _coordinate + (409, 289) + + + _rotation + 180 + + + + wxgui_fftsink2 + + id + resamp_fft + + + _enabled + True + + + type + complex + + + title + Resampled Spectrum + + + samp_rate + new_rate + + + baseband_freq + 0 + + + y_per_div + 10 + + + y_divs + 10 + + + ref_level + 30 + + + fft_size + 1024 + + + fft_rate + 30 + + + peak_hold + True + + + average + False + + + avg_alpha + 0 + + + grid_pos + 2, 0, 1, 3 + + + notebook + + + + _coordinate + (640, 256) + + + _rotation + 180 + + + + gr_sig_source_x + + id + tri_source + + + _enabled + True + + + type + float + + + samp_rate + samp_rate + + + waveform + gr.GR_TRI_WAVE + + + freq + 0.05 + + + amp + 2.0 + + + offset + 0 + + + _coordinate + (21, 271) + + + _rotation + 0 + + + + gr_frequency_modulator_fc + + id + fm_mod + + + _enabled + True + + + sensitivity + math.pi + + + _coordinate + (411, 493) + + + _rotation + 0 + + + + blks2_pfb_arb_resampler_ccf + + id + resampler + + + _enabled + True + + + rate + float(new_rate)/samp_rate + + + taps + rs_taps + + + size + nphases + + + _coordinate + (641, 477) + + + _rotation + 0 + + + + variable + + id + nphases + + + _enabled + True + + + value + 32 + + + _coordinate + (185, 153) + + + _rotation + 0 + + + + variable_static_text + + id + samp_rate + + + _enabled + True + + + label + Sample Rate + + + value + 44100 + + + converver + float_converter + + + formatter + None + + + grid_pos + 0, 0, 1, 1 + + + notebook + + + + _coordinate + (179, 14) + + + _rotation + 0 + + + + variable_static_text + + id + new_rate + + + _enabled + True + + + label + Resampled Rate + + + value + 48000 + + + converver + float_converter + + + formatter + None + + + grid_pos + 0, 1, 1, 1 + + + notebook + + + + _coordinate + (328, 15) + + + _rotation + 0 + + + + variable_static_text + + id + frac_bw + + + _enabled + True + + + label + Fractional Bandwidth + + + value + 0.45 + + + converver + float_converter + + + formatter + lambda x: "%0.2f"%x + + + grid_pos + 0,2,1,1 + + + notebook + + + + _coordinate + (473, 14) + + + _rotation + 0 + + + + tri_source + adder + 0 + 0 + + + adder + throttle + 0 + 0 + + + resampler + resamp_fft + 0 + 0 + + + fm_mod + resampler + 0 + 0 + + + fm_mod + orig_fft + 0 + 0 + + + throttle + fm_mod + 0 + 0 + + diff --git a/gnuradio-core/src/examples/pfb/synth_filter.py b/gnuradio-core/src/examples/pfb/synth_filter.py new file mode 100755 index 000000000..a91edfebf --- /dev/null +++ b/gnuradio-core/src/examples/pfb/synth_filter.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python +# +# Copyright 2010 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys + +try: + import scipy +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +def main(): + N = 1000000 + fs = 8000 + + freqs = [100, 200, 300, 400, 500] + nchans = 7 + + sigs = list() + for fi in freqs: + s = gr.sig_source_c(fs, gr.GR_SIN_WAVE, fi, 1) + sigs.append(s) + + taps = gr.firdes.low_pass_2(len(freqs), fs, fs/float(nchans)/2, 100, 100) + print "Num. Taps = %d (taps per filter = %d)" % (len(taps), + len(taps)/nchans) + filtbank = gr.pfb_synthesizer_ccf(nchans, taps) + + head = gr.head(gr.sizeof_gr_complex, N) + snk = gr.vector_sink_c() + + tb = gr.top_block() + tb.connect(filtbank, head, snk) + + for i,si in enumerate(sigs): + tb.connect(si, (filtbank, i)) + + tb.run() + + if 1: + f1 = pylab.figure(1) + s1 = f1.add_subplot(1,1,1) + s1.plot(snk.data()[1000:]) + + fftlen = 2048 + f2 = pylab.figure(2) + s2 = f2.add_subplot(1,1,1) + winfunc = scipy.blackman + s2.psd(snk.data()[10000:], NFFT=fftlen, + Fs = nchans*fs, + noverlap=fftlen/4, + window = lambda d: d*winfunc(fftlen)) + + pylab.show() + +if __name__ == "__main__": + main() diff --git a/gnuradio-core/src/examples/pfb/synth_to_chan.py b/gnuradio-core/src/examples/pfb/synth_to_chan.py new file mode 100755 index 000000000..c6c80b2f8 --- /dev/null +++ b/gnuradio-core/src/examples/pfb/synth_to_chan.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python +# +# Copyright 2010 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +from gnuradio import gr, blks2 +import sys + +try: + import scipy +except ImportError: + print "Error: Program requires scipy (see: www.scipy.org)." + sys.exit(1) + +try: + import pylab +except ImportError: + print "Error: Program requires matplotlib (see: matplotlib.sourceforge.net)." + sys.exit(1) + +def main(): + N = 1000000 + fs = 8000 + + freqs = [100, 200, 300, 400, 500] + nchans = 7 + + sigs = list() + fmtx = list() + for fi in freqs: + s = gr.sig_source_f(fs, gr.GR_SIN_WAVE, fi, 1) + fm = blks2.nbfm_tx (fs, 4*fs, max_dev=10000, tau=75e-6) + sigs.append(s) + fmtx.append(fm) + + syntaps = gr.firdes.low_pass_2(len(freqs), fs, fs/float(nchans)/2, 100, 100) + print "Synthesis Num. Taps = %d (taps per filter = %d)" % (len(syntaps), + len(syntaps)/nchans) + chtaps = gr.firdes.low_pass_2(len(freqs), fs, fs/float(nchans)/2, 100, 100) + print "Channelizer Num. Taps = %d (taps per filter = %d)" % (len(chtaps), + len(chtaps)/nchans) + filtbank = gr.pfb_synthesizer_ccf(nchans, syntaps) + channelizer = blks2.pfb_channelizer_ccf(nchans, chtaps) + + noise_level = 0.01 + head = gr.head(gr.sizeof_gr_complex, N) + noise = gr.noise_source_c(gr.GR_GAUSSIAN, noise_level) + addnoise = gr.add_cc() + snk_synth = gr.vector_sink_c() + + tb = gr.top_block() + + tb.connect(noise, (addnoise,0)) + tb.connect(filtbank, head, (addnoise, 1)) + tb.connect(addnoise, channelizer) + tb.connect(addnoise, snk_synth) + + snk = list() + for i,si in enumerate(sigs): + tb.connect(si, fmtx[i], (filtbank, i)) + + for i in xrange(nchans): + snk.append(gr.vector_sink_c()) + tb.connect((channelizer, i), snk[i]) + + tb.run() + + if 1: + channel = 1 + data = snk[channel].data()[1000:] + + f1 = pylab.figure(1) + s1 = f1.add_subplot(1,1,1) + s1.plot(data[10000:10200] ) + s1.set_title(("Output Signal from Channel %d" % channel)) + + fftlen = 2048 + winfunc = scipy.blackman + #winfunc = scipy.hamming + + f2 = pylab.figure(2) + s2 = f2.add_subplot(1,1,1) + s2.psd(data, NFFT=fftlen, + Fs = nchans*fs, + noverlap=fftlen/4, + window = lambda d: d*winfunc(fftlen)) + s2.set_title(("Output PSD from Channel %d" % channel)) + + f3 = pylab.figure(3) + s3 = f3.add_subplot(1,1,1) + s3.psd(snk_synth.data()[1000:], NFFT=fftlen, + Fs = nchans*fs, + noverlap=fftlen/4, + window = lambda d: d*winfunc(fftlen)) + s3.set_title("Output of Synthesis Filter") + + pylab.show() + +if __name__ == "__main__": + main() -- cgit