Added examples for Viterbi Equalization in gnuradio-examples/python/channel-coding

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

2 files changed, 195 insertions, 0 deletions

diff --git a/gnuradio-examples/python/channel-coding/test_viterbi_equalization.py b/gnuradio-examples/python/channel-coding/test_viterbi_equalization.py
new file mode 100755
index 000000000..24545e3b0
--- /dev/null
+++ b/gnuradio-examples/python/channel-coding/test_viterbi_equalization.py
@@ -0,0 +1,94 @@
+#!/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 fsm_utils
+
+def run_test (f,Kb,bitspersymbol,K,dimensionality,tot_constellation,N0,seed):
+    fg = gr.flow_graph ()
+
+    # TX
+    src = gr.lfsr_32k_source_s()
+    src_head = gr.head (gr.sizeof_short,Kb/16) # packet size in shorts
+    s2fsmi = gr.packed_to_unpacked_ss(bitspersymbol,gr.GR_MSB_FIRST) # unpack shorts to symbols compatible with the FSM input cardinality
+    enc = trellis.encoder_ss(f,0) # initial state = 0
+    # essentially here we implement the combination of modulation and channel as a memoryless modulation (the memory induced by the channel is hidden in the FSM)
+    mod = gr.chunks_to_symbols_sf(tot_constellation,dimensionality)
+
+    # CHANNEL
+    add = gr.add_ff()
+    noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
+    
+    # RX
+    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,0,-1) # Put -1 if the Initial/Final states are not set.
+    fsmi2s = gr.unpacked_to_packed_ss(bitspersymbol,gr.GR_MSB_FIRST) # pack FSM input symbols to shorts
+    dst = gr.check_lfsr_32k_s(); 
+    
+    fg.connect (src,src_head,s2fsmi,enc,mod)
+    fg.connect (mod,(add,0))
+    fg.connect (noise,(add,1))
+    fg.connect (add,metrics)
+    fg.connect (metrics,va,fsmi2s,dst)
+    
+    fg.run()
+
+    ntotal = dst.ntotal ()
+    nright = dst.nright ()
+    runlength = dst.runlength ()
+    #print ntotal,nright,runlength 
+    
+    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_equalization.py Es/No_db  repetitions\n')
+        sys.exit (1)
+
+    # system parameters
+    Kb=128*16  # packet size in bits (multiple of 16)
+    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,dimensionality,tot_constellation,N0,-long(666+i)) # run experiment with different seed to get different 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 bit) 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/gnuradio-examples/python/channel-coding/test_viterbi_equalization1.py b/gnuradio-examples/python/channel-coding/test_viterbi_equalization1.py
new file mode 100755
index 000000000..133e550db
--- /dev/null
+++ b/gnuradio-examples/python/channel-coding/test_viterbi_equalization1.py
@@ -0,0 +1,101 @@
+#!/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)
+    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
+    #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+2*L,-1,0) # Put -1 if the Initial/Final states are not set. Better if we could skip the first L symbols and start with a 0 state... don't know how to "skip" samples in gnuradio
+    va = trellis.viterbi_combined_s(f,dimensionality,tot_constellation,K+2*L,-1,0,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,va,dst)
+
+    fg.run()
+
+    data = dst.data() 
+    ntotal = len(data) - 2*L
+    nright=0
+    for i in range(ntotal):
+        if packet[i+L]==data[i+L]:
+            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=128*16  # packet size in bits (multiple of 16)
+    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 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:])