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

   1 #!/usr/bin/env python
   2 #
   3 # Copyright 2004 Free Software Foundation, Inc.
   4 #
   5 # This file is part of GNU Radio
   6 #
   7 # GNU Radio is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation; either version 3, or (at your option)
  10 # any later version.
  11 #
  12 # GNU Radio is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with GNU Radio; see the file COPYING.  If not, write to
  19 # the Free Software Foundation, Inc., 51 Franklin Street,
  20 # Boston, MA 02110-1301, USA.
  21 #
  22
  23
  24 import re
  25 import math
  26 import sys
  27 import operator
  28 import numpy
  29 import scipy.linalg
  30
  31 from gnuradio import trellis
  32
  33
  34
  35 ######################################################################
  36 # Decimal to any base conversion.
  37 # Convert 'num' to a list of 'l' numbers representing 'num'
  38 # to base 'base' (most significant symbol first).
  39 ######################################################################
  40 def dec2base(num,base,l):
  41     s=range(l)
  42     n=num
  43     for i in range(l):
  44         s[l-i-1]=n%base
  45         n=int(n/base)
  46     if n!=0:
  47         print 'Number ', num, ' requires more than ', l, 'digits.'
  48     return s
  49
  50
  51 ######################################################################
  52 # Conversion from any base to decimal.
  53 # Convert a list 's' of symbols to a decimal number
  54 # (most significant symbol first)
  55 ######################################################################
  56 def base2dec(s,base):
  57     num=0
  58     for i in range(len(s)):
  59         num=num*base+s[i]
  60     return num
  61
  62
  63
  64 ######################################################################
  65 # Generate a new FSM representing n stages through the original FSM
  66 ######################################################################
  67 def fsm_radix(f,n):
  68     I=f.I()**n
  69     S=f.S()
  70     O=f.O()**n
  71     nsm=list([0]*I*S)
  72     osm=list([0]*I*S)
  73     for s in range(f.S()):
  74         for i in range(I):
  75             ii=dec2base(i,f.I(),n)
  76             oo=list([0]*n)
  77             ns=s
  78             for k in range(n):
  79                 oo[k]=f.OS()[ns*f.I()+ii[k]]
  80                 ns=f.NS()[ns*f.I()+ii[k]]
  81
  82             nsm[s*I+i]=ns
  83             osm[s*I+i]=base2dec(oo,f.O())
  84
  85
  86     f=trellis.fsm(I,S,O,nsm,osm)
  87     return f
  88
  89
  90
  91
  92 ######################################################################
  93 # Automatically generate the lookup table that maps the FSM outputs
  94 # to channel inputs corresponding to a channel 'channel' and a modulation
  95 # 'mod'. Optional normalization of channel to unit energy.
  96 # This table is used by the 'metrics' block to translate
  97 # channel outputs to metrics for use with the Viterbi algorithm.
  98 # Limitations: currently supports only one-dimensional modulations.
  99 ######################################################################
 100 def make_isi_lookup(mod,channel,normalize):
 101     dim=mod[0]
 102     constellation = mod[1]
 103
 104     if normalize:
 105         p = 0
 106         for i in range(len(channel)):
 107             p = p + channel[i]**2
 108         for i in range(len(channel)):
 109             channel[i] = channel[i]/math.sqrt(p)
 110
 111     lookup=range(len(constellation)**len(channel))
 112     for o in range(len(constellation)**len(channel)):
 113         ss=dec2base(o,len(constellation),len(channel))
 114         ll=0
 115         for i in range(len(channel)):
 116             ll=ll+constellation[ss[i]]*channel[i]
 117         lookup[o]=ll
 118     return (1,lookup)
 119
 120
 121
 122
 123
 124
 125 ######################################################################
 126 # Automatically generate the signals appropriate for CPM
 127 # decomposition.
 128 # This decomposition is based on the paper by B. Rimoldi
 129 # "A decomposition approach to CPM", IEEE Trans. Info Theory, March 1988
 130 # See also my own notes at http://www.eecs.umich.edu/~anastas/docs/cpm.pdf
 131 ######################################################################
 132 def make_cpm_signals(K,P,M,L,q,frac):
 133
 134     Q=numpy.size(q)/L
 135     h=(1.0*K)/P
 136     f0=-h*(M-1)/2
 137     dt=0.0; # maybe start at t=0.5
 138     t=(dt+numpy.arange(0,Q))/Q
 139     qq=numpy.zeros(Q)
 140     for m in range(L):
 141        qq=qq + q[m*Q:m*Q+Q]
 142     w=math.pi*h*(M-1)*t-2*math.pi*h*(M-1)*qq+math.pi*h*(L-1)*(M-1)
 143
 144     X=(M**L)*P
 145     PSI=numpy.empty((X,Q))
 146     for x in range(X):
 147        xv=dec2base(x/P,M,L)
 148        xv=numpy.append(xv, x%P)
 149        qq1=numpy.zeros(Q)
 150        for m in range(L):
 151           qq1=qq1+xv[m]*q[m*Q:m*Q+Q]
 152        psi=2*math.pi*h*xv[-1]+4*math.pi*h*qq1+w
 153        #print psi
 154        PSI[x]=psi
 155     PSI = numpy.transpose(PSI)
 156     SS=numpy.exp(1j*PSI) # contains all signals as columns
 157     #print SS
 158
 159
 160     # Now we need to orthogonalize the signals
 161     F = scipy.linalg.orth(SS) # find an orthonormal basis for SS
 162     #print numpy.dot(numpy.transpose(F.conjugate()),F) # check for orthonormality
 163     S = numpy.dot(numpy.transpose(F.conjugate()),SS)
 164     #print F
 165     #print S
 166
 167     # We only want to keep those dimensions that contain most
 168     # of the energy of the overall constellation (eg, frac=0.9 ==> 90%)
 169     # evaluate mean energy in each dimension
 170     E=numpy.sum(numpy.absolute(S)**2,axis=1)/Q
 171     E=E/numpy.sum(E)
 172     #print E
 173     Es = -numpy.sort(-E)
 174     Esi = numpy.argsort(-E)
 175     #print Es
 176     #print Esi
 177     Ecum=numpy.cumsum(Es)
 178     #print Ecum
 179     v0=numpy.searchsorted(Ecum,frac)
 180     N = v0+1
 181     #print v0
 182     #print Esi[0:v0+1]
 183     Ff=numpy.transpose(numpy.transpose(F)[Esi[0:v0+1]])
 184     #print Ff
 185     Sf = S[Esi[0:v0+1]]
 186     #print Sf
 187
 188
 189     return (f0,SS,S,F,Sf,Ff,N)
 190     #return f0
 191
 192
 193
 194
 195 ######################################################################
 196 # A list of common modulations.
 197 # Format: (dimensionality,constellation)
 198 ######################################################################
 199 pam2 = (1,[-1, 1])
 200 pam4 = (1,[-3, -1, 3, 1])               # includes Gray mapping
 201 pam8 = (1,[-7, -5, -3, -1, 1, 3, 5, 7])
 202
 203 psk4=(2,[1, 0, \
 204          0, 1, \
 205          0, -1,\
 206         -1, 0])                         # includes Gray mapping
 207 psk8=(2,[math.cos(2*math.pi*0/8), math.sin(2*math.pi*0/8),  \
 208          math.cos(2*math.pi*1/8), math.sin(2*math.pi*1/8),  \
 209          math.cos(2*math.pi*2/8), math.sin(2*math.pi*2/8),  \
 210          math.cos(2*math.pi*3/8), math.sin(2*math.pi*3/8),  \
 211          math.cos(2*math.pi*4/8), math.sin(2*math.pi*4/8),  \
 212          math.cos(2*math.pi*5/8), math.sin(2*math.pi*5/8),  \
 213          math.cos(2*math.pi*6/8), math.sin(2*math.pi*6/8),  \
 214          math.cos(2*math.pi*7/8), math.sin(2*math.pi*7/8)])
 215
 216 orth2 = (2,[1, 0, \
 217             0, 1])
 218 orth4=(4,[1, 0, 0, 0, \
 219           0, 1, 0, 0, \
 220           0, 0, 1, 0, \
 221           0, 0, 0, 1])
 222
 223 ######################################################################
 224 # A list of channels to be tested
 225 ######################################################################
 226
 227 # C test channel (J. Proakis, Digital Communications, McGraw-Hill Inc., 2001)
 228 c_channel = [0.227, 0.460, 0.688, 0.460, 0.227]
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239 if __name__ == '__main__':
 240     f1=trellis.fsm('fsm_files/awgn1o2_4.fsm')
 241     #f2=trellis.fsm('fsm_files/awgn2o3_4.fsm')
 242     #print f1.I(), f1.S(), f1.O()
 243     #print f1.NS()
 244     #print f1.OS()
 245     #print f2.I(), f2.S(), f2.O()
 246     #print f2.NS()
 247     #print f2.OS()
 248     ##f1.write_trellis_svg('f1.svg',4)
 249     #f2.write_trellis_svg('f2.svg',4)
 250     #f=fsm_concatenate(f1,f2)
 251     #f=fsm_radix(f1,2)
 252
 253     #print "----------\n"
 254     #print f.I(), f.S(), f.O()
 255     #print f.NS()
 256     #print f.OS()
 257     #f.write_trellis_svg('f.svg',4)
 258
 259     q=numpy.arange(0,8)/(2.0*8)
 260     (f0,SS,S,F,Sf,Ff,N) = make_cpm_signals(1,2,2,1,q,0.99)
 261