git.gag.com Git - debian/gnuradio/blob - gr-trellis/src/examples/test_cpm.py

   1 #!/usr/bin/env python
   2 ##################################################
   3 # Gnuradio Python Flow Graph
   4 # Title: CPM test
   5 # Author: Achilleas Anastasopoulos
   6 # Description: gnuradio flow graph
   7 # Generated: Thu Feb 19 23:16:23 2009
   8 ##################################################
   9
  10 from gnuradio import gr
  11 from gnuradio import trellis
  12 from gnuradio.gr import firdes
  13 from grc_gnuradio import blks2 as grc_blks2
  14 import math
  15 import numpy
  16 import scipy.stats
  17 import fsm_utils
  18 from gnuradio import trellis
  19
  20 def run_test(seed,blocksize):
  21         tb = gr.top_block()
  22
  23         ##################################################
  24         # Variables
  25         ##################################################
  26         M = 2
  27         K = 1
  28         P = 2
  29         h = (1.0*K)/P
  30         L = 3
  31         Q = 4
  32         frac = 0.99
  33         f = trellis.fsm(P,M,L)
  34
  35         # CPFSK signals
  36         #p = numpy.ones(Q)/(2.0)
  37         #q = numpy.cumsum(p)/(1.0*Q)
  38
  39         # GMSK signals
  40         BT=0.3;
  41         tt=numpy.arange(0,L*Q)/(1.0*Q)-L/2.0;
  42         #print tt
  43         p=(0.5*scipy.stats.erfc(2*math.pi*BT*(tt-0.5)/math.sqrt(math.log(2.0))/math.sqrt(2.0))-0.5*scipy.stats.erfc(2*math.pi*BT*(tt+0.5)/math.sqrt(math.log(2.0))/math.sqrt(2.0)))/2.0;
  44         p=p/sum(p)*Q/2.0;
  45         #print p
  46         q=numpy.cumsum(p)/Q;
  47         q=q/q[-1]/2.0;
  48         #print q
  49
  50         (f0T,SS,S,F,Sf,Ff,N) = fsm_utils.make_cpm_signals(K,P,M,L,q,frac)
  51         #print N
  52         #print Ff
  53         Ffa = numpy.insert(Ff,Q,numpy.zeros(N),axis=0)
  54         #print Ffa
  55         MF = numpy.fliplr(numpy.transpose(Ffa))
  56         #print MF
  57         E = numpy.sum(numpy.abs(Sf)**2,axis=0)
  58         Es = numpy.sum(E)/f.O()
  59         #print Es
  60
  61         constellation = numpy.reshape(numpy.transpose(Sf),N*f.O())
  62         #print Ff
  63         #print Sf
  64         #print constellation
  65         #print numpy.max(numpy.abs(SS - numpy.dot(Ff , Sf)))
  66
  67         EsN0_db = 10.0
  68         N0 =  Es * 10.0**(-(1.0*EsN0_db)/10.0)
  69         #N0 = 0.0
  70         #print N0
  71         head = 4
  72         tail = 4
  73         numpy.random.seed(seed*666)
  74         data = numpy.random.randint(0, M, head+blocksize+tail+1)
  75         #data = numpy.zeros(blocksize+1+head+tail,'int')
  76         for i in range(head):
  77             data[i]=0
  78         for i in range(tail+1):
  79             data[-i]=0
  80
  81
  82
  83         ##################################################
  84         # Blocks
  85         ##################################################
  86         random_source_x_0 = gr.vector_source_b(data, False)
  87         gr_chunks_to_symbols_xx_0 = gr.chunks_to_symbols_bf((-1, 1), 1)
  88         gr_interp_fir_filter_xxx_0 = gr.interp_fir_filter_fff(Q, p)
  89         gr_frequency_modulator_fc_0 = gr.frequency_modulator_fc(2*math.pi*h*(1.0/Q))
  90
  91         gr_add_vxx_0 = gr.add_vcc(1)
  92         gr_noise_source_x_0 = gr.noise_source_c(gr.GR_GAUSSIAN, (N0/2.0)**0.5, -long(seed))
  93
  94         gr_multiply_vxx_0 = gr.multiply_vcc(1)
  95         gr_sig_source_x_0 = gr.sig_source_c(Q, gr.GR_COS_WAVE, -f0T, 1, 0)
  96         # only works for N=2, do it manually for N>2...
  97         gr_fir_filter_xxx_0_0 = gr.fir_filter_ccc(Q, MF[0].conjugate())
  98         gr_fir_filter_xxx_0_0_0 = gr.fir_filter_ccc(Q, MF[1].conjugate())
  99         gr_streams_to_stream_0 = gr.streams_to_stream(gr.sizeof_gr_complex*1, N)
 100         gr_skiphead_0 = gr.skiphead(gr.sizeof_gr_complex*1, N*(1+0))
 101         viterbi = trellis.viterbi_combined_cb(f, head+blocksize+tail, 0, -1, N, constellation, trellis.TRELLIS_EUCLIDEAN)
 102
 103         gr_vector_sink_x_0 = gr.vector_sink_b()
 104
 105         ##################################################
 106         # Connections
 107         ##################################################
 108         tb.connect((random_source_x_0, 0), (gr_chunks_to_symbols_xx_0, 0))
 109         tb.connect((gr_chunks_to_symbols_xx_0, 0), (gr_interp_fir_filter_xxx_0, 0))
 110         tb.connect((gr_interp_fir_filter_xxx_0, 0), (gr_frequency_modulator_fc_0, 0))
 111         tb.connect((gr_frequency_modulator_fc_0, 0), (gr_add_vxx_0, 0))
 112         tb.connect((gr_noise_source_x_0, 0), (gr_add_vxx_0, 1))
 113         tb.connect((gr_add_vxx_0, 0), (gr_multiply_vxx_0, 0))
 114         tb.connect((gr_sig_source_x_0, 0), (gr_multiply_vxx_0, 1))
 115         tb.connect((gr_multiply_vxx_0, 0), (gr_fir_filter_xxx_0_0, 0))
 116         tb.connect((gr_multiply_vxx_0, 0), (gr_fir_filter_xxx_0_0_0, 0))
 117         tb.connect((gr_fir_filter_xxx_0_0, 0), (gr_streams_to_stream_0, 0))
 118         tb.connect((gr_fir_filter_xxx_0_0_0, 0), (gr_streams_to_stream_0, 1))
 119         tb.connect((gr_streams_to_stream_0, 0), (gr_skiphead_0, 0))
 120         tb.connect((gr_skiphead_0, 0), (viterbi, 0))
 121         tb.connect((viterbi, 0), (gr_vector_sink_x_0, 0))
 122
 123
 124         tb.run()
 125         dataest = gr_vector_sink_x_0.data()
 126         #print data
 127         #print numpy.array(dataest)
 128         perr = 0
 129         err = 0
 130         for i in range(blocksize):
 131           if data[head+i] != dataest[head+i]:
 132             #print i
 133             err += 1
 134         if err != 0 :
 135           perr = 1
 136         return (err,perr)
 137
 138 if __name__ == '__main__':
 139         blocksize = 1000
 140         ss=0
 141         ee=0
 142         for i in range(10000):
 143             (s,e) = run_test(i,blocksize)
 144             ss += s
 145             ee += e
 146             if (i+1) % 100 == 0:
 147                 print i+1,ss,ee,(1.0*ss)/(i+1)/(1.0*blocksize),(1.0*ee)/(i+1)
 148         print i+1,ss,ee,(1.0*ss)/(i+1)/(1.0*blocksize),(1.0*ee)/(i+1)