git.gag.com Git - debian/gnuradio/blob - gr-trellis/src/lib/fsm.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002 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 #include <cstdio>
  24 #include <string>
  25 #include <iostream>
  26 #include <fstream>
  27 #include <stdexcept>
  28 #include <cmath>
  29 #include <stdlib.h>
  30 #include "base.h"
  31 #include "fsm.h"
  32
  33
  34 fsm::fsm()
  35 {
  36   d_I=0;
  37   d_S=0;
  38   d_O=0;
  39   d_NS.resize(0);
  40   d_OS.resize(0);
  41   d_PS.resize(0);
  42   d_PI.resize(0);
  43   d_TMi.resize(0);
  44   d_TMl.resize(0);
  45 }
  46
  47 fsm::fsm(const fsm &FSM)
  48 {
  49   d_I=FSM.I();
  50   d_S=FSM.S();
  51   d_O=FSM.O();
  52   d_NS=FSM.NS();
  53   d_OS=FSM.OS();
  54   d_PS=FSM.PS(); // is this going to make a deep copy?
  55   d_PI=FSM.PI();
  56   d_TMi=FSM.TMi();
  57   d_TMl=FSM.TMl();
  58 }
  59
  60 fsm::fsm(int I, int S, int O, const std::vector<int> &NS, const std::vector<int> &OS)
  61 {
  62   d_I=I;
  63   d_S=S;
  64   d_O=O;
  65   d_NS=NS;
  66   d_OS=OS;
  67
  68   generate_PS_PI();
  69   generate_TM();
  70 }
  71
  72 //######################################################################
  73 //# Read an FSM specification from a file.
  74 //# Format (hopefully will become more flexible in the future...):
  75 //# I S O (in the first line)
  76 //# blank line
  77 //# Next state matrix (S lines, each with I integers separated by spaces)
  78 //# blank line
  79 //# output symbol matrix (S lines, each with I integers separated by spaces)
  80 //# optional comments
  81 //######################################################################
  82 fsm::fsm(const char *name)
  83 {
  84   FILE *fsmfile;
  85
  86   if((fsmfile=fopen(name,"r"))==NULL)
  87     throw std::runtime_error ("fsm::fsm(const char *name): file open error\n");
  88     //printf("file open error in fsm()\n");
  89
  90   if(fscanf(fsmfile,"%d %d %d\n",&d_I,&d_S,&d_O) == EOF) {
  91     if(ferror(fsmfile) != 0)
  92       throw std::runtime_error ("fsm::fsm(const char *name): file read error\n");
  93   }
  94
  95   d_NS.resize(d_I*d_S);
  96   d_OS.resize(d_I*d_S);
  97
  98   for(int i=0;i<d_S;i++) {
  99     for(int j=0;j<d_I;j++) {
 100       if(fscanf(fsmfile,"%d",&(d_NS[i*d_I+j])) == EOF) {
 101         if(ferror(fsmfile) != 0)
 102           throw std::runtime_error ("fsm::fsm(const char *name): file read error\n");
 103       }
 104     }
 105   }
 106   for(int i=0;i<d_S;i++) {
 107     for(int j=0;j<d_I;j++) {
 108       if(fscanf(fsmfile,"%d",&(d_OS[i*d_I+j])) == EOF) {
 109         if(ferror(fsmfile) != 0)
 110           throw std::runtime_error ("fsm::fsm(const char *name): file read error\n");
 111       }
 112     }
 113   }
 114
 115   generate_PS_PI();
 116   generate_TM();
 117 }
 118
 119
 120
 121 //######################################################################
 122 //# Automatically generate the FSM from the generator matrix
 123 //# of a (n,k) binary convolutional code
 124 //######################################################################
 125 fsm::fsm(int k, int n, const std::vector<int> &G)
 126 {
 127
 128   // calculate maximum memory requirements for each input stream
 129   std::vector<int> max_mem_x(k,-1);
 130   int max_mem = -1;
 131   for(int i=0;i<k;i++) {
 132     for(int j=0;j<n;j++) {
 133       int mem = -1;
 134       if(G[i*n+j]!=0)
 135         mem=(int)(log(G[i*n+j])/log(2.0));
 136       if(mem>max_mem_x[i])
 137         max_mem_x[i]=mem;
 138       if(mem>max_mem)
 139         max_mem=mem;
 140     }
 141   }
 142
 143 //printf("max_mem_x\n");
 144 //for(int j=0;j<max_mem_x.size();j++) printf("%d ",max_mem_x[j]); printf("\n");
 145
 146   // calculate total memory requirements to set S
 147   int sum_max_mem = 0;
 148   for(int i=0;i<k;i++)
 149     sum_max_mem += max_mem_x[i];
 150
 151 //printf("sum_max_mem = %d\n",sum_max_mem);
 152
 153   d_I=1<<k;
 154   d_S=1<<sum_max_mem;
 155   d_O=1<<n;
 156
 157   // binary representation of the G matrix
 158   std::vector<std::vector<int> > Gb(k*n);
 159   for(int j=0;j<k*n;j++) {
 160     Gb[j].resize(max_mem+1);
 161     dec2base(G[j],2,Gb[j]);
 162 //printf("Gb\n");
 163 //for(int m=0;m<Gb[j].size();m++) printf("%d ",Gb[j][m]); printf("\n");
 164   }
 165
 166   // alphabet size of each shift register
 167   std::vector<int> bases_x(k);
 168   for(int j=0;j<k ;j++)
 169     bases_x[j] = 1 << max_mem_x[j];
 170 //printf("bases_x\n");
 171 //for(int j=0;j<max_mem_x.size();j++) printf("%d ",max_mem_x[j]); printf("\n");
 172
 173   d_NS.resize(d_I*d_S);
 174   d_OS.resize(d_I*d_S);
 175
 176   std::vector<int> sx(k);
 177   std::vector<int> nsx(k);
 178   std::vector<int> tx(k);
 179   std::vector<std::vector<int> > tb(k);
 180   for(int j=0;j<k;j++)
 181     tb[j].resize(max_mem+1);
 182   std::vector<int> inb(k);
 183   std::vector<int> outb(n);
 184
 185
 186   for(int s=0;s<d_S;s++) {
 187     dec2bases(s,bases_x,sx); // split s into k values, each representing one of the k shift registers
 188 //printf("state = %d \nstates = ",s);
 189 //for(int j=0;j<sx.size();j++) printf("%d ",sx[j]); printf("\n");
 190     for(int i=0;i<d_I;i++) {
 191       dec2base(i,2,inb); // input in binary
 192 //printf("input = %d \ninputs = ",i);
 193 //for(int j=0;j<inb.size();j++) printf("%d ",inb[j]); printf("\n");
 194
 195       // evaluate next state
 196       for(int j=0;j<k;j++)
 197         nsx[j] = (inb[j]*bases_x[j]+sx[j])/2; // next state (for each shift register) MSB first
 198       d_NS[s*d_I+i]=bases2dec(nsx,bases_x); // collect all values into the new state
 199
 200       // evaluate transitions
 201       for(int j=0;j<k;j++)
 202         tx[j] = inb[j]*bases_x[j]+sx[j]; // transition (for each shift register)MSB first
 203       for(int j=0;j<k;j++) {
 204         dec2base(tx[j],2,tb[j]); // transition in binary
 205 //printf("transition = %d \ntransitions = ",tx[j]);
 206 //for(int m=0;m<tb[j].size();m++) printf("%d ",tb[j][m]); printf("\n");
 207       }
 208
 209       // evaluate outputs
 210       for(int nn=0;nn<n;nn++) {
 211         outb[nn] = 0;
 212         for(int j=0;j<k;j++) {
 213           for(int m=0;m<max_mem+1;m++)
 214             outb[nn] = (outb[nn] + Gb[j*n+nn][m]*tb[j][m]) % 2; // careful: polynomial 1+D ir represented as 110, not as 011
 215 //printf("output %d equals %d\n",nn,outb[nn]);
 216         }
 217       }
 218       d_OS[s*d_I+i] = base2dec(outb,2);
 219     }
 220   }
 221
 222   generate_PS_PI();
 223   generate_TM();
 224 }
 225
 226
 227
 228
 229 //######################################################################
 230 //# Automatically generate an FSM specification describing the
 231 //# ISI for a channel
 232 //# of length ch_length and a modulation of size mod_size
 233 //######################################################################
 234 fsm::fsm(int mod_size, int ch_length)
 235 {
 236   d_I=mod_size;
 237   d_S=(int) (pow(1.0*d_I,1.0*ch_length-1)+0.5);
 238   d_O=d_S*d_I;
 239
 240   d_NS.resize(d_I*d_S);
 241   d_OS.resize(d_I*d_S);
 242
 243   for(int s=0;s<d_S;s++) {
 244     for(int i=0;i<d_I;i++) {
 245       int t=i*d_S+s;
 246       d_NS[s*d_I+i] = t/d_I;
 247       d_OS[s*d_I+i] = t;
 248     }
 249   }
 250
 251   generate_PS_PI();
 252   generate_TM();
 253 }
 254
 255
 256
 257
 258 //######################################################################
 259 //# Automatically generate an FSM specification describing the
 260 //# the trellis for a CPM with h=K/P (relatively prime),
 261 //# alphabet size M, and frequency pulse duration L symbols
 262 //#
 263 //# This FSM is based on the paper by B. Rimoldi
 264 //# "A decomposition approach to CPM", IEEE Trans. Info Theory, March 1988
 265 //# See also my own notes at http://www.eecs.umich.edu/~anastas/docs/cpm.pdf
 266 //######################################################################
 267 fsm::fsm(int P, int M, int L)
 268 {
 269   d_I=M;
 270   d_S=(int)(pow(1.0*M,1.0*L-1)+0.5)*P;
 271   d_O=(int)(pow(1.0*M,1.0*L)+0.5)*P;
 272
 273   d_NS.resize(d_I*d_S);
 274   d_OS.resize(d_I*d_S);
 275   int nv;
 276   for(int s=0;s<d_S;s++) {
 277     for(int i=0;i<d_I;i++) {
 278       int s1=s/P;
 279       int v=s%P;
 280       int ns1= (i*(int)(pow(1.0*M,1.0*(L-1))+0.5)+s1)/M;
 281       if (L==1)
 282         nv=(i+v)%P;
 283       else
 284         nv=(s1%M+v)%P;
 285       d_NS[s*d_I+i] = ns1*P+nv;
 286       d_OS[s*d_I+i] = i*d_S+s;
 287     }
 288   }
 289
 290   generate_PS_PI();
 291   generate_TM();
 292 }
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303 //######################################################################
 304 //# Automatically generate an FSM specification describing the
 305 //# the joint trellis of fsm1 and fsm2
 306 //######################################################################
 307 fsm::fsm(const fsm &FSM1, const fsm &FSM2)
 308 {
 309   d_I=FSM1.I()*FSM2.I();
 310   d_S=FSM1.S()*FSM2.S();
 311   d_O=FSM1.O()*FSM2.O();
 312
 313   d_NS.resize(d_I*d_S);
 314   d_OS.resize(d_I*d_S);
 315
 316   for(int s=0;s<d_S;s++) {
 317     for(int i=0;i<d_I;i++) {
 318       int s1=s/FSM2.S();
 319       int s2=s%FSM2.S();
 320       int i1=i/FSM2.I();
 321       int i2=i%FSM2.I();
 322       d_NS[s*d_I+i] = FSM1.NS()[s1 * FSM1.I() + i1] * FSM2.S() + FSM2.NS()[s2 * FSM2.I() + i2];
 323       d_OS[s*d_I+i] = FSM1.OS()[s1 * FSM1.I() + i1] * FSM2.O() + FSM2.OS()[s2 * FSM2.I() + i2];
 324     }
 325   }
 326
 327   generate_PS_PI();
 328   generate_TM();
 329 }
 330
 331
 332
 333
 334 //######################################################################
 335 //# Generate a new FSM representing n stages through the original FSM
 336 //# AKA radix-n FSM
 337 //######################################################################
 338 fsm::fsm(const fsm &FSM, int n)
 339 {
 340   d_I=(int) (pow(1.0*FSM.I(),1.0*n)+0.5);
 341   d_S=FSM.S();
 342   d_O=(int) (pow(1.0*FSM.O(),1.0*n)+0.5);
 343
 344   d_NS.resize(d_I*d_S);
 345   d_OS.resize(d_I*d_S);
 346
 347   for(int s=0;s<d_S;s++ ) {
 348     for(int i=0;i<d_I;i++ ) {
 349       std::vector<int> ii(n);
 350       dec2base(i,FSM.I(),ii);
 351       std::vector<int> oo(n);
 352       int ns=s;
 353       for(int k=0;k<n;k++) {
 354         oo[k]=FSM.OS()[ns*FSM.I()+ii[k]];
 355         ns=FSM.NS()[ns*FSM.I()+ii[k]];
 356       }
 357       d_NS[s*d_I+i]=ns;
 358       d_OS[s*d_I+i]=base2dec(oo,FSM.O());
 359     }
 360   }
 361
 362   generate_PS_PI();
 363   generate_TM();
 364 }
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374 //######################################################################
 375 //# generate the PS and PI tables for later use
 376 //######################################################################
 377 void fsm::generate_PS_PI()
 378 {
 379   d_PS.resize(d_S);
 380   d_PI.resize(d_S);
 381
 382   for(int i=0;i<d_S;i++) {
 383     d_PS[i].resize(d_I*d_S); // max possible size
 384     d_PI[i].resize(d_I*d_S);
 385     int j=0;
 386     for(int ii=0;ii<d_S;ii++) for(int jj=0;jj<d_I;jj++) {
 387       if(d_NS[ii*d_I+jj]!=i) continue;
 388       d_PS[i][j]=ii;
 389       d_PI[i][j]=jj;
 390       j++;
 391     }
 392     d_PS[i].resize(j);
 393     d_PI[i].resize(j);
 394   }
 395 }
 396
 397
 398 //######################################################################
 399 //# generate the termination matrices TMl and TMi for later use
 400 //######################################################################
 401 void fsm::generate_TM()
 402 {
 403   d_TMi.resize(d_S*d_S);
 404   d_TMl.resize(d_S*d_S);
 405
 406   for(int i=0;i<d_S*d_S;i++) {
 407     d_TMi[i] = -1; // no meaning
 408     d_TMl[i] = d_S; //infinity: you need at most S-1 steps
 409     if (i/d_S == i%d_S)
 410       d_TMl[i] = 0;
 411   }
 412
 413   for(int s=0;s<d_S;s++) {
 414     bool done = false;
 415     int attempts = 0;
 416     while (done == false && attempts < d_S-1) {
 417       done = find_es(s);
 418       attempts ++;
 419     }
 420     if (done == false) {
 421       //throw std::runtime_error ("fsm::generate_TM(): FSM appears to be disconnected\n");
 422       printf("fsm::generate_TM(): FSM appears to be disconnected\n");
 423       printf("state %d cannot be reached from all other states\n",s);
 424     }
 425   }
 426 }
 427
 428
 429 // find a path from any state to the ending state "es"
 430 bool fsm::find_es(int es)
 431 {
 432   bool done = true;
 433   for(int s=0;s<d_S;s++) {
 434     if(d_TMl[s*d_S+es] < d_S)
 435       continue;
 436     int minl=d_S;
 437     int mini=-1;
 438     for(int i=0;i<d_I;i++) {
 439       if( 1 + d_TMl[d_NS[s*d_I+i]*d_S+es] < minl) {
 440         minl = 1 + d_TMl[d_NS[s*d_I+i]*d_S+es];
 441         mini = i;
 442       }
 443     }
 444     if (mini != -1) {
 445       d_TMl[s*d_S+es]=minl;
 446       d_TMi[s*d_S+es]=mini;
 447     }
 448     else
 449       done = false;
 450   }
 451   return done;
 452 }
 453
 454
 455
 456
 457
 458 //######################################################################
 459 //#  generate trellis representation of FSM as an SVG file
 460 //######################################################################
 461 void fsm::write_trellis_svg( std::string filename ,int number_stages)
 462 {
 463    std::ofstream trellis_fname (filename.c_str());
 464    if (!trellis_fname) {std::cout << "file not found " << std::endl ; exit(-1);}
 465    const int TRELLIS_Y_OFFSET = 30;
 466    const int TRELLIS_X_OFFSET = 20;
 467    const int STAGE_LABEL_Y_OFFSET = 25;
 468    const int STAGE_LABEL_X_OFFSET = 20;
 469    const int STATE_LABEL_Y_OFFSET = 30;
 470    const int STATE_LABEL_X_OFFSET = 5;
 471    const int STAGE_STATE_OFFSETS = 10;
 472 //   std::cout << "################## BEGIN SVG TRELLIS PIC #####################" << std::endl;
 473    trellis_fname << "<svg viewBox = \"0 0 200 200\" version = \"1.1\">" << std::endl;
 474
 475     for( int stage_num = 0;stage_num < number_stages;stage_num ++){
 476     // draw states
 477       for ( int state_num = 0;state_num < d_S ; state_num ++ ) {
 478         trellis_fname << "<circle cx = \"" << stage_num * STAGE_STATE_OFFSETS + TRELLIS_X_OFFSET <<
 479         "\" cy = \"" << state_num * STAGE_STATE_OFFSETS + TRELLIS_Y_OFFSET << "\" r = \"1\"/>" << std::endl;
 480       //draw branches
 481         if(stage_num != number_stages-1){
 482           for( int branch_num = 0;branch_num < d_I; branch_num++){
 483             trellis_fname << "<line x1 =\"" << STAGE_STATE_OFFSETS * stage_num+ TRELLIS_X_OFFSET  << "\" ";
 484             trellis_fname << "y1 =\"" << state_num * STAGE_STATE_OFFSETS + TRELLIS_Y_OFFSET<< "\" ";
 485             trellis_fname << "x2 =\"" <<  STAGE_STATE_OFFSETS *stage_num + STAGE_STATE_OFFSETS+ TRELLIS_X_OFFSET << "\" ";
 486             trellis_fname << "y2 =\"" << d_NS[d_I * state_num + branch_num] * STAGE_STATE_OFFSETS + TRELLIS_Y_OFFSET << "\" ";
 487             trellis_fname << " stroke-dasharray = \"3," <<  branch_num << "\" ";
 488             trellis_fname << " stroke = \"black\" stroke-width = \"0.3\"/>" << std::endl;
 489           }
 490         }
 491       }
 492     }
 493   // label the stages
 494   trellis_fname << "<g font-size = \"4\" font= \"times\" fill = \"black\">" << std::endl;
 495   for( int stage_num = 0;stage_num < number_stages ;stage_num ++){
 496     trellis_fname << "<text x = \"" << stage_num * STAGE_STATE_OFFSETS + STAGE_LABEL_X_OFFSET <<
 497       "\" y = \""  << STAGE_LABEL_Y_OFFSET  << "\" >" << std::endl;
 498     trellis_fname << stage_num <<  std::endl;
 499     trellis_fname << "</text>" << std::endl;
 500   }
 501   trellis_fname << "</g>" << std::endl;
 502
 503   // label the states
 504   trellis_fname << "<g font-size = \"4\" font= \"times\" fill = \"black\">" << std::endl;
 505   for( int state_num = 0;state_num < d_S ; state_num ++){
 506     trellis_fname << "<text y = \"" << state_num * STAGE_STATE_OFFSETS + STATE_LABEL_Y_OFFSET <<
 507       "\" x = \""  << STATE_LABEL_X_OFFSET  << "\" >" << std::endl;
 508     trellis_fname << state_num <<  std::endl;
 509     trellis_fname << "</text>" << std::endl;
 510   }
 511   trellis_fname << "</g>" << std::endl;
 512
 513
 514   trellis_fname << "</svg>" << std::endl;
 515 //  std::cout << "################## END SVG TRELLIS PIC ##################### " << std::endl;
 516   trellis_fname.close();
 517 }
 518
 519
 520
 521
 522
 523
 524 //######################################################################
 525 //# Write trellis specification to a text file,
 526 //# in the same format used when reading FSM files
 527 //######################################################################
 528 void fsm::write_fsm_txt(std::string filename)
 529 {
 530    std::ofstream trellis_fname (filename.c_str());
 531    if (!trellis_fname) {std::cout << "file not found " << std::endl ; exit(-1);}
 532    trellis_fname << d_I << ' ' << d_S << ' ' << d_O << std::endl;
 533    trellis_fname << std::endl;
 534    for(int i=0;i<d_S;i++) {
 535      for(int j=0;j<d_I;j++)  trellis_fname << d_NS[i*d_I+j] << ' ';
 536      trellis_fname << std::endl;
 537    }
 538    trellis_fname << std::endl;
 539    for(int i=0;i<d_S;i++) {
 540      for(int j=0;j<d_I;j++) trellis_fname << d_OS[i*d_I+j] << ' ';
 541      trellis_fname << std::endl;
 542    }
 543    trellis_fname << std::endl;
 544    trellis_fname.close();
 545 }
 546