1 <?xml version="1.0" encoding="ISO-8859-1"?>
2 <!DOCTYPE article PUBLIC "-//OASIS//DTD DocBook XML V4.2//EN"
4 <!ENTITY test_tcm_listing SYSTEM "test_tcm.py">
10 <title>Trellis-based algorithms for GNU Radio</title>
12 <firstname>Achilleas</firstname>
13 <surname>Anastasopoulos</surname>
16 <email>anastas@umich.edu</email>
23 <revnumber>v0.0</revnumber>
24 <date>2006-08-03</date>
32 <abstract><para>This document provides a description of the
33 Finite State Machine (FSM) implementation and the related
34 trellis-based encoding and decoding algorithms
43 <!--=====================================================-->
44 <sect1 id="intro"><title>Introduction</title>
49 The basic goal of the implementation is to have a generic way of
50 describing an FSM that is decoupled from whether it describes a
52 code (CC), a trellis code (TC), an inter-symbol interference (ISI)
54 other communication system that can be modeled with an FSM.
55 To achieve this goal, we need to separate the pure FSM descrition from the
56 rest of the model details. For instance, in the case of a rate 2/3 TC,
57 the FSM should not involve details about the modulation used (it can
58 be an 8-ary PAM, or 8-PSK, etc). Similarly, when attempting maximum likelihood
59 sequence detection (MLSD)--using for instance the Viterbi algorithm (VA)--
60 the VA implementation should not be concerned with the channel details
61 (such as modulations, channel type, hard or soft inputs, etc).
62 Clearly, having generality as the primary goal implies some penalty
63 on the code efficiency, as compared to fully custom implementations.
67 We will now describe the implementation of the basic ingedient, the FSM.
73 <!--=====================================================-->
74 <sect1 id="fsm"><title>The FSM class</title>
76 <para>An FSM describes the evolution of a system with inputs
77 x<subscript>k</subscript>, states s<subscript>k</subscript> and outputs y<subscript>k</subscript>. At time k the FSM state is s<subscript>k</subscript>.
78 Upon reception of a new input symbol x<subscript>k</subscript>, it outputs an output symbol
79 y<subscript>k</subscript> which is a function of both x<subscript>k</subscript> and s<subscript>k</subscript>.
80 It will then move to a next state s<subscript>k+1</subscript>.
81 An FSM has a finite number of states, input and output symbols.
82 All these are formally described as follows:
86 <listitem><para>The input alphabet A<subscript>I</subscript>={0,1,2,...,I-1}, with cardinality I, so that x<subscript>k</subscript> takes values in A<subscript>I</subscript>.</para></listitem>
87 <listitem><para>The state alphabet A<subscript>S</subscript>={0,1,2,...,S-1}, with cardinality S, so that s<subscript>k</subscript> takes values in A<subscript>S</subscript>.</para></listitem>
88 <listitem><para>The output alphabet A<subscript>O</subscript>={0,1,2,...,O-1}, with cardinality O, so that y<subscript>k</subscript> takes values in A<subscript>O</subscript></para></listitem>
89 <listitem><para>The "next-state" function NS: A<subscript>S</subscript> x A<subscript>I</subscript> --> A<subscript>S</subscript>,
91 s<subscript>k+1</subscript> = NS(s<subscript>k</subscript>, x<subscript>k</subscript>)</para></listitem>
92 <listitem><para>The "output-symbol" function OS: A<subscript>S</subscript> x A<subscript>I</subscript> --> A<subscript>S</subscript>,
94 y<subscript>k</subscript> = OS(s<subscript>k</subscript>, x<subscript>k</subscript>)</para></listitem>
98 Thus, a complete description of the FSM is given by the
99 the five-tuple (I,S,O,NS,OS).
100 Observe that implementation details are hidden
101 in how the outside world interprets these input and output
103 Here is an example of an FSM describing the (2,1) CC
104 with constraint length 3 and generator polynomial matrix
105 (1+D+D<superscript>2</superscript> , 1+D<superscript>2</superscript>)
106 from Proakis-Salehi pg. 779.
110 <example id="cc_ex"><title>(2,1) CC with generator polynomials (1+D+D<superscript>2</superscript> , 1+D<superscript>2</superscript>)</title>
113 This CC accepts 1 bit at a time, and outputs 2 bits at a time.
114 It has a shift register storing the last two input bits.
116 b<subscript>k</subscript>(0)=x<subscript>k</subscript>+
117 x<subscript>k-1</subscript>+x<subscript>k-2</subscript>, and
118 b<subscript>k</subscript>(1)=x<subscript>k</subscript>+
119 x<subscript>k-2</subscript>, where addition is mod-2.
120 We can represent the state of this system
121 as s<subscript>k</subscript> = (x<subscript>k-1</subscript> x<subscript>k-2</subscript>)<subscript>10</subscript>. In addition we can represent its
122 output symbol as y<subscript>k</subscript> = (b<subscript>k</subscript>(1) b<subscript>k</subscript>(0))<subscript>10</subscript>.
123 Based on the above assumptions, the input alphabet A<subscript>I</subscript>={0,1}, so I=2;
124 the state alphabet A<subscript>S</subscript>={0,1,2,3}, so S=4; and
125 the output alphabet A<subscript>O</subscript>={0,1,2,3}, so O=4.
126 The "next-state" function NS(,) is given by
128 s<subscript>k</subscript> x<subscript>k</subscript> s<subscript>k+1</subscript>
138 The "output-symbol" function OS(,) can be given by
140 s<subscript>k</subscript> x<subscript>k</subscript> y<subscript>k</subscript>
153 Note that although the CC outputs 2 bits per time period, following
154 our approach, there is only one (quaternary) output symbol per
155 time period (for instance, here we use the decimal representation
156 of the 2-bits). Also note that the modulation used is not part of
157 the FSM description: it can be BPSK, OOK, BFSK, QPSK with or without Gray mapping, etc;
158 it is up to the rest of the program to interpret the meaning of
159 the symbol y<subscript>k</subscript>.
166 The C++ implementation of the FSM class keeps private information about
167 I,S,O,NS,OS and public methods to read and write them. The NS
168 and OS matrices are implemented as STL 1-dimensional vectors.
177 std::vector<int> d_NS;
178 std::vector<int> d_OS;
179 std::vector<int> d_PS;
180 std::vector<int> d_PI;
183 fsm(const fsm &FSM);
184 fsm(const int I, const int S, const int O, const std::vector<int> &NS, const std::vector<int> &OS);
185 fsm(const char *name);
186 fsm(const int mod_size, const int ch_length);
187 int I () const { return d_I; }
188 int S () const { return d_S; }
189 int O () const { return d_O; }
190 const std::vector<int> & NS () const { return d_NS; }
191 const std::vector<int> & OS () const { return d_OS; }
192 const std::vector<int> & PS () const { return d_PS; }
193 const std::vector<int> & PI () const { return d_PI; }
198 As can be seen, other than the trivial and the copy constructor,
199 there are three additional
200 ways to construct an FSM.
205 <para>Supplying the parameters I,S,O,NS,OS:</para>
207 fsm(const int I, const int S, const int O, const std::vector<int> &NS, const std::vector<int> &OS);
212 <para>Giving a filename containing all the FSM information:</para>
214 fsm(const char *name);
217 This information has to be in the following format:
221 NS(0,0) NS(0,1) ... NS(0,I-1)
222 NS(1,0) NS(1,1) ... NS(1,I-1)
224 NS(S-1,0) NS(S-1,1) ... NS(S-1,I-1)
226 OS(0,0) OS(0,1) ... OS(0,I-1)
227 OS(1,0) OS(1,1) ... OS(1,I-1)
229 OS(S-1,0) OS(S-1,1) ... OS(S-1,I-1)
233 For instance, the file containing the information for the example mentioned above is of the form:
251 <para>The third way is specific to FSMs resulting from shift registers, and the output symbol being the entire transition (ie, current_state and current_input). These FSMs are usefull when describibg ISI channels. In particular the state is comprised of the.....
254 fsm(const int mod_size, const int ch_length);
262 Finally, as can be seen from the above description, there are
263 two more variables included in the FSM class implementation,
264 the PS and the PI matrices. These are computed internally
265 when an FSM is instantiated and their meaning is as follows.
266 Sometimes (eg, in the traceback operation of the VA) we need
267 to trace the history of the state or the input sequence.
268 To do this we would like to know for a given state s<subscript>k</subscript>, what are the possible previous states s<subscript>k-1</subscript>
269 and what input symbols x<subscript>k-1</subscript> will get us from
270 s<subscript>k-1</subscript> to s<subscript>k</subscript>. This information can be derived from NS; however we want to have it ready in a
272 In the following we assume that for any state,
273 the number of incoming transitions is the same as the number of
274 outgoing transitions, ie, equal to I. All applications of interest
275 have FSMs satisfying this requirement.
277 If we arbitrarily index the incoming transitions to the current state
278 by "i", then as i goes from 0 to I-1, PS(s<subscript>k</subscript>,i)
279 gives all previous states s<subscript>k-1</subscript>,
280 and PI(s<subscript>k</subscript>,i) gives all previous inputs x<subscript>k-1</subscript>.
281 In other words, for any given s<subscript>k</subscript> and any index i=0,1,...I-1, starting from
282 s<subscript>k-1</subscript>=PS(s<subscript>k</subscript>,i)
284 x<subscript>k-1</subscript>=PI(s<subscript>k</subscript>,i)
285 will get us to the state s<subscript>k</subscript>.
286 More formally, for any i=0,1,...I-1 we have
287 s<subscript>k</subscript> = NS(PS(s<subscript>k</subscript>,i),PI(s<subscript>k</subscript>,i)).
298 <!--=====================================================-->
299 <sect1 id="tcm"><title>TCM: A Complete Example</title>
302 We now discuss through a concrete example how
303 the above FSM model can be used in GNU Radio.
305 The communication system that we want to simulate
306 consists of a source generating the
307 input information in packets, a CC encoding each packet separately,
308 a memoryless modulator,
309 an additive white Gaussian noise (AWGN) channel, and
310 the VA performing MLSD.
311 The program source is as follows.
317 #!/usr/bin/env python
319 from gnuradio import gr
320 from gnuradio import audio
321 from gnuradio import trellis
322 from gnuradio import eng_notation
328 def run_test (f,Kb,bitspersymbol,K,dimensionality,constellation,N0,seed):
329 fg = gr.flow_graph ()
332 src = gr.lfsr_32k_source_s()
333 src_head = gr.head (gr.sizeof_short,Kb/16) # packet size in shorts
334 s2fsmi = gr.packed_to_unpacked_ss(bitspersymbol,gr.GR_MSB_FIRST) # unpack shorts to symbols compatible with the FSM input cardinality
335 enc = trellis.encoder_ss(f,0) # initial state = 0
336 mod = gr.chunks_to_symbols_sf(constellation,dimensionality)
340 noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
343 metrics = trellis.metrics_f(f.O(),dimensionality,constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
344 va = trellis.viterbi_s(f,K,0,-1) # Put -1 if the Initial/Final states are not set.
345 fsmi2s = gr.unpacked_to_packed_ss(bitspersymbol,gr.GR_MSB_FIRST) # pack FSM input symbols to shorts
346 dst = gr.check_lfsr_32k_s();
348 fg.connect (src,src_head,s2fsmi,enc,mod)
349 fg.connect (mod,(add,0))
350 fg.connect (noise,(add,1))
351 fg.connect (add,metrics)
352 fg.connect (metrics,va,fsmi2s,dst)
356 # A bit of cheating: run the program once and print the
357 # final encoder state..
358 # Then put it as the last argument in the viterbi block
359 #print "final state = " , enc.ST()
361 ntotal = dst.ntotal ()
362 nright = dst.nright ()
363 runlength = dst.runlength ()
364 return (ntotal,ntotal-nright)
371 esn0_db=float(args[1]) # Es/No in dB
372 rep=int(args[2]) # number of times the experiment is run to collect enough errors
374 sys.stderr.write ('usage: test_tcm.py fsm_fname Es/No_db repetitions\n')
378 f=trellis.fsm(fname) # get the FSM specification from a file (will hopefully be automated in the future...)
379 Kb=1024*16 # packet size in bits (make it multiple of 16 so it can be packed in a short)
380 bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
381 K=Kb/bitspersymbol # packet size in trellis steps
382 modulation = fsm_utils.psk4 # see fsm_utlis.py for available predefined modulations
383 dimensionality = modulation[0]
384 constellation = modulation[1]
385 if len(constellation)/dimensionality != f.O():
386 sys.stderr.write ('Incompatible FSM output cardinality and modulation size.\n')
388 # calculate average symbol energy
390 for i in range(len(constellation)):
391 Es = Es + constellation[i]**2
392 Es = Es / (len(constellation)/dimensionality)
393 N0=Es/pow(10.0,esn0_db/10.0); # noise variance
398 (s,e)=run_test(f,Kb,bitspersymbol,K,dimensionality,constellation,N0,-long(666+i)) # run experiment with different seed to get different noise realizations
402 print i,s,e,tot_s,terr_s, '%e' % ((1.0*terr_s)/tot_s)
403 # estimate of the (short) error rate
404 print tot_s,terr_s, '%e' % ((1.0*terr_s)/tot_s)
407 if __name__ == '__main__':
412 The program is called by
415 ./test_tcm.py fsm_fname Es/No_db repetitions
418 where "fsm_fname" is the file containing the FSM specification of the
419 tested TCM code, "Es/No_db" is the SNR in dB, and "repetitions"
420 are the number of packets to be transmitted and received in order to
421 collect sufficient number of errors for an accurate estimate of the
426 The FSM is first intantiated in "main" by
439 Each packet has size Kb bits
440 (we choose Kb to be a multiple of 16 so that all bits fit nicely into shorts and can be generated by the lfsr GNU Radio).
441 Assuming that the FSM input has cardinality I, each input symbol consists
442 of bitspersymbol=log<subscript>2</subscript>( I ). The Kb/16 shorts are now
443 unpacked to K=Kb/bitspersymbol input
444 symbols that will drive the FSM encoder.
447 Kb=1024*16 # packet size in bits (make it multiple of 16 so it can be packed in a short)
448 bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
449 K=Kb/bitspersymbol # packet size in trellis steps
455 The FSM will produce K output symbols (remeber the FSM produces always one output symbol for each input symbol). Each of these symbols needs to be modulated. Since we are simulating the communication system, we need not simulate the actual waveforms. An M-ary, N-dimensional
456 modulation is completely specified by a set of M, N-dimensional real vectors. In "fsm_utils.py" file we give a number of useful modulations with the following format: modulation = (N,constellation), where
457 constellation=[c11,c12,...,c1N,c21,c22,...,c2N,...,cM1,cM2,...cMN].
458 The meaning of the above is that every constellation point c_i
459 is an N-dimnsional vector c_i=(ci1,ci2,...,ciN)
460 For instance, 4-ary PAM is represented as
461 (1,[-3, -1, 1, 3]), while QPSK is represented as
462 (2,[1, 0, 0, 1, 0, -1, -1, 0]). In our example we choose QPSK modulation.
463 Clearly, M should be equal to the cardinality of the FSM output, O.
464 Finally the average symbol energy and noise variance are calculated.
467 modulation = fsm_utils.psk4 # see fsm_utlis.py for available predefined modulations
468 dimensionality = modulation[0]
469 constellation = modulation[1]
470 if len(constellation)/dimensionality != f.O():
471 sys.stderr.write ('Incompatible FSM output cardinality and modulation size.\n')
473 # calculate average symbol energy
475 for i in range(len(constellation)):
476 Es = Es + constellation[i]**2
477 Es = Es / (len(constellation)/dimensionality)
478 N0=Es/pow(10.0,esn0_db/10.0); # noise variance
484 Then, "run_test" is called with a different "seed" so that we get
485 different noise realizations.
488 (s,e)=run_test(f,Kb,bitspersymbol,K,dimensionality,constellation,N0,-long(666+i)) # run experiment with different seed to get different noise realizations
494 Let us examine now the "run_test" function.
495 First we set up the transmitter blocks.
496 The Kb/16 shorts are first unpacked to
497 symbols consistent with the FSM input alphabet.
498 The FSm encoder requires the FSM specification,
499 and an initial state (which is set to 0 in this example).
503 src = gr.lfsr_32k_source_s()
504 src_head = gr.head (gr.sizeof_short,Kb/16) # packet size in shorts
505 s2fsmi = gr.packed_to_unpacked_ss(bitspersymbol,gr.GR_MSB_FIRST) # unpack shorts to symbols compatible with the FSM input cardinality
506 enc = trellis.encoder_ss(f,0) # initial state = 0
514 The "chunks_to_symbols_sf" is essentially a memoryless mapper which
515 for each input symbol y_k
516 outputs a sequence of N numbers ci1,ci2,...,ciN representing the
517 coordianates of the constellation symbol c_i with i=y_k.
520 mod = gr.chunks_to_symbols_sf(constellation,dimensionality)
524 The channel is AWGN with appropriate noise variance.
525 For each transmitted symbol c_k=(ck1,ck2,...,ckN) we receive a noisy version
526 r_k=(rk1,rk2,...,rkN).
531 noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
537 Part of the design methodology was to decouple the FSM and VA from
538 the details of the modulation, channel, receiver front-end etc.
539 In order for the VA to run, we only need to provide it with
540 a number representing a cost associated with each transition
541 in the trellis. Then the VA will find the sequence with
542 the smallest total cost through the trellis.
543 The cost associated with a transition (s_k,x_k) is only a function
544 of the output y_k = OS(s_k,x_k), and the observation
545 vector r_k. Thus, for each time period, k,
546 we need to label each of the SxI transitions with such a cost.
547 This means that for each time period we need to evaluate
548 O such numbers (one for each possible output symbol y_k).
550 in "metrics_f". In particular, metrics_f is a memoryless device
551 taking N inputs at a time and producing O outputs. The N inputs are
554 are the costs associated with observations rk1,rk2,...,rkN and
555 hypothesized output symbols c_1,c_2,...,c_M. For instance,
556 if we choose to perform soft-input VA, we need to evaluate
557 the Euclidean distance between r_k and each of c_1,c_2,...,c_M,
558 for each of the K transmitted symbols.
559 Other options are available as well; for instance, we can
560 do hard decision demodulation and feed the VA with
561 symbol Hamming distances, or even bit Hamming distances, etc.
562 These are all implemented in "metrics_f".
566 metrics = trellis.metrics_f(f.O(),dimensionality,constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
570 Now the VA can run once it is supplied by the initial and final states.
571 The initial state is known to be 0; the final state is usually
572 forced to some value by padding the information sequence appropriately.
573 In this example, we always send the the same info sequence (we only randomize noise) so we can evaluate off line the final state and then provide it to the VA (a value of -1 signifies that there is no fixed initial
574 or final state). The VA outputs the estimates of the symbols x_k which
575 are then packed to shorts and compared with the transmitted sequence.
578 va = trellis.viterbi_s(f,K,0,-1) # Put -1 if the Initial/Final states are not set.
579 fsmi2s = gr.unpacked_to_packed_ss(bitspersymbol,gr.GR_MSB_FIRST) # pack FSM input symbols to shorts
580 dst = gr.check_lfsr_32k_s();
587 The function returns the number of shorts and the number of shorts in error. Observe that this way the estimated error rate refers to
588 16-bit-symbol error rate.
591 return (ntotal,ntotal-nright)
599 <!--=====================================================-->
600 <sect1 id="future"><title>Future Work</title>
608 Improve the documentation :-)
614 automate fsm generation from generator polynomials
615 (feedforward or feedback form).
621 Optimize the VA code.
627 Provide implementation of soft-input soft-output (SISO) decoders for
628 potential use in concatenated systems. Also a host of suboptimal
629 decoders, eg, sphere decoding, M- and T- algorithms, sequential decoding, etc.
637 Although turbo decoding is rediculously slow in software,
638 we can design it in principle. One question is, whether we should
639 use the encoder, and SISO blocks and connect them
640 through GNU radio or should we implement turbo-decoding
641 as a single block (issues with buffering between blocks).