Several enhancements to gr-trellis and gnuradio-examples/python/channel-coding:
[debian/gnuradio] / gr-trellis / src / lib / fsm.cc
index 05d1589fe6aea40ed15e67e5499ddb07cfb32eae..f343088343fc31c6791ede0531a07330fe48b163 100644 (file)
 #include <cstdio>\r
 #include <stdexcept>\r
 #include <cmath>\r
+#include "base.h"\r
 #include "fsm.h"\r
 \r
+\r
 fsm::fsm()\r
 {\r
   d_I=0;\r
@@ -34,6 +36,8 @@ fsm::fsm()
   d_OS.resize(0);\r
   d_PS.resize(0);\r
   d_PI.resize(0);\r
+  d_TMi.resize(0);\r
+  d_TMl.resize(0);\r
 }\r
 \r
 fsm::fsm(const fsm &FSM)\r
@@ -45,28 +49,20 @@ fsm::fsm(const fsm &FSM)
   d_OS=FSM.OS();\r
   d_PS=FSM.PS();\r
   d_PI=FSM.PI();\r
+  d_TMi=FSM.TMi();\r
+  d_TMl=FSM.TMl();\r
 }\r
 \r
-fsm::fsm(const int I, const int S, const int O, const std::vector<int> &NS, const std::vector<int> &OS)\r
+fsm::fsm(int I, int S, int O, const std::vector<int> &NS, const std::vector<int> &OS)\r
 {\r
   d_I=I;\r
   d_S=S;\r
   d_O=O;\r
   d_NS=NS;\r
   d_OS=OS;\r
-  d_PS.resize(d_I*d_S);\r
-  d_PI.resize(d_I*d_S);\r
-  \r
-  // generate the PS, PI tables for later use\r
-  for(int i=0;i<d_S;i++) {\r
-    int j=0;\r
-    for(int ii=0;ii<d_S;ii++) for(int jj=0;jj<d_I;jj++) {\r
-      if(d_NS[ii*d_I+jj]!=i) continue;\r
-      d_PS[i*d_I+j]=ii;\r
-      d_PI[i*d_I+j]=jj;\r
-      j++;\r
-    }\r
-  }\r
\r
+  generate_PS_PI();\r
+  generate_TM();\r
 }\r
 \r
 //######################################################################\r
@@ -84,14 +80,12 @@ fsm::fsm(const char *name)
   FILE *fsmfile;\r
 \r
   if((fsmfile=fopen(name,"r"))==NULL) \r
-    throw std::runtime_error ("file open error in fsm()");\r
+    throw std::runtime_error ("fsm::fsm(const char *name): file open error\n");\r
     //printf("file open error in fsm()\n");\r
   \r
   fscanf(fsmfile,"%d %d %d\n",&d_I,&d_S,&d_O);\r
   d_NS.resize(d_I*d_S);\r
   d_OS.resize(d_I*d_S);\r
-  d_PS.resize(d_I*d_S);\r
-  d_PI.resize(d_I*d_S);\r
 \r
   for(int i=0;i<d_S;i++) {\r
     for(int j=0;j<d_I;j++) fscanf(fsmfile,"%d",&(d_NS[i*d_I+j]));\r
@@ -99,25 +93,128 @@ fsm::fsm(const char *name)
   for(int i=0;i<d_S;i++) {\r
     for(int j=0;j<d_I;j++) fscanf(fsmfile,"%d",&(d_OS[i*d_I+j]));\r
   }\r
\r
+  generate_PS_PI();\r
+  generate_TM();\r
+}\r
+\r
+\r
+\r
+\r
+//######################################################################\r
+//# Automatically generate the FSM from the generator matrix\r
+//# of a (n,k) binary convolutional code\r
+//######################################################################\r
+fsm::fsm(int k, int n, const std::vector<int> &G)\r
+{\r
+\r
+  // calculate maximum memory requirements for each input stream\r
+  std::vector<int> max_mem_x(k,-1);\r
+  int max_mem = -1;\r
+  for(int i=0;i<k;i++) {\r
+    for(int j=0;j<n;j++) {\r
+      int mem = -1;\r
+      if(G[i*n+j]!=0)\r
+        mem=(int)(log(G[i*n+j])/log(2.0));\r
+      if(mem>max_mem_x[i])\r
+        max_mem_x[i]=mem;\r
+      if(mem>max_mem)\r
+        max_mem=mem;\r
+    }\r
+  }\r
   \r
-  // generate the PS, PI tables for later use\r
-  for(int i=0;i<d_S;i++) {\r
-    int j=0;\r
-    for(int ii=0;ii<d_S;ii++) for(int jj=0;jj<d_I;jj++) {\r
-      if(d_NS[ii*d_I+jj]!=i) continue;\r
-      d_PS[i*d_I+j]=ii;\r
-      d_PI[i*d_I+j]=jj;\r
-      j++;\r
+//printf("max_mem_x\n");\r
+//for(int j=0;j<max_mem_x.size();j++) printf("%d ",max_mem_x[j]); printf("\n");\r
+\r
+  // calculate total memory requirements to set S\r
+  int sum_max_mem = 0;\r
+  for(int i=0;i<k;i++)\r
+    sum_max_mem += max_mem_x[i];\r
+\r
+//printf("sum_max_mem = %d\n",sum_max_mem);\r
+\r
+  d_I=1<<k;\r
+  d_S=1<<sum_max_mem;\r
+  d_O=1<<n;\r
\r
+  // binary representation of the G matrix\r
+  std::vector<std::vector<int> > Gb(k*n);\r
+  for(int j=0;j<k*n;j++) {\r
+    Gb[j].resize(max_mem+1);\r
+    dec2base(G[j],2,Gb[j]);\r
+//printf("Gb\n");\r
+//for(int m=0;m<Gb[j].size();m++) printf("%d ",Gb[j][m]); printf("\n");\r
+  }\r
+\r
+  // alphabet size of each shift register \r
+  std::vector<int> bases_x(k);\r
+  for(int j=0;j<k ;j++) \r
+    bases_x[j] = 1 << max_mem_x[j];\r
+//printf("bases_x\n");\r
+//for(int j=0;j<max_mem_x.size();j++) printf("%d ",max_mem_x[j]); printf("\n");\r
+\r
+  d_NS.resize(d_I*d_S);\r
+  d_OS.resize(d_I*d_S);\r
+\r
+  std::vector<int> sx(k);\r
+  std::vector<int> nsx(k);\r
+  std::vector<int> tx(k);\r
+  std::vector<std::vector<int> > tb(k);\r
+  for(int j=0;j<k;j++)\r
+    tb[j].resize(max_mem+1);\r
+  std::vector<int> inb(k);\r
+  std::vector<int> outb(n);\r
+\r
+\r
+  for(int s=0;s<d_S;s++) {\r
+    dec2bases(s,bases_x,sx); // split s into k values, each representing on of the k shift registers\r
+//printf("state = %d \nstates = ",s);\r
+//for(int j=0;j<sx.size();j++) printf("%d ",sx[j]); printf("\n");\r
+    for(int i=0;i<d_I;i++) {\r
+      dec2base(i,2,inb); // input in binary\r
+//printf("input = %d \ninputs = ",i);\r
+//for(int j=0;j<inb.size();j++) printf("%d ",inb[j]); printf("\n");\r
+\r
+      // evaluate next state\r
+      for(int j=0;j<k;j++)\r
+        nsx[j] = (inb[j]*bases_x[j]+sx[j])/2; // next state (for each shift register) MSB first\r
+      d_NS[s*d_I+i]=bases2dec(nsx,bases_x); // collect all values into the new state\r
+\r
+      // evaluate transitions\r
+      for(int j=0;j<k;j++)\r
+        tx[j] = inb[j]*bases_x[j]+sx[j]; // transition (for each shift register)MSB first\r
+      for(int j=0;j<k;j++) {\r
+        dec2base(tx[j],2,tb[j]); // transition in binary\r
+//printf("transition = %d \ntransitions = ",tx[j]);\r
+//for(int m=0;m<tb[j].size();m++) printf("%d ",tb[j][m]); printf("\n");\r
+      }\r
+\r
+      // evaluate outputs\r
+      for(int nn=0;nn<n;nn++) {\r
+        outb[nn] = 0;\r
+        for(int j=0;j<k;j++) {\r
+          for(int m=0;m<max_mem+1;m++)\r
+            outb[nn] = (outb[nn] + Gb[j*n+nn][m]*tb[j][m]) % 2; // careful: polynomial 1+D ir represented as 110, not as 011\r
+//printf("output %d equals %d\n",nn,outb[nn]);\r
+        }\r
+      }\r
+      d_OS[s*d_I+i] = base2dec(outb,2);\r
     }\r
   }\r
+\r
+  generate_PS_PI();\r
+  generate_TM();\r
 }\r
 \r
+\r
+\r
+\r
 //######################################################################\r
 //# Automatically generate an FSM specification describing the \r
 //# ISI for a channel\r
 //# of length ch_length and a modulation of size mod_size\r
 //######################################################################\r
-fsm::fsm(const int mod_size, const int ch_length)\r
+fsm::fsm(int mod_size, int ch_length)\r
 {\r
   d_I=mod_size;\r
   d_S=(int) (pow(1.0*d_I,1.0*ch_length-1)+0.5);\r
@@ -125,8 +222,6 @@ fsm::fsm(const int mod_size, const int ch_length)
 \r
   d_NS.resize(d_I*d_S);\r
   d_OS.resize(d_I*d_S);\r
-  d_PS.resize(d_I*d_S);\r
-  d_PI.resize(d_I*d_S);\r
 \r
   for(int s=0;s<d_S;s++) {\r
     for(int i=0;i<d_I;i++) { \r
@@ -135,8 +230,20 @@ fsm::fsm(const int mod_size, const int ch_length)
       d_OS[s*d_I+i] = t;\r
     }\r
   }\r
-  \r
-  // generate the PS, PI tables for later use\r
\r
+  generate_PS_PI();\r
+  generate_TM();\r
+}\r
+\r
+\r
+//######################################################################\r
+//# generate the PS and PI tables for later use\r
+//######################################################################\r
+void fsm::generate_PS_PI()\r
+{\r
+  d_PS.resize(d_I*d_S);\r
+  d_PI.resize(d_I*d_S);\r
+\r
   for(int i=0;i<d_S;i++) {\r
     int j=0;\r
     for(int ii=0;ii<d_S;ii++) for(int jj=0;jj<d_I;jj++) {\r
@@ -147,3 +254,70 @@ fsm::fsm(const int mod_size, const int ch_length)
     }\r
   }\r
 }\r
+\r
+\r
+//######################################################################\r
+//# generate the termination matrices TMl and TMi for later use\r
+//######################################################################\r
+void fsm::generate_TM()\r
+{\r
+  d_TMi.resize(d_S*d_S);\r
+  d_TMl.resize(d_S*d_S);\r
+\r
+  for(int i=0;i<d_S*d_S;i++) {\r
+    d_TMi[i] = -1; // no meaning\r
+    d_TMl[i] = d_S; //infinity: you need at most S-1 steps\r
+    if (i/d_S == i%d_S)\r
+      d_TMl[i] = 0;\r
+  }\r
+\r
+  for(int s=0;s<d_S;s++) {\r
+    bool done = false;\r
+    int attempts = 0;\r
+    while (done == false && attempts < d_S-1) {\r
+      done = find_es(s);\r
+      attempts ++;\r
+    }\r
+    if (done == false)\r
+      //throw std::runtime_error ("fsm::generate_TM(): FSM appears to be disconnected\n");\r
+      printf("fsm::generate_TM(): FSM appears to be disconnected\n");\r
+  }\r
+}\r
+\r
+\r
+// find a path from any state to the ending state "es"\r
+bool fsm::find_es(int es)\r
+{\r
+  bool done = true;\r
+  for(int s=0;s<d_S;s++) {\r
+    if(d_TMl[s*d_S+es] < d_S) \r
+      continue;\r
+    int minl=d_S;\r
+    int mini=-1;\r
+    for(int i=0;i<d_I;i++) {\r
+      if( 1 + d_TMl[d_NS[s*d_I+i]*d_S+es] < minl) {\r
+        minl = 1 + d_TMl[d_NS[s*d_I+i]*d_S+es];\r
+        mini = i;\r
+      }\r
+    }\r
+    if (mini != -1) {\r
+      d_TMl[s*d_S+es]=minl;\r
+      d_TMi[s*d_S+es]=mini;\r
+    }\r
+    else\r
+      done = false;\r
+  }\r
+  return done;\r
+}\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r
+\r