Merged r9564:9568 from jcorgan/bert into trunk. Adds BPSK bit error testing examples and support blocks.

git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@9569 221aa14e-8319-0410-a670-987f0aec2ac5

6 files changed, 548 insertions, 0 deletions

diff --git a/gnuradio-examples/python/digital-bert/Makefile.am b/gnuradio-examples/python/digital-bert/Makefile.am
new file mode 100644
index 000000000..886033792
--- /dev/null
+++ b/gnuradio-examples/python/digital-bert/Makefile.am
@@ -0,0 +1,33 @@
+#
+# Copyright 2008 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 $(top_srcdir)/Makefile.common
+
+EXTRA_DIST =			\
+	benchmark_rx.py		\
+	benchmark_tx.py		\
+	receive_path.py		\
+	transmit_path.py
+
+ourdatadir = $(exampledir)/digital-bert
+ourdata_DATA = $(EXTRA_DIST)
+
+MOSTLYCLEANFILES = *.pyc *.pyo *~
diff --git a/gnuradio-examples/python/digital-bert/README b/gnuradio-examples/python/digital-bert/README
new file mode 100644
index 000000000..74d85b703
--- /dev/null
+++ b/gnuradio-examples/python/digital-bert/README
@@ -0,0 +1,62 @@
+BERT testing example scripts
+
+
+benchmark_tx.py
+
+This sets up a BPSK transmitter that is modulated with a pseudorandom
+sequence of bits.  The PN code is generated by sending an all 1s
+sequence through a 7-bit scrambler. The transmitter performs the BPSK
+modulation, then passes the complex baseband waveform through a
+root-raised-cosine filter and onto the USRP.
+
+The --sps parameter controls how many baseband samples per symbol
+are created and passed through the RRC filter, prior to going to the
+USRP over the USB for interpolation to the final DAC rate.
+
+The baseband bit rate is controlled by -r or --rate.  This value, when
+multiplied by the --sps parameter, must result in valid interpolation
+rate for the USRP.  For example, if the baseband rate is 250k bits/sec,
+and the samples per symbol is 4, then the final rate is 1M samples/sec,
+which results in an interpolation rate of 128.  The valid interpolation
+rates for the USRP are multiples of 4 between 16 and 512.
+
+Finally, the RRC excess bandwidth may be specified by --excess-bw.
+(See ./benchmark_tx.py -h for additional parameters.)
+
+
+benchmark_rx.py
+
+This sets up a BPSK receiver to demodulate the received waveform.  It
+accepts a similar set of parameters as the transmitter, except that one
+specifies the USRP decimation rate desired.  The resulting sample stream
+rate must be an integral number of baseband symbols.  For example, the
+parameters corresponding to the above transmitter would be to use a
+decimation rate of 8 (32 sps), 16 (16 sps), 32 (8 sps), 64, (4 sps), or
+128 (2 sps).  The lower the USRP decimation, the more CPU is required to
+demodulate the signal, so not all valid decimation rates will work.
+
+The baseband signal from the USRP is first passed through an AGC to
+establish an average power of 1.0.  It is then passed through a matched
+filter (another RRC), a Costas phase-locked loop, and an M&M bit timing
+recovery loop.  The resulting constellation has an SNR estimation probe
+attached, and is then sliced into a bit stream.
+
+The recovered bits are then passed through a 7-bit descrambler.  If
+there are no channel errors, the all 1s sequence is recovered. In the 
+event of a channel error, there will be a 0 in the bit stream for each
+feedback tap in the descrambler.  In this case, the CCSDS descrambler is
+using 3 feedback taps.
+
+Finally, the signal is passed into a bit density measurement probe. The
+channel BER is measured by dividing the 0s density by three.  This
+measurement is inaccurate at high BER rates (>10%) as the error 0s
+begin to overlap.
+
+The benchmark script will, once per second, output the mean symbol power,
+the average noise variance, the estimated SNR, and the average BER.
+
+NOTE: The particular SNR estimator used is inaccurate below about 7dB,
+and will report erroneously high values even for random noise.
+
+There are a variety of Costas and M&M loop parameters one can adjust.
+See ./benchmark_rx.py -h for the full set.
diff --git a/gnuradio-examples/python/digital-bert/benchmark_rx.py b/gnuradio-examples/python/digital-bert/benchmark_rx.py
new file mode 100755
index 000000000..96bf7fa12
--- /dev/null
+++ b/gnuradio-examples/python/digital-bert/benchmark_rx.py
@@ -0,0 +1,164 @@
+#!/usr/bin/env python
+#
+# Copyright 2008 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, eng_notation, usrp
+from optparse import OptionParser
+from gnuradio.eng_option import eng_option
+from receive_path import receive_path
+import gnuradio.gr.gr_threading as _threading
+import sys, time
+
+n2s = eng_notation.num_to_str
+
+class status_thread(_threading.Thread):
+    def __init__(self, tb):
+        _threading.Thread.__init__(self)
+        self.setDaemon(1)
+        self.tb = tb
+        self.done = False
+        self.start()
+
+    def run(self):
+        while not self.done:
+            print "Magnitude: %5.3f  Variance: %5.3f  Estimated SNR: %4.1f dB  BER: %g" % (
+                tb.mag(), tb.var(), tb.snr(), tb.ber())
+            try:
+                time.sleep(1.0)
+            except KeyboardInterrupt:
+                self.done = True
+
+class rx_bpsk_block(gr.top_block):
+    def __init__(self, options):    
+
+	gr.top_block.__init__(self, "rx_mpsk")
+
+        print "USRP decimation rate", options.decim_rate
+
+        # Create a USRP source at desired board, sample rate, frequency, and gain
+        self._setup_usrp(options.which,
+                         options.decim_rate,
+                         options.rx_subdev_spec,
+                         options.freq,
+                         options.gain)
+
+        # Create the BERT receiver
+        if_rate = self._usrp.adc_rate()/options.decim_rate
+        self._receiver = receive_path(if_rate,
+                                      options.rate,
+                                      options.excess_bw,
+                                      options.costas_alpha,
+                                      options.costas_beta,
+                                      options.costas_max,
+                                      options.mm_gain_mu,
+                                      options.mm_gain_omega,
+                                      options.mm_omega_limit)
+        
+        self.connect(self._usrp, self._receiver)
+
+
+    def _setup_usrp(self, which, decim, subdev_spec, freq, gain):
+        self._usrp = usrp.source_c(which=which, decim_rate=decim)
+	if subdev_spec is None:
+	    subdev_spec = usrp.pick_rx_subdevice(self._usrp)
+	self._subdev = usrp.selected_subdev(self._usrp, subdev_spec)
+        mux = usrp.determine_rx_mux_value(self._usrp, subdev_spec)
+	self._usrp.set_mux(mux)
+	tr = self._usrp.tune(0, self._subdev, freq)
+        if not (tr):
+            print "Failed to tune to center frequency!"
+        else:
+            print "Center frequency:", n2s(freq)
+        if gain is None:
+            g = self._subdev.gain_range();
+            gain = float(g[0]+g[1])/2.0
+	self._subdev.set_gain(gain)
+        print "RX d'board:", self._subdev.side_and_name()
+        
+    def snr(self):
+        return self._receiver.snr()
+
+    def mag(self):
+        return self._receiver.signal_mean()
+        
+    def var(self):
+        return self._receiver.noise_variance()
+
+    def ber(self):
+        return self._receiver.ber()
+            
+def get_options():
+    parser = OptionParser(option_class=eng_option)
+    parser.add_option("-w", "--which", type="int", default=0,
+                      help="select which USRP (0, 1, ...) (default is %default)",
+                      metavar="NUM")
+    parser.add_option("-R", "--rx-subdev-spec", type="subdev", default=None,
+                      help="select USRP Rx side A or B (default=first one with a daughterboard)")
+    parser.add_option("-f", "--freq", type="eng_float", default=None,
+                      help="set frequency to FREQ", metavar="FREQ")
+    parser.add_option("-g", "--gain", type="eng_float", default=None,
+                      help="set Rx gain (default is mid-point)")
+    parser.add_option("-r", "--rate", type="eng_float", default=250e3,
+                      help="Select modulation symbol rate (default=%default)")
+    parser.add_option("-d", "--decim-rate", type="int", default=8,
+                      help="Select USRP decimation rate (default=%default)")
+    parser.add_option("", "--excess-bw", type="eng_float", default=0.35,
+                      help="Select RRC excess bandwidth (default=%default)")
+    parser.add_option("", "--costas-alpha", type="eng_float", default=0.2,
+                      help="set Costas loop 1st order gain, (default=%default)")
+    parser.add_option("", "--costas-beta", type="eng_float", default=0.01,
+                      help="set Costas loop 2nd order gain, (default=%default)")
+    parser.add_option("", "--costas-max", type="eng_float", default=0.005,
+                      help="set Costas loop max freq (rad/sample) (default=%default)")
+    parser.add_option("", "--mm-gain-mu", type="eng_float", default=0.2,
+                      help="set M&M loop 1st order gain, (default=%default)")
+    parser.add_option("", "--mm-gain-omega", type="eng_float", default=0.01,
+                      help="set M&M loop 2nd order gain, (default=%default)")
+    parser.add_option("", "--mm-omega-limit", type="eng_float", default=0.005,
+                      help="set M&M max timing error, (default=%default)")
+
+		      
+    (options, args) = parser.parse_args()
+    if len(args) != 0:
+        parser.print_help()
+        sys.exit(1)
+	
+    if options.freq == None:
+        print "You must supply a frequency with -f or --freq"
+        sys.exit(1)
+
+    return (options, args)
+
+
+if __name__ == "__main__":
+    (options, args) = get_options()
+
+    tb = rx_bpsk_block(options)
+
+    print "\n*** SNR estimator is inaccurate below about 7dB"
+    print "*** BER estimator is inaccurate above about 10%\n"
+    updater = status_thread(tb)
+
+    try:
+        tb.run()
+    except KeyboardInterrupt:
+        updater.done = True
+        updater = None
diff --git a/gnuradio-examples/python/digital-bert/benchmark_tx.py b/gnuradio-examples/python/digital-bert/benchmark_tx.py
new file mode 100755
index 000000000..1778a74b5
--- /dev/null
+++ b/gnuradio-examples/python/digital-bert/benchmark_tx.py
@@ -0,0 +1,111 @@
+#!/usr/bin/env python
+#
+# Copyright 2008 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, eng_notation, usrp
+from gnuradio.eng_option import eng_option
+from optparse import OptionParser
+from transmit_path import transmit_path
+import sys
+
+_dac_rate = 128e6
+
+n2s = eng_notation.num_to_str
+
+class tx_bpsk_block(gr.top_block):
+    def __init__(self, options):    
+	gr.top_block.__init__(self, "tx_mpsk")
+
+        self._transmitter = transmit_path(options.sps,
+                                          options.excess_bw,
+                                          options.amplitude)
+
+        if_rate = options.rate*options.sps
+        interp = int(_dac_rate/if_rate)
+
+        print "Modulation:", n2s(options.rate), "bits/sec"
+        print "TX IF rate:", n2s(if_rate), "samples/sec"
+        print "USRP interpolation:", interp
+        print "DAC amplitude:", options.amplitude
+
+        self._setup_usrp(options.which,
+                         interp,
+                         options.tx_subdev_spec,
+                         options.freq)
+	
+	self.connect(self._transmitter, self._usrp)
+
+        
+    def _setup_usrp(self, which, interp, subdev_spec, freq):
+        self._usrp = usrp.sink_c(which=which, interp_rate=interp)
+	if subdev_spec is None:
+	    subdev_spec = usrp.pick_tx_subdevice(self._usrp)
+	self._usrp.set_mux(usrp.determine_tx_mux_value(self._usrp, subdev_spec))
+	self._subdev = usrp.selected_subdev(self._usrp, subdev_spec)
+	tr = usrp.tune(self._usrp, self._subdev._which, self._subdev, freq)
+        if not (tr):
+            print "Failed to tune to center frequency!"
+        else:
+            print "Center frequency:", n2s(freq)
+        gain = float(self._subdev.gain_range()[1]) # Max TX gain
+	self._subdev.set_gain(gain)
+	self._subdev.set_enable(True)
+        print "TX d'board:", self._subdev.side_and_name()
+        
+
+def get_options():
+    parser = OptionParser(option_class=eng_option)
+    parser.add_option("-w", "--which", type="int", default=0,
+                      help="select which USRP (0, 1, ...) default is %default",
+                      metavar="NUM")
+    parser.add_option("-T", "--tx-subdev-spec", type="subdev", default=None,
+                      help="select USRP Tx side A or B (default=first one with a daughterboard)")
+    parser.add_option("-f", "--freq", type="eng_float", default=None,
+                      help="set frequency to FREQ", metavar="FREQ")
+    parser.add_option("-a", "--amplitude", type="eng_float", default=2000,
+                      help="set Tx amplitude (0-32767) (default=%default)")
+    parser.add_option("-r", "--rate", type="eng_float", default=250e3,
+                      help="Select modulation symbol rate (default=%default)")
+    parser.add_option("", "--sps", type="int", default=2,
+                      help="Select samples per symbol (default=%default)")
+    parser.add_option("", "--excess-bw", type="eng_float", default=0.35,
+                      help="Select RRC excess bandwidth (default=%default)")
+		      
+    (options, args) = parser.parse_args()
+    if len(args) != 0:
+        parser.print_help()
+        sys.exit(1)
+	
+    if options.freq == None:
+        print "Must supply frequency as -f or --freq"
+        sys.exit(1)
+
+    return (options, args)
+
+if __name__ == "__main__":
+    (options, args) = get_options()
+    
+    tb = tx_bpsk_block(options)
+
+    try:
+        tb.run()
+    except KeyboardInterrupt:
+        pass
diff --git a/gnuradio-examples/python/digital-bert/receive_path.py b/gnuradio-examples/python/digital-bert/receive_path.py
new file mode 100644
index 000000000..84ff32916
--- /dev/null
+++ b/gnuradio-examples/python/digital-bert/receive_path.py
@@ -0,0 +1,116 @@
+#
+# Copyright 2008 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, eng_notation
+import math
+
+n2s = eng_notation.num_to_str
+
+class receive_path(gr.hier_block2):
+    def __init__(self,
+                 if_rate,        # Incoming sample rate
+                 symbol_rate,    # Original symbol rate
+                 excess_bw,      # RRC excess bandwidth, typically 0.35-0.5
+                 costas_alpha,   # Costas loop 1st order gain, typically 0.01-0.2
+                 costas_beta,    # Costas loop 2nd order gain, typically alpha^2/4.0
+                 costas_max,     # Costas loop max frequency offset in radians/sample
+                 mm_gain_mu,     # M&M loop 1st order gain, typically 0.001-0.2
+                 mm_gain_omega,  # M&M loop 2nd order gain, typically alpha^2/4.0
+                 mm_omega_limit, # M&M loop max timing error
+                 ):
+        
+        gr.hier_block2.__init__(self, "receive_path",
+                                gr.io_signature(1, 1, gr.sizeof_gr_complex), # Input signature
+                                gr.io_signature(0, 0, 0))                    # Output signature
+
+        sps = int(if_rate/symbol_rate)
+        print "IF sample rate:", n2s(if_rate)
+        print "Symbol rate:", n2s(symbol_rate)
+        print "Samples/symbol:", sps
+        print "RRC bandwidth:", excess_bw
+        
+        # Create AGC to scale input to unity
+        self._agc = gr.agc_cc(1e-5, 1.0, 1.0, 1.0)
+
+	# Create RRC with specified excess bandwidth
+	taps = gr.firdes.root_raised_cosine(1.0,       # Gain
+					    sps,       # Sampling rate
+					    1.0,       # Symbol rate
+					    excess_bw, # Roll-off factor
+					    11*sps)    # Number of taps
+
+	self._rrc = gr.fir_filter_ccf(1, taps)
+        
+        # Create a Costas loop frequency/phase recovery block
+
+        print "Costas alpha:", costas_alpha
+        print "Costas beta:", costas_beta
+        print "Costas max:", costas_max
+        
+        self._costas = gr.costas_loop_cc(costas_alpha,  # PLL first order gain
+                                         costas_beta,   # PLL second order gain
+                                         costas_max,    # Max frequency offset rad/sample
+                                         -costas_max,   # Min frequency offset rad/sample
+                                         2)             # BPSK
+
+        # Create a M&M bit synchronization retiming block
+        mm_mu = 0.5
+        mm_omega = sps
+
+        print "MM gain mu:", mm_gain_mu
+        print "MM gain omega:", mm_gain_omega
+        print "MM omega limit:", mm_omega_limit
+        
+        self._mm = gr.clock_recovery_mm_cc(mm_omega,       # Initial samples/symbol
+                                           mm_gain_omega,  # Second order gain
+                                           mm_mu,          # Initial symbol phase
+                                           mm_gain_mu,     # First order gain
+                                           mm_omega_limit) # Maximum timing offset
+
+        # Add an SNR probe on the demodulated constellation
+        self._snr_probe = gr.probe_mpsk_snr_c(10.0/symbol_rate)
+        self.connect(self._mm, self._snr_probe)
+        
+        # Slice the resulting constellation into bits.
+        # Get inphase channel and make decision about 0
+        self._c2r = gr.complex_to_real()
+        self._slicer = gr.binary_slicer_fb() 
+        
+        # Descramble BERT sequence.  A channel error will create 3 incorrect bits
+        self._descrambler = gr.descrambler_bb(0x8A, 0x7F, 7) # CCSDS 7-bit descrambler
+
+        # Measure BER by the density of 0s in the stream
+        self._ber = gr.probe_density_b(1.0/symbol_rate)
+
+        self.connect(self, self._agc, self._rrc, self._costas, self._mm, 
+                     self._c2r, self._slicer, self._descrambler, self._ber)
+
+    def snr(self):
+        return self._snr_probe.snr()
+
+    def signal_mean(self):
+        return self._snr_probe.signal_mean()
+        
+    def noise_variance(self):
+        return self._snr_probe.noise_variance()
+
+    def ber(self):
+        return (1.0-self._ber.density())/3.0
diff --git a/gnuradio-examples/python/digital-bert/transmit_path.py b/gnuradio-examples/python/digital-bert/transmit_path.py
new file mode 100644
index 000000000..47b6ab229
--- /dev/null
+++ b/gnuradio-examples/python/digital-bert/transmit_path.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+#
+# Copyright 2008 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
+from math import pi, log10
+import cmath
+
+class transmit_path(gr.hier_block2):
+    """
+    This transmits a BERT sequence of bits using filtered BPSK and
+    outputs the complex baseband waveform.
+    """
+    def __init__(self,
+                 sps,          # Samples per symbol
+                 excess_bw,    # RRC filter excess bandwidth (typically 0.35-0.5)
+                 amplitude     # DAC output level, 0-32767, typically 2000-8000
+                 ):
+
+        gr.hier_block2.__init__(self, "transmit_path",
+                                gr.io_signature(0, 0, 0),                    # Input signature
+                                gr.io_signature(1, 1, gr.sizeof_gr_complex)) # Output signature
+
+        # Create BERT data bit stream	
+	self._bits = gr.vector_source_b([1,], True)      # Infinite stream of ones
+        self._scrambler = gr.scrambler_bb(0x8A, 0x7F, 7) # CCSDS 7-bit scrambler
+
+	# Map to constellation
+	self._constellation = [-1+0j, 1+0j]
+	self._mapper = gr.chunks_to_symbols_bc(self._constellation)	
+
+	# Create RRC with specified excess bandwidth
+	taps = gr.firdes.root_raised_cosine(sps*amplitude,  # Gain
+					    sps,	    # Sampling rate
+					    1.0,	    # Symbol rate
+					    excess_bw,	    # Roll-off factor
+					    11*sps)         # Number of taps
+
+	self._rrc = gr.interp_fir_filter_ccf(sps,	    # Interpolation rate
+					     taps)	    # FIR taps
+
+        # Wire block inputs and outputs
+        self.connect(self._bits, self._scrambler, self._mapper, self._rrc, self)
+