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);
216 <para>This information has to be in the following format</para>
220 NS(0,0) NS(0,1) ... NS(0,I-1)
221 NS(1,0) NS(1,1) ... NS(1,I-1)
223 NS(S-1,0) NS(S-1,1) ... NS(S-1,I-1)
225 OS(0,0) OS(0,1) ... OS(0,I-1)
226 OS(1,0) OS(1,1) ... OS(1,I-1)
228 OS(S-1,0) OS(S-1,1) ... OS(S-1,I-1)
230 <para>For instance, the file containing the information for the example mentioned above is of the form</para>
247 <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.....
250 fsm(const int mod_size, const int ch_length);
258 Finally, as can be seen from the above description, there are
259 two more variables included in the FSM class implementation,
260 the PS and the PI matrices. These are computed internally
261 when an FSM is instantiated and their meaning is as follows.
262 Sometimes (eg, in the traceback operation of the VA) we need
263 to trace the history of the state or the input sequence.
264 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>
265 and what input symbols x<subscript>k-1</subscript> will get us from
266 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
268 In the following we assume that for any state,
269 the number of incoming transitions is the same as the number of
270 outgoing transitions, ie, equal to I. All applications of interest
271 have FSMs satisfying this requirement.
273 If we arbitrarily index the incoming transitions to the current state
274 by "i", then as i goes from 0 to I-1, PS(s<subscript>k</subscript>,i)
275 gives all previous states s<subscript>k-1</subscript>,
276 and PI(s<subscript>k</subscript>,i) gives all previous inputs x<subscript>k-1</subscript>.
277 In other words, for any given s<subscript>k</subscript> and any index i=0,1,...I-1, starting from
278 s<subscript>k-1</subscript>=PS(s<subscript>k</subscript>,i)
280 x<subscript>k-1</subscript>=PI(s<subscript>k</subscript>,i)
281 will get us to the state s<subscript>k</subscript>.
282 More formally, for any i=0,1,...I-1 we have
283 s<subscript>k</subscript> = NS(PS(s<subscript>k</subscript>,i),PI(s<subscript>k</subscript>,i)).
294 <!--=====================================================-->
295 <sect1 id="tcm"><title>TCM: A Complete Example</title>
298 We now discuss through a concrete example how
299 the above FSM model can be used in GNU Radio.
301 The communication system that we want to simulate
302 consists of a source generating the
303 input information in packets, a CC encoding each packet separately,
304 a memoryless modulator,
305 an additive white Gaussian noise (AWGN) channel, and
306 the VA performing MLSD.
307 The program source is as follows.
313 #!/usr/bin/env python
315 from gnuradio import gr
316 from gnuradio import audio
317 from gnuradio import trellis
318 from gnuradio import eng_notation
324 def run_test (f,Kb,bitspersymbol,K,dimensionality,constellation,N0,seed):
325 fg = gr.flow_graph ()
329 src = gr.lfsr_32k_source_s()
330 src_head = gr.head (gr.sizeof_short,Kb/16) # packet size in shorts
331 s2fsmi = gr.packed_to_unpacked_ss(bitspersymbol,gr.GR_MSB_FIRST) # unpack shorts to symbols compatible with the FSM input cardinality
332 enc = trellis.encoder_ss(f,0) # initial state = 0
333 mod = gr.chunks_to_symbols_sf(constellation,dimensionality)
338 noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
342 metrics = trellis.metrics_f(f.O(),dimensionality,constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
343 va = trellis.viterbi_s(f,K,0,-1) # Put -1 if the Initial/Final states are not set.
344 fsmi2s = gr.unpacked_to_packed_ss(bitspersymbol,gr.GR_MSB_FIRST) # pack FSM input symbols to shorts
345 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)
357 # A bit of cheating: run the program once and print the
358 # final encoder state..
359 # Then put it as the last argument in the viterbi block
360 #print "final state = " , enc.ST()
362 ntotal = dst.ntotal ()
363 nright = dst.nright ()
364 runlength = dst.runlength ()
365 return (ntotal,ntotal-nright)
374 esn0_db=float(args[1]) # Es/No in dB
375 rep=int(args[2]) # number of times the experiment is run to collect enough errors
377 sys.stderr.write ('usage: test_tcm.py fsm_fname Es/No_db repetitions\n')
381 f=trellis.fsm(fname) # get the FSM specification from a file (will hopefully be automated in the future...)
382 Kb=1024*16 # packet size in bits (make it multiple of 16 so it can be packed in a short)
383 bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
384 K=Kb/bitspersymbol # packet size in trellis steps
385 modulation = fsm_utils.psk4 # see fsm_utlis.py for available predefined modulations
386 dimensionality = modulation[0]
387 constellation = modulation[1]
388 if len(constellation)/dimensionality != f.O():
389 sys.stderr.write ('Incompatible FSM output cardinality and modulation size.\n')
391 # calculate average symbol energy
393 for i in range(len(constellation)):
394 Es = Es + constellation[i]**2
395 Es = Es / (len(constellation)/dimensionality)
396 N0=Es/pow(10.0,esn0_db/10.0); # noise variance
403 (s,e)=run_test(f,Kb,bitspersymbol,K,dimensionality,constellation,N0,-long(666+i)) # run experiment with different seed to get different noise realizations
407 print i,s,e,tot_s,terr_s, '%e' % ((1.0*terr_s)/tot_s)
408 # estimate of the (short) error rate
409 print tot_s,terr_s, '%e' % ((1.0*terr_s)/tot_s)
412 if __name__ == '__main__':
417 The program is called by
420 ./test_tcm.py fsm_fname Es/No_db repetitions
423 where "fsm_fname" is the file containing the FSM specification of the
424 tested TCM code, "Es/No_db" is the SNR in dB, and "repetitions"
425 are the number of packets to be transmitted and received in order to
426 collect sufficient number of errors for an accurate estimate of the
431 The FSM is first intantiated in "main" by
444 Each packet has size Kb bits
445 (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).
446 Assuming that the FSM input has cardinality I, each input symbol consists
447 of bitspersymbol=log<subscript>2</subscript>( I ). The Kb/16 shorts are now
448 unpacked to K=Kb/bitspersymbol input
449 symbols that will drive the FSM encoder.
452 Kb=1024*16 # packet size in bits (make it multiple of 16 so it can be packed in a short)
453 bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
454 K=Kb/bitspersymbol # packet size in trellis steps
460 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
461 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
462 constellation=[c11,c12,...,c1N,c21,c22,...,c2N,...,cM1,cM2,...cMN].
463 The meaning of the above is that every constellation point c_i
464 is an N-dimnsional vector c_i=(ci1,ci2,...,ciN)
465 For instance, 4-ary PAM is represented as
466 (1,[-3, -1, 1, 3]), while QPSK is represented as
467 (2,[1, 0, 0, 1, 0, -1, -1, 0]). In our example we choose QPSK modulation.
468 Clearly, M should be equal to the cardinality of the FSM output, O.
469 Finally the average symbol energy and noise variance are calculated.
472 modulation = fsm_utils.psk4 # see fsm_utlis.py for available predefined modulations
473 dimensionality = modulation[0]
474 constellation = modulation[1]
475 if len(constellation)/dimensionality != f.O():
476 sys.stderr.write ('Incompatible FSM output cardinality and modulation size.\n')
478 # calculate average symbol energy
480 for i in range(len(constellation)):
481 Es = Es + constellation[i]**2
482 Es = Es / (len(constellation)/dimensionality)
483 N0=Es/pow(10.0,esn0_db/10.0); # noise variance
489 Then, "run_test" is called with a different "seed" so that we get
490 different noise realizations.
493 (s,e)=run_test(f,Kb,bitspersymbol,K,dimensionality,constellation,N0,-long(666+i)) # run experiment with different seed to get different noise realizations
499 Let us examine now the "run_test" function.
500 First we set up the transmitter blocks.
501 The Kb/16 shorts are first unpacked to
502 symbols consistent with the FSM input alphabet.
503 The FSm encoder requires the FSM specification,
504 and an initial state (which is set to 0 in this example).
508 src = gr.lfsr_32k_source_s()
509 src_head = gr.head (gr.sizeof_short,Kb/16) # packet size in shorts
510 s2fsmi = gr.packed_to_unpacked_ss(bitspersymbol,gr.GR_MSB_FIRST) # unpack shorts to symbols compatible with the FSM input cardinality
511 enc = trellis.encoder_ss(f,0) # initial state = 0
519 The "chunks_to_symbols_sf" is essentially a memoryless mapper which
520 for each input symbol y_k
521 outputs a sequence of N numbers ci1,ci2,...,ciN representing the
522 coordianates of the constellation symbol c_i with i=y_k.
525 mod = gr.chunks_to_symbols_sf(constellation,dimensionality)
529 The channel is AWGN with appropriate noise variance.
530 For each transmitted symbol c_k=(ck1,ck2,...,ckN) we receive a noisy version
531 r_k=(rk1,rk2,...,rkN).
535 noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
541 Part of the design methodology was to decouple the FSM and VA from
542 the details of the modulation, channel, receiver front-end etc.
543 In order for the VA to run, we only need to provide it with
544 a number representing a cost associated with each transition
545 in the trellis. Then the VA will find the sequence with
546 the smallest total cost through the trellis.
547 The cost associated with a transition (s_k,x_k) is only a function
548 of the output y_k = OS(s_k,x_k), and the observation
549 vector r_k. Thus, for each time period, k,
550 we need to label each of the SxI transitions with such a cost.
551 This means that for each time period we need to evaluate
552 O such numbers (one for each possible output symbol y_k).
554 in "metrics_f". In particular, metrics_f is a memoryless device
555 taking N inputs at a time and producing O outputs. The N inputs are
558 are the costs associated with observations rk1,rk2,...,rkN and
559 hypothesized output symbols c_1,c_2,...,c_M. For instance,
560 if we choose to perform soft-input VA, we need to evaluate
561 the Euclidean distance between r_k and each of c_1,c_2,...,c_M,
562 for each of the K transmitted symbols.
563 Other options are available as well; for instance, we can
564 do hard decision demodulation and feed the VA with
565 symbol Hamming distances, or even bit Hamming distances, etc.
566 These are all implemented in "metrics_f".
570 metrics = trellis.metrics_f(f.O(),dimensionality,constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
574 Now the VA can run once it is supplied by the initial and final states.
575 The initial state is known to be 0; the final state is usually
576 forced to some value by padding the information sequence appropriately.
577 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
578 or final state). The VA outputs the estimates of the symbols x_k which
579 are then packed to shorts and compared with the transmitted sequence.
582 va = trellis.viterbi_s(f,K,0,-1) # Put -1 if the Initial/Final states are not set.
583 fsmi2s = gr.unpacked_to_packed_ss(bitspersymbol,gr.GR_MSB_FIRST) # pack FSM input symbols to shorts
584 dst = gr.check_lfsr_32k_s();
591 The function returns the number of shorts and the number of shorts in error. Observe that this way the estimated error rate refers to
592 16-bit-symbol error rate.
595 return (ntotal,ntotal-nright)
603 <!--=====================================================-->
604 <sect1 id="future"><title>Future Work</title>
612 Improve the documentation :-)
618 automate fsm generation from generator polynomials
619 (feedforward or feedback form).
625 Optimize the VA code.
631 Provide implementation of soft-input soft-output (SISO) decoders for
632 potential use in concatenated systems. Also a host of suboptimal
633 decoders, eg, sphere decoding, M- and T- algorithms, sequential decoding, etc.
640 Although turbo decoding is rediculously slow in software,
641 we can design it in pronciple. The question is, should
642 we use the FSM and SISO abstractions and cnnect them
643 through GNU radio or should we implement turbo-decoding
644 as a single block (issues with buffering between blocks).