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