Updated documentation in gr-trellis/doc

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

3 files changed, 436 insertions, 8 deletions

diff --git a/gr-trellis/doc/gr-trellis.xml b/gr-trellis/doc/gr-trellis.xml
index e33a94a2d..ed53715a8 100644
--- a/gr-trellis/doc/gr-trellis.xml
+++ b/gr-trellis/doc/gr-trellis.xml
@@ -2,6 +2,7 @@
 <!DOCTYPE article PUBLIC "-//OASIS//DTD DocBook XML V4.2//EN"
           "docbookx.dtd" [
   <!ENTITY test_tcm_listing SYSTEM "test_tcm.py.xml">
+  <!ENTITY test_viterbi_equalization1_listing SYSTEM "test_viterbi_equalization1.py.xml">
 ]>
 
 <article>
@@ -18,6 +19,7 @@
      </affiliation>
   </author>
 
+<!--
 <revhistory>
   <revision>
   <revnumber>v0.0</revnumber>
@@ -26,8 +28,8 @@
     First cut.
   </revremark>
   </revision>
-
 </revhistory>
+-->
 
 <abstract><para>This document provides a description of the 
 Finite State Machine (FSM) implementation and the related 
@@ -43,8 +45,6 @@ for GNU Radio.
 <!--=====================================================-->
 <sect1 id="intro"><title>Introduction</title>
 
-<para>....</para>
-
 <para>
 The basic goal of the implementation is to have a generic way of 
 describing an FSM that is decoupled from whether it describes a 
@@ -342,14 +342,14 @@ In this section we give a brief description of the basic blocks implemented that
 <!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -->
 <sect2 id="encoder"><title>Trellis Encoder</title>
 <para>
-The trellis.encoder_XX(FSM, ST) block instantiates an FSM encoder corresponding to the fsm FSM and having initial state ST. The input and output is a sequence of bytes, shorts or integers. 
+The <methodname>trellis.encoder_XX(FSM, ST)</methodname> block instantiates an FSM encoder corresponding to the fsm FSM and having initial state ST. The input and output is a sequence of bytes, shorts or integers. 
 </para>
 </sect2>
 
 <!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -->
 <sect2 id="decoder"><title>Viterbi Decoder</title>
 <para>
-The trellis.viterbi_X(FSM, K, S0, SK) block instantiates a Viterbi decoder 
+The <methodname>trellis.viterbi_X(FSM, K, S0, SK)</methodname> block instantiates a Viterbi decoder 
 for a sequence of K trellis steps generated by the given FSM and with initial and final states set to S0 and SK, respectively (S0 and/or SK are set to -1
 if the corresponding states are not fixed/known at the receiver side).
 The output of this block is a sequence of K bytes, shorts or integers representing the estimated input (i.e., uncoded) sequence.
@@ -363,7 +363,7 @@ Observe that these inputs are generated externally and thus the Viterbi block is
 <!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -->
 <sect2 id="metrics"><title>Metrics Calculator</title>
 <para>
-The trellis.metrics_X(O,D,TABLE,TYPE) block is responsible for 
+The <methodname>trellis.metrics_X(O,D,TABLE,TYPE)</methodname> block is responsible for 
 transforming the channel output to the stream of metrics appropriate as 
 inputs to the Viterbi block described above. For each D input bytes/shorts/integers/floats/complexes it produces O output floats 
 
@@ -442,7 +442,7 @@ output of the metric block/input of the Viterbi block is FSM.O( ) floats for
 each trellis step. Sometimes this results in buffer overflow even for
 moderate sequence lengths.
 To overcome this problem we provide a block that incorporates the metric calculation and Viterbi algorithm into a single GNU Radio block, namely
-trellis.viterbi_combined_X( FSM, K, S0, SK, D, TABLE, TYPE) where the arguments are exactly those used in the aforementioned two blocks.
+<methodname>trellis.viterbi_combined_X( FSM, K, S0, SK, D, TABLE, TYPE)</methodname> where the arguments are exactly those used in the aforementioned two blocks.
 </para>
 </sect2>
 
@@ -454,7 +454,7 @@ trellis.viterbi_combined_X( FSM, K, S0, SK, D, TABLE, TYPE) where the arguments
 
 
 <!--=====================================================-->
-<sect1 id="tcm"><title>TCM: A Complete Example</title>
+<sect1 id="tcm"><title>A Complete Example: Trellis Coded Modulation (TCM)</title>
 
 <para>
 We now discuss through a concrete example how
@@ -652,11 +652,232 @@ The function returns the number of shorts and the number of shorts in error. Obs
  48      return (ntotal,ntotal-nright)
 </programlisting>
 
+</sect1>
+
+
+
+
+
+
+
+
+
+
+
+
+
+<!--=====================================================-->
+<sect1 id="isi"><title>Another Complete Example: Viterbi Equalization</title>
+
+<para>
+We now discuss through another concrete example how
+the above FSM model can be used in GNU Radio.
+
+The communication system that we want to simulate
+consists of a source generating the
+input information in packets, an ISI channel with
+additive white Gaussian noise (AWGN), and
+the VA performing MLSD.
+The program source is as follows.
+</para>
+
+&test_viterbi_equalization1_listing;
+
+<para>
+The program is called by
+<literallayout>
+./test_viterbi_equalization1.py Es/No_db repetitions
+</literallayout>
+where
+"Es/No_db" is the SNR in dB, and "repetitions" 
+are the number of packets to be transmitted and received in order to
+collect sufficient number of errors for an accurate estimate of the
+error rate.
+</para>
+
+
+<para>
+Each packet has size Kb bits.
+The modulation is chosen to be 4-PAM in this example and the channel is chosen 
+to be one of the test channels  defined in fsm_utils.py
+</para>
+<programlisting>
+ 71      Kb=2048  # packet size in bits
+ 72      modulation = fsm_utils.pam4 # see fsm_utlis.py for available predefined modulations
+ 73      channel = fsm_utils.c_channel # see fsm_utlis.py for available predefined test channels
+</programlisting>
+
+<para>
+The FSM is instantiated in
+</para>
+<programlisting>
+ 74      f=trellis.fsm(len(modulation[1]),len(channel)) # generate the FSM automatically
+</programlisting>
+<para>
+and generated automatically given the channel length and the modulation size.
+Since in this example the channel has length 5 and the modulation is 4-ary, the corresponding FSM has 4<superscript>5-1</superscript>=256 states and
+4<superscript>5</superscript>=1024 outputs (see the documentation on FSM for more explanation).
+</para>
+
+<para>
+Assuming that the FSM input has cardinality I, each input symbol consists 
+of bitspersymbol=log<subscript>2</subscript>( I ) bits, and thus correspond to K = Kb/bitspersymbol symbols. 
+</para>
+<programlisting>
+ 75      bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
+ 76      K=Kb/bitspersymbol # packet size in trellis steps
+</programlisting>
+
+
+
+<para>
+The overall system with input the 4-ary input symbols 
+x<subscript>k</subscript>, modulated to the
+4-PAM symbols s<subscript>k</subscript> and passed through the ISI channel to produce the
+noise-free observations 
+z<subscript>k</subscript> = 
+sum<subscript>j=0</subscript><superscript>L-1</superscript> c<subscript>j</subscript> s<subscript>k-j</subscript> (where L is the channel length)
+can be modeled as a FSM followed by a memoryless modulation. 
+In particular, the FSM input is the sequence 
+x<subscript>k</subscript>
+and its output is the "combined" symbol 
+y<subscript>k</subscript>=
+(x<subscript>k</subscript>,x<subscript>k-1</subscript>,...,x<subscript>k-L+1</subscript>) (actually its decimal representation). 
+The memoryless modulator maps every "combined" symbol
+y<subscript>k</subscript> to 
+z<subscript>k</subscript> = 
+sum<subscript>j=0</subscript><superscript>L-1</superscript> c<subscript>j</subscript> s<subscript>k-j</subscript> 
+Clearly this modulation is memoryless since 
+each z<subscript>k</subscript> depends only on y<subscript>k</subscript>; the memory inherent in the ISI is hidden in the FSM structure.
+This memoryless modulator is automatically generated by a helper function in 
+fsm_utils.py given the channel and modulation as in line 78, and has the 
+familiar format tot_channel=(dimensionality,tot_constellation) as described in the TCM example.
+This is exactly what the metrics block (or the viterbi_combined block) require in order to evaluate the VA metrics from the noisy observations.
+</para>
+<programlisting>
+ 78      tot_channel = fsm_utils.make_isi_lookup(modulation,channel,True) # generate the lookup table (normalize energy to 1)
+ 79      dimensionality = tot_channel[0]
+ 80      tot_constellation = tot_channel[1]
+ 81      N0=pow(10.0,-esn0_db/10.0); # noise variance
+ 82      if len(tot_constellation)/dimensionality != f.O():
+ 83          sys.stderr.write (&apos;Incompatible FSM output cardinality and lookup table size.\n&apos;)
+ 84          sys.exit (1)
+</programlisting>
+
+
+
+<para>
+Then, "run_test" is called with a different "seed" so that we get
+different data and noise realizations.
+</para>
+<programlisting>
+ 91          (s,e)=run_test(f,Kb,bitspersymbol,K,channel,modulation,dimensionality,tot_constellation,N0,-long(666+i)) # run experiment with different seed to get different data and noise realizations
+</programlisting>
+
+
+
+<para>
+Let us examine now the "run_test" function. 
+First we set up the transmitter blocks.
+We generate a packet of K random symbols and add a head and a tail of L zeros, 
+L being the channel length. This is sufficient to drive the initial and final states to 0.
+</para>
+<programlisting>
+ 14      L = len(channel)
+ 15  
+ 16      # TX
+ 17      # this for loop is TOO slow in python!!!
+ 18      packet = [0]*(K+2*L)
+ 19      random.seed(seed)
+ 20      for i in range(len(packet)):
+ 21          packet[i] = random.randint(0, 2**bitspersymbol - 1) # random symbols
+ 22      for i in range(L): # first/last L symbols set to 0
+ 23          packet[i] = 0
+ 24          packet[len(packet)-i-1] = 0
+ 25      src = gr.vector_source_s(packet,False)
+ 26      mod = gr.chunks_to_symbols_sf(modulation[1],modulation[0])
+</programlisting>
+
+
+<para>
+The modulated symbols are filtered by the ISI channel and AWGN with appropriate noise variance is added.
+</para>
+<programlisting>
+ 28      # CHANNEL
+ 29      isi = gr.fir_filter_fff(1,channel)
+ 30      add = gr.add_ff()
+ 31      noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
+</programlisting>
+
+
+
+<para>
+Since the output alphabet of the equivalent FSM is quite large (1024) we chose to utilize the combined metrics calculator and Viterbi algorithm block.
+also note that the first L observations are irrelevant and tus can be skipped.
+</para>
+<programlisting>
+ 33      # RX
+ 34      skip = gr.skiphead(gr.sizeof_float, L) # skip the first L samples since you know they are coming from the L zero symbols
+ 35      #metrics = trellis.metrics_f(f.O(),dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
+ 36      #va = trellis.viterbi_s(f,K+L,0,0) # Put -1 if the Initial/Final states are not set.
+ 37      va = trellis.viterbi_combined_s(f,K+L,0,0,dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # using viterbi_combined_s instead of metrics_f/viterbi_s allows larger packet lengths because metrics_f is complaining for not being able to allocate large buffers. This is due to the large f.O() in this application...                                    
+ 38      dst = gr.vector_sink_s()
+</programlisting>
+
+<para>
+Now the VA can run once it is supplied by the initial and final states.
+In this example both the initial and final states are set to 0.
+The VA outputs the estimates of the input symbols which
+are then compared with the transmitted sequence.
+</para>
+<programlisting>
+ 49      data = dst.data()
+ 50      ntotal = len(data) - L
+ 51      nright=0
+ 52      for i in range(ntotal):
+ 53          if packet[i+L]==data[i]:
+ 54              nright=nright+1
+ 55          #else:
+ 56              #print &quot;Error in &quot;, i
+</programlisting>
+
+
+<para>
+The function returns the number of symbols and the number of symbols in error. Observe that this way the estimated error rate refers to 
+2-bit-symbol error rate.
+</para>
+<programlisting>
+ 58      return (ntotal,ntotal-nright)
+</programlisting>
+
 
 
 </sect1>
 
 
+
+
+
+
+<!--=====================================================-->
+<sect1 id="turbo"><title>Support for Concatenated Coding and Turbo Decoding</title>
+
+<para>
+To do...
+</para>
+
+
+
+
+</sect1>
+
+
+
+
+
+
+
+
 <!--====================n================================-->
 <sect1 id="future"><title>Future Work</title>
 
diff --git a/gr-trellis/doc/test_viterbi_equalization1.py b/gr-trellis/doc/test_viterbi_equalization1.py
new file mode 100755
index 000000000..57d617aa0
--- /dev/null
+++ b/gr-trellis/doc/test_viterbi_equalization1.py
@@ -0,0 +1,102 @@
+#!/usr/bin/env python
+
+from gnuradio import gr
+from gnuradio import audio
+from gnuradio import trellis
+from gnuradio import eng_notation
+import math
+import sys
+import random
+import fsm_utils
+
+def run_test (f,Kb,bitspersymbol,K,channel,modulation,dimensionality,tot_constellation,N0,seed):
+    fg = gr.flow_graph ()
+    L = len(channel)
+
+    # TX
+    # this for loop is TOO slow in python!!!
+    packet = [0]*(K+2*L)
+    random.seed(seed)
+    for i in range(len(packet)):
+        packet[i] = random.randint(0, 2**bitspersymbol - 1) # random symbols
+    for i in range(L): # first/last L symbols set to 0
+        packet[i] = 0
+        packet[len(packet)-i-1] = 0
+    src = gr.vector_source_s(packet,False)
+    mod = gr.chunks_to_symbols_sf(modulation[1],modulation[0])
+
+    # CHANNEL
+    isi = gr.fir_filter_fff(1,channel)
+    add = gr.add_ff()
+    noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
+    
+    # RX
+    skip = gr.skiphead(gr.sizeof_float, L) # skip the first L samples since you know they are coming from the L zero symbols
+    #metrics = trellis.metrics_f(f.O(),dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
+    #va = trellis.viterbi_s(f,K+L,0,0) # Put -1 if the Initial/Final states are not set.
+    va = trellis.viterbi_combined_s(f,K+L,0,0,dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # using viterbi_combined_s instead of metrics_f/viterbi_s allows larger packet lengths because metrics_f is complaining for not being able to allocate large buffers. This is due to the large f.O() in this application...
+    dst = gr.vector_sink_s()
+
+    fg.connect (src,mod)
+    fg.connect (mod,isi,(add,0))
+    fg.connect (noise,(add,1))
+    #fg.connect (add,metrics)
+    #fg.connect (metrics,va,dst)
+    fg.connect (add,skip,va,dst)
+
+    fg.run()
+
+    data = dst.data() 
+    ntotal = len(data) - L
+    nright=0
+    for i in range(ntotal):
+        if packet[i+L]==data[i]:
+            nright=nright+1
+        #else:
+            #print "Error in ", i
+    
+    return (ntotal,ntotal-nright)
+
+
+def main(args):
+    nargs = len (args)
+    if nargs == 2:
+        esn0_db=float(args[0])
+        rep=int(args[1])
+    else:
+        sys.stderr.write ('usage: test_viterbi_equalization1.py Es/No_db  repetitions\n')
+        sys.exit (1)
+
+    # system parameters
+    Kb=2048  # packet size in bits
+    modulation = fsm_utils.pam4 # see fsm_utlis.py for available predefined modulations
+    channel = fsm_utils.c_channel # see fsm_utlis.py for available predefined test channels
+    f=trellis.fsm(len(modulation[1]),len(channel)) # generate the FSM automatically
+    bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
+    K=Kb/bitspersymbol # packet size in trellis steps
+
+    tot_channel = fsm_utils.make_isi_lookup(modulation,channel,True) # generate the lookup table (normalize energy to 1)
+    dimensionality = tot_channel[0]
+    tot_constellation = tot_channel[1]
+    N0=pow(10.0,-esn0_db/10.0); # noise variance
+    if len(tot_constellation)/dimensionality != f.O():
+        sys.stderr.write ('Incompatible FSM output cardinality and lookup table size.\n')
+        sys.exit (1)
+
+    tot_s=0 # total number of transmitted shorts
+    terr_s=0 # total number of shorts in error
+    terr_p=0 # total number of packets in error
+
+    for i in range(rep):
+        (s,e)=run_test(f,Kb,bitspersymbol,K,channel,modulation,dimensionality,tot_constellation,N0,-long(666+i)) # run experiment with different seed to get different data and noise realizations
+        tot_s=tot_s+s
+        terr_s=terr_s+e
+        terr_p=terr_p+(terr_s!=0)
+        if ((i+1)%100==0) : # display progress
+            print i+1,terr_p, '%.2e' % ((1.0*terr_p)/(i+1)),tot_s,terr_s, '%.2e' % ((1.0*terr_s)/tot_s)
+    # estimate of the (short or symbol) error rate
+    print rep,terr_p, '%.2e' % ((1.0*terr_p)/(i+1)),tot_s,terr_s, '%.2e' % ((1.0*terr_s)/tot_s)
+
+
+if __name__ == '__main__':
+    main (sys.argv[1:])
diff --git a/gr-trellis/doc/test_viterbi_equalization1.py.xml b/gr-trellis/doc/test_viterbi_equalization1.py.xml
new file mode 100644
index 000000000..cb13772fc
--- /dev/null
+++ b/gr-trellis/doc/test_viterbi_equalization1.py.xml
@@ -0,0 +1,105 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<programlisting>
+  1  #!/usr/bin/env python
+  2  
+  3  from gnuradio import gr
+  4  from gnuradio import audio
+  5  from gnuradio import trellis
+  6  from gnuradio import eng_notation
+  7  import math
+  8  import sys
+  9  import random
+ 10  import fsm_utils
+ 11  
+ 12  def run_test (f,Kb,bitspersymbol,K,channel,modulation,dimensionality,tot_constellation,N0,seed):
+ 13      fg = gr.flow_graph ()
+ 14      L = len(channel)
+ 15  
+ 16      # TX
+ 17      # this for loop is TOO slow in python!!!
+ 18      packet = [0]*(K+2*L)
+ 19      random.seed(seed)
+ 20      for i in range(len(packet)):
+ 21          packet[i] = random.randint(0, 2**bitspersymbol - 1) # random symbols
+ 22      for i in range(L): # first/last L symbols set to 0
+ 23          packet[i] = 0
+ 24          packet[len(packet)-i-1] = 0
+ 25      src = gr.vector_source_s(packet,False)
+ 26      mod = gr.chunks_to_symbols_sf(modulation[1],modulation[0])
+ 27  
+ 28      # CHANNEL
+ 29      isi = gr.fir_filter_fff(1,channel)
+ 30      add = gr.add_ff()
+ 31      noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
+ 32      
+ 33      # RX
+ 34      skip = gr.skiphead(gr.sizeof_float, L) # skip the first L samples since you know they are coming from the L zero symbols
+ 35      #metrics = trellis.metrics_f(f.O(),dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
+ 36      #va = trellis.viterbi_s(f,K+L,-1,0) # Put -1 if the Initial/Final states are not set.
+ 37      va = trellis.viterbi_combined_s(f,K+L,0,0,dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # using viterbi_combined_s instead of metrics_f/viterbi_s allows larger packet lengths because metrics_f is complaining for not being able to allocate large buffers. This is due to the large f.O() in this application...
+ 38      dst = gr.vector_sink_s()
+ 39  
+ 40      fg.connect (src,mod)
+ 41      fg.connect (mod,isi,(add,0))
+ 42      fg.connect (noise,(add,1))
+ 43      #fg.connect (add,metrics)
+ 44      #fg.connect (metrics,va,dst)
+ 45      fg.connect (add,skip,va,dst)
+ 46  
+ 47      fg.run()
+ 48  
+ 49      data = dst.data() 
+ 50      ntotal = len(data) - L
+ 51      nright=0
+ 52      for i in range(ntotal):
+ 53          if packet[i+L]==data[i]:
+ 54              nright=nright+1
+ 55          #else:
+ 56              #print &quot;Error in &quot;, i
+ 57      
+ 58      return (ntotal,ntotal-nright)
+ 59  
+ 60  
+ 61  def main(args):
+ 62      nargs = len (args)
+ 63      if nargs == 2:
+ 64          esn0_db=float(args[0])
+ 65          rep=int(args[1])
+ 66      else:
+ 67          sys.stderr.write (&apos;usage: test_viterbi_equalization1.py Es/No_db  repetitions\n&apos;)
+ 68          sys.exit (1)
+ 69  
+ 70      # system parameters
+ 71      Kb=2048  # packet size in bits
+ 72      modulation = fsm_utils.pam4 # see fsm_utlis.py for available predefined modulations
+ 73      channel = fsm_utils.c_channel # see fsm_utlis.py for available predefined test channels
+ 74      f=trellis.fsm(len(modulation[1]),len(channel)) # generate the FSM automatically
+ 75      bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
+ 76      K=Kb/bitspersymbol # packet size in trellis steps
+ 77  
+ 78      tot_channel = fsm_utils.make_isi_lookup(modulation,channel,True) # generate the lookup table (normalize energy to 1)
+ 79      dimensionality = tot_channel[0]
+ 80      tot_constellation = tot_channel[1]
+ 81      N0=pow(10.0,-esn0_db/10.0); # noise variance
+ 82      if len(tot_constellation)/dimensionality != f.O():
+ 83          sys.stderr.write (&apos;Incompatible FSM output cardinality and lookup table size.\n&apos;)
+ 84          sys.exit (1)
+ 85  
+ 86      tot_s=0 # total number of transmitted shorts
+ 87      terr_s=0 # total number of shorts in error
+ 88      terr_p=0 # total number of packets in error
+ 89  
+ 90      for i in range(rep):
+ 91          (s,e)=run_test(f,Kb,bitspersymbol,K,channel,modulation,dimensionality,tot_constellation,N0,-long(666+i)) # run experiment with different seed to get different data and noise realizations
+ 92          tot_s=tot_s+s
+ 93          terr_s=terr_s+e
+ 94          terr_p=terr_p+(terr_s!=0)
+ 95          if ((i+1)%100==0) : # display progress
+ 96              print i+1,terr_p, &apos;%.2e&apos; % ((1.0*terr_p)/(i+1)),tot_s,terr_s, &apos;%.2e&apos; % ((1.0*terr_s)/tot_s)
+ 97      # estimate of the (short or symbol) error rate
+ 98      print rep,terr_p, &apos;%.2e&apos; % ((1.0*terr_p)/(i+1)),tot_s,terr_s, &apos;%.2e&apos; % ((1.0*terr_s)/tot_s)
+ 99  
+100  
+101  if __name__ == &apos;__main__&apos;:
+102      main (sys.argv[1:])
+</programlisting>