Added support for Continuous Phase Modulation in gr-trellis + an example

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

6 files changed, 285 insertions, 37 deletions

diff --git a/gr-trellis/src/examples/fsm_utils.py b/gr-trellis/src/examples/fsm_utils.py
index 1b011246c..ab7b4e946 100755
--- a/gr-trellis/src/examples/fsm_utils.py
+++ b/gr-trellis/src/examples/fsm_utils.py
@@ -25,6 +25,8 @@ import re
 import math
 import sys
 import operator
+import numpy
+import scipy.linalg
 
 from gnuradio import trellis
 
@@ -58,33 +60,6 @@ def base2dec(s,base):
     return num
 
 
-######################################################################
-# Generate a new FSM representing the concatenation of two FSMs
-######################################################################
-def fsm_concatenate(f1,f2):
-    if f1.O() > f2.I():
-        print "Not compatible FSMs\n"
-    I=f1.I()
-    S=f1.S()*f2.S() 
-    O=f2.O()
-    nsm=list([0]*I*S)
-    osm=list([0]*I*S)
-    for s1 in range(f1.S()):
-        for s2 in range(f2.S()):
-            for i in range(f1.I()):
-                ns1=f1.NS()[s1*f1.I()+i]
-                o1=f1.OS()[s1*f1.I()+i]
-                ns2=f2.NS()[s2*f2.I()+o1]
-                o2=f2.OS()[s2*f2.I()+o1]
-
-                s=s1*f2.S()+s2
-                ns=ns1*f2.S()+ns2
-                nsm[s*I+i]=ns
-                osm[s*I+i]=o2
-
-
-    f=trellis.fsm(I,S,O,nsm,osm)
-    return f
 
 ######################################################################
 # Generate a new FSM representing n stages through the original FSM
@@ -143,6 +118,76 @@ def make_isi_lookup(mod,channel,normalize):
     return (1,lookup)
 
 
+
+
+
+
+######################################################################
+# Automatically generate the signals appropriate for CPM
+# decomposition. 
+# This decomposition is based on the paper by B. Rimoldi
+# "A decomposition approach to CPM", IEEE Trans. Info Theory, March 1988
+# See also my own notes at http://www.eecs.umich.edu/~anastas/docs/cpm.pdf
+######################################################################
+def make_cpm_signals(K,P,M,L,q,frac):
+
+    Q=numpy.size(q)/L
+    h=(1.0*K)/P
+    f0=-h*(M-1)/2
+    dt=0.0; # maybe start at t=0.5
+    t=(dt+numpy.arange(0,Q))/Q
+    qq=numpy.zeros(Q)
+    for m in range(L):
+       qq=qq + q[m*Q:m*Q+Q]
+    w=math.pi*h*(M-1)*t-2*math.pi*h*(M-1)*qq+math.pi*h*(L-1)*(M-1)
+    
+    X=(M**L)*P
+    PSI=numpy.empty((X,Q))
+    for x in range(X):
+       xv=dec2base(x/P,M,L)
+       xv=numpy.append(xv, x%P)
+       qq1=numpy.zeros(Q)
+       for m in range(L):
+          qq1=qq1+xv[m]*q[m*Q:m*Q+Q]
+       psi=2*math.pi*h*xv[-1]+4*math.pi*h*qq1+w
+       #print psi
+       PSI[x]=psi
+    PSI = numpy.transpose(PSI)
+    SS=numpy.exp(1j*PSI) # contains all signals as columns
+    #print SS
+   
+
+    # Now we need to orthogonalize the signals 
+    F = scipy.linalg.orth(SS) # find an orthonormal basis for SS
+    #print numpy.dot(numpy.transpose(F.conjugate()),F) # check for orthonormality
+    S = numpy.dot(numpy.transpose(F.conjugate()),SS)
+    #print F
+    #print S
+
+    # We only want to keep those dimensions that contain most
+    # of the energy of the overall constellation (eg, frac=0.9 ==> 90%)
+    # evaluate mean energy in each dimension
+    E=numpy.sum(numpy.absolute(S)**2,axis=1)/Q
+    E=E/numpy.sum(E)
+    #print E
+    Es = -numpy.sort(-E)
+    Esi = numpy.argsort(-E)
+    #print Es
+    #print Esi
+    Ecum=numpy.cumsum(Es)
+    #print Ecum
+    v0=numpy.searchsorted(Ecum,frac)
+    N = v0+1
+    #print v0
+    #print Esi[0:v0+1]
+    Ff=numpy.transpose(numpy.transpose(F)[Esi[0:v0+1]])
+    #print Ff
+    Sf = S[Esi[0:v0+1]]
+    #print Sf
+    
+
+    return (f0,SS,S,F,Sf,Ff,N)
+    #return f0
     
 
 
@@ -194,20 +239,23 @@ c_channel = [0.227, 0.460, 0.688, 0.460, 0.227]
 if __name__ == '__main__':
     f1=trellis.fsm('fsm_files/awgn1o2_4.fsm')
     #f2=trellis.fsm('fsm_files/awgn2o3_4.fsm')
-    print f1.I(), f1.S(), f1.O()
-    print f1.NS()
-    print f1.OS()
+    #print f1.I(), f1.S(), f1.O()
+    #print f1.NS()
+    #print f1.OS()
     #print f2.I(), f2.S(), f2.O()
     #print f2.NS()
     #print f2.OS()
     ##f1.write_trellis_svg('f1.svg',4)
     #f2.write_trellis_svg('f2.svg',4)
     #f=fsm_concatenate(f1,f2)
-    f=fsm_radix(f1,2)
+    #f=fsm_radix(f1,2)
 
-    print "----------\n"
-    print f.I(), f.S(), f.O()
-    print f.NS()
-    print f.OS()
+    #print "----------\n"
+    #print f.I(), f.S(), f.O()
+    #print f.NS()
+    #print f.OS()
     #f.write_trellis_svg('f.svg',4)
 
+    q=numpy.arange(0,8)/(2.0*8)
+    (f0,SS,S,F,Sf,Ff,N) = make_cpm_signals(1,2,2,1,q,0.99)
+
diff --git a/gr-trellis/src/examples/test_cpm.py b/gr-trellis/src/examples/test_cpm.py
new file mode 100755
index 000000000..ec432d4ff
--- /dev/null
+++ b/gr-trellis/src/examples/test_cpm.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+##################################################
+# Gnuradio Python Flow Graph
+# Title: CPM test
+# Author: Achilleas Anastasopoulos
+# Description: gnuradio flow graph
+# Generated: Thu Feb 19 23:16:23 2009
+##################################################
+
+from gnuradio import gr
+from gnuradio import trellis
+from gnuradio.gr import firdes
+from grc_gnuradio import blks2 as grc_blks2
+import math
+import numpy
+import scipy.stats
+import fsm_utils
+from gnuradio import trellis
+
+def run_test(seed,blocksize):
+        tb = gr.top_block()
+
+	##################################################
+	# Variables
+	##################################################
+	M = 2
+	K = 1
+	P = 2
+	h = (1.0*K)/P
+	L = 3
+	Q = 4
+        frac = 0.99
+        f = trellis.fsm(P,M,L)
+
+        # CPFSK signals
+        #p = numpy.ones(Q)/(2.0)
+        #q = numpy.cumsum(p)/(1.0*Q)
+
+        # GMSK signals
+        BT=0.3;
+        tt=numpy.arange(0,L*Q)/(1.0*Q)-L/2.0;
+        #print tt
+        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;
+        p=p/sum(p)*Q/2.0;
+        #print p
+        q=numpy.cumsum(p)/Q;
+        q=q/q[-1]/2.0;
+        #print q
+
+        (f0T,SS,S,F,Sf,Ff,N) = fsm_utils.make_cpm_signals(K,P,M,L,q,frac)
+        #print N
+        #print Ff
+        Ffa = numpy.insert(Ff,Q,numpy.zeros(N),axis=0)
+        #print Ffa
+        MF = numpy.fliplr(numpy.transpose(Ffa))
+        #print MF
+        E = numpy.sum(numpy.abs(Sf)**2,axis=0)
+        Es = numpy.sum(E)/f.O()
+        #print Es
+
+        constellation = numpy.reshape(numpy.transpose(Sf),N*f.O())
+        #print Ff
+        #print Sf
+        #print constellation
+        #print numpy.max(numpy.abs(SS - numpy.dot(Ff , Sf)))
+
+	EsN0_db = 10.0
+        N0 =  Es * 10.0**(-(1.0*EsN0_db)/10.0)
+        #N0 = 0.0
+        #print N0
+        head = 4
+        tail = 4
+        numpy.random.seed(seed*666)
+        data = numpy.random.randint(0, M, head+blocksize+tail+1)
+        #data = numpy.zeros(blocksize+1+head+tail,'int')
+        for i in range(head):
+            data[i]=0
+        for i in range(tail+1):
+            data[-i]=0
+      
+
+
+	##################################################
+	# Blocks
+	##################################################
+	random_source_x_0 = gr.vector_source_b(data, False)
+	gr_chunks_to_symbols_xx_0 = gr.chunks_to_symbols_bf((-1, 1), 1)
+	gr_interp_fir_filter_xxx_0 = gr.interp_fir_filter_fff(Q, p)
+	gr_frequency_modulator_fc_0 = gr.frequency_modulator_fc(2*math.pi*h*(1.0/Q))
+
+	gr_add_vxx_0 = gr.add_vcc(1)
+	gr_noise_source_x_0 = gr.noise_source_c(gr.GR_GAUSSIAN, (N0/2.0)**0.5, -long(seed))
+
+	gr_multiply_vxx_0 = gr.multiply_vcc(1)
+	gr_sig_source_x_0 = gr.sig_source_c(Q, gr.GR_COS_WAVE, -f0T, 1, 0)
+        # only works for N=2, do it manually for N>2...
+	gr_fir_filter_xxx_0_0 = gr.fir_filter_ccc(Q, MF[0].conjugate())
+	gr_fir_filter_xxx_0_0_0 = gr.fir_filter_ccc(Q, MF[1].conjugate())
+	gr_streams_to_stream_0 = gr.streams_to_stream(gr.sizeof_gr_complex*1, N)
+	gr_skiphead_0 = gr.skiphead(gr.sizeof_gr_complex*1, N*(1+0))
+	viterbi = trellis.viterbi_combined_cb(f, head+blocksize+tail, 0, -1, N, constellation, trellis.TRELLIS_EUCLIDEAN)
+
+        gr_vector_sink_x_0 = gr.vector_sink_b()
+
+	##################################################
+	# Connections
+	##################################################
+	tb.connect((random_source_x_0, 0), (gr_chunks_to_symbols_xx_0, 0))
+	tb.connect((gr_chunks_to_symbols_xx_0, 0), (gr_interp_fir_filter_xxx_0, 0))
+	tb.connect((gr_interp_fir_filter_xxx_0, 0), (gr_frequency_modulator_fc_0, 0))
+	tb.connect((gr_frequency_modulator_fc_0, 0), (gr_add_vxx_0, 0))
+	tb.connect((gr_noise_source_x_0, 0), (gr_add_vxx_0, 1))
+	tb.connect((gr_add_vxx_0, 0), (gr_multiply_vxx_0, 0))
+	tb.connect((gr_sig_source_x_0, 0), (gr_multiply_vxx_0, 1))
+	tb.connect((gr_multiply_vxx_0, 0), (gr_fir_filter_xxx_0_0, 0))
+	tb.connect((gr_multiply_vxx_0, 0), (gr_fir_filter_xxx_0_0_0, 0))
+	tb.connect((gr_fir_filter_xxx_0_0, 0), (gr_streams_to_stream_0, 0))
+	tb.connect((gr_fir_filter_xxx_0_0_0, 0), (gr_streams_to_stream_0, 1))
+	tb.connect((gr_streams_to_stream_0, 0), (gr_skiphead_0, 0))
+	tb.connect((gr_skiphead_0, 0), (viterbi, 0))
+	tb.connect((viterbi, 0), (gr_vector_sink_x_0, 0))
+        
+
+        tb.run()
+        dataest = gr_vector_sink_x_0.data()
+        #print data
+        #print numpy.array(dataest)
+        perr = 0
+        err = 0
+        for i in range(blocksize):
+          if data[head+i] != dataest[head+i]:
+            #print i
+            err += 1
+        if err != 0 :
+          perr = 1
+        return (err,perr)
+
+if __name__ == '__main__':
+        blocksize = 1000
+        ss=0
+        ee=0
+        for i in range(10000):
+            (s,e) = run_test(i,blocksize)
+            ss += s
+            ee += e
+            if (i+1) % 100 == 0:
+                print i+1,ss,ee,(1.0*ss)/(i+1)/(1.0*blocksize),(1.0*ee)/(i+1)
+        print i+1,ss,ee,(1.0*ss)/(i+1)/(1.0*blocksize),(1.0*ee)/(i+1)
diff --git a/gr-trellis/src/lib/fsm.cc b/gr-trellis/src/lib/fsm.cc
index 50598a94c..550674ad5 100644
--- a/gr-trellis/src/lib/fsm.cc
+++ b/gr-trellis/src/lib/fsm.cc
@@ -170,7 +170,7 @@ fsm::fsm(int k, int n, const std::vector<int> &G)
 
 
   for(int s=0;s<d_S;s++) {
-    dec2bases(s,bases_x,sx); // split s into k values, each representing on of the k shift registers
+    dec2bases(s,bases_x,sx); // split s into k values, each representing one of the k shift registers
 //printf("state = %d \nstates = ",s);
 //for(int j=0;j<sx.size();j++) printf("%d ",sx[j]); printf("\n");
     for(int i=0;i<d_I;i++) {
@@ -239,6 +239,53 @@ fsm::fsm(int mod_size, int ch_length)
 }
 
 
+
+
+//######################################################################
+//# Automatically generate an FSM specification describing the 
+//# the trellis for a CPM with h=K/P (relatively prime), 
+//# alphabet size M, and frequency pulse duration L symbols
+//#
+//# This FSM is based on the paper by B. Rimoldi
+//# "A decomposition approach to CPM", IEEE Trans. Info Theory, March 1988
+//# See also my own notes at http://www.eecs.umich.edu/~anastas/docs/cpm.pdf
+//######################################################################
+fsm::fsm(int P, int M, int L)
+{
+  d_I=M;
+  d_S=(int)(pow(1.0*M,1.0*L-1)+0.5)*P;
+  d_O=(int)(pow(1.0*M,1.0*L)+0.5)*P;
+
+  d_NS.resize(d_I*d_S);
+  d_OS.resize(d_I*d_S);
+  int nv;
+  for(int s=0;s<d_S;s++) {
+    for(int i=0;i<d_I;i++) {
+      int s1=s/P;
+      int v=s%P;
+      int ns1= (i*(int)(pow(1.0*M,1.0*(L-1))+0.5)+s1)/M;
+      if (L==1)
+        nv=(i+v)%P;
+      else
+        nv=(s1%M+v)%P;
+      d_NS[s*d_I+i] = ns1*P+nv;
+      d_OS[s*d_I+i] = i*d_S+s;
+    }
+  }
+
+  generate_PS_PI();
+  generate_TM();
+}
+
+
+
+
+
+
+
+
+
+
 //######################################################################
 //# Automatically generate an FSM specification describing the 
 //# the joint trellis of fsm1 and fsm2
@@ -419,7 +466,7 @@ void fsm::write_trellis_svg( std::string filename ,int number_stages)
 
 
 //######################################################################
-//# Write trellis specification to a text files,
+//# Write trellis specification to a text file,
 //# in the same format used when reading FSM files
 //######################################################################
 void fsm::write_fsm_txt(std::string filename)
diff --git a/gr-trellis/src/lib/fsm.h b/gr-trellis/src/lib/fsm.h
index e9ebb91ec..b91dc2ab1 100644
--- a/gr-trellis/src/lib/fsm.h
+++ b/gr-trellis/src/lib/fsm.h
@@ -50,6 +50,7 @@ public:
   fsm(const char *name);
   fsm(int k, int n, const std::vector<int> &G);
   fsm(int mod_size, int ch_length);
+  fsm(int P, int M, int L);
   fsm(const fsm &FSM1, const fsm &FSM2);
   int I () const { return d_I; }
   int S () const { return d_S; }
diff --git a/gr-trellis/src/lib/fsm.i b/gr-trellis/src/lib/fsm.i
index ac00d40e3..e9db804a0 100644
--- a/gr-trellis/src/lib/fsm.i
+++ b/gr-trellis/src/lib/fsm.i
@@ -40,6 +40,7 @@ public:
   fsm(const char *name);
   fsm(int k, int n, const std::vector<int> &G);
   fsm(int mod_size, int ch_length);
+  fsm(int P, int M, int L);
   fsm(const fsm &FSM1, const fsm &FSM2);
   int I () const { return d_I; }
   int S () const { return d_S; }
diff --git a/gr-trellis/src/lib/trellis_calc_metric.cc b/gr-trellis/src/lib/trellis_calc_metric.cc
index 676f53fd7..0d03fd1a7 100644
--- a/gr-trellis/src/lib/trellis_calc_metric.cc
+++ b/gr-trellis/src/lib/trellis_calc_metric.cc
@@ -31,6 +31,7 @@ void calc_metric(int O, int D, const std::vector<T> &TABLE, const T *in, float *
 {
   float minm = FLT_MAX;
   int minmi = 0;
+  
 
   switch (type){
   case TRELLIS_EUCLIDEAN:
@@ -213,6 +214,7 @@ void calc_metric(int O, int D, const std::vector<gr_complex> &TABLE, const gr_co
   float minm = FLT_MAX;
   int minmi = 0;
 
+
   switch (type){
   case TRELLIS_EUCLIDEAN:
     for(int o=0;o<O;o++) {
@@ -222,6 +224,7 @@ void calc_metric(int O, int D, const std::vector<gr_complex> &TABLE, const gr_co
         metric[o]+=s.real()*s.real()+s.imag()*s.imag();
       }
     }
+    break;
   case TRELLIS_HARD_SYMBOL:
     for(int o=0;o<O;o++) {
       metric[o]=0.0;