Imported Upstream version 3.0
[debian/gnuradio] / gr-trellis / doc / gr-trellis.xml
1 <?xml version="1.0" encoding="ISO-8859-1"?>
2 <!DOCTYPE article PUBLIC "-//OASIS//DTD DocBook XML V4.2//EN"
3           "docbookx.dtd" [
4   <!ENTITY test_tcm_listing SYSTEM "test_tcm.py.xml">
5   <!ENTITY test_viterbi_equalization1_listing SYSTEM "test_viterbi_equalization1.py.xml">
6 ]>
7
8 <article>
9
10 <articleinfo>
11   <title>Trellis-based algorithms for GNU Radio</title>
12   <author>
13      <firstname>Achilleas</firstname>
14      <surname>Anastasopoulos</surname>
15      <affiliation>
16         <address>
17            <email>anastas@umich.edu</email>
18         </address>
19      </affiliation>
20   </author>
21
22 <!--
23 <revhistory>
24   <revision>
25   <revnumber>v0.0</revnumber>
26   <date>2006-08-03</date>
27   <revremark>
28     First cut.
29   </revremark>
30   </revision>
31 </revhistory>
32 -->
33
34 <abstract><para>This document provides a description of the 
35 Finite State Machine (FSM) implementation and the related 
36 trellis-based encoding and decoding algorithms
37 for GNU Radio.
38 </para></abstract>
39
40 </articleinfo>
41
42
43
44
45 <!--=====================================================-->
46 <sect1 id="intro"><title>Introduction</title>
47
48 <para>
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 
51 convolutional 
52 code (CC), a trellis code (TC), an inter-symbol interference (ISI) 
53 channel, or any
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. 
64 </para>
65
66 <para>
67 We will now describe the implementation of the basic ingedient, the FSM.
68 </para>
69
70 </sect1>
71
72
73 <!--=====================================================-->
74 <sect1 id="fsm"><title>The FSM class</title>
75
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:
83 </para>
84
85 <itemizedlist>
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>, 
90 with the meaning 
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>, 
93 with the meaning 
94 y<subscript>k</subscript> = OS(s<subscript>k</subscript>, x<subscript>k</subscript>)</para></listitem>
95 </itemizedlist>
96
97 <para>
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 
102 integer symbols.
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.
107 </para>
108
109
110 <example id="cc_ex"><title>(2,1) CC with generator polynomials (1+D+D<superscript>2</superscript> , 1+D<superscript>2</superscript>)</title>
111
112 <para>
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.
115 In particular, 
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
127 <programlisting>
128 s<subscript>k</subscript>       x<subscript>k</subscript>       s<subscript>k+1</subscript>
129 0       0       0
130 0       1       2
131 1       0       0
132 1       1       2
133 2       0       1
134 2       1       3
135 3       0       1
136 3       1       3
137 </programlisting>
138 The "output-symbol" function OS(,) can be given by
139 <programlisting>
140 s<subscript>k</subscript>       x<subscript>k</subscript>       y<subscript>k</subscript>
141 0       0       0
142 0       1       3
143 1       0       3
144 1       1       0
145 2       0       1
146 2       1       2
147 3       0       2
148 3       1       1
149 </programlisting>
150 </para>
151
152 <para>
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>.
160 </para>
161
162 </example>
163
164
165 <para>
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.
169 </para>
170
171 <programlisting>
172 class fsm {
173 private:
174   int d_I;
175   int d_S;
176   int d_O;
177   std::vector&lt;int&gt; d_NS;
178   std::vector&lt;int&gt; d_OS;
179   std::vector&lt;int&gt; d_PS;
180   std::vector&lt;int&gt; d_PI;
181   std::vector&lt;int&gt; d_TMi;
182   std::vector&lt;int&gt; d_TMl;
183   void generate_PS_PI ();
184   void generate_TM ();
185   bool find_es(int es);
186 public:
187   fsm();
188   fsm(const fsm &amp;FSM);
189   fsm(int I, int S, int O, const std::vector&lt;int&gt; &amp;NS, const std::vector&lt;int&gt; &amp;OS);
190   fsm(const char *name);
191   fsm(int k, int n, const std::vector&lt;int&gt; &amp;G);
192   fsm(int mod_size, int ch_length);
193   int I () const { return d_I; }
194   int S () const { return d_S; }
195   int O () const { return d_O; }
196   const std::vector&lt;int&gt; &amp; NS () const { return d_NS; }
197   const std::vector&lt;int&gt; &amp; OS () const { return d_OS; }
198   const std::vector&lt;int&gt; &amp; PS () const { return d_PS; }
199   const std::vector&lt;int&gt; &amp; PI () const { return d_PI; }
200   const std::vector&lt;int&gt; &amp; TMi () const { return d_TMi; }
201   const std::vector&lt;int&gt; &amp; TMl () const { return d_TMl; }
202 };
203 </programlisting>
204
205 <para>
206 As can be seen, other than the trivial and the copy constructor, 
207 there are three additional
208 ways to construct an FSM. 
209 </para>
210
211 <itemizedlist>
212 <listitem>
213 <para>Supplying the parameters I,S,O,NS,OS:</para>
214 <programlisting>
215   fsm(const int I, const int S, const int O, const std::vector&lt;int&gt; &amp;NS, const std::vector&lt;int&gt; &amp;OS);
216 </programlisting>
217 </listitem>
218
219 <listitem>
220 <para>Giving a filename containing all the FSM information:</para>
221 <programlisting>
222   fsm(const char *name);
223 </programlisting>
224 <para>
225 This information has to be in the following format:
226 <programlisting>
227 I S O
228
229 NS(0,0)   NS(0,1)   ...  NS(0,I-1)
230 NS(1,0)   NS(1,1)   ...  NS(1,I-1)
231 ...
232 NS(S-1,0) NS(S-1,1) ...  NS(S-1,I-1)
233
234 OS(0,0)   OS(0,1)   ...  OS(0,I-1)
235 OS(1,0)   OS(1,1)   ...  OS(1,I-1)
236 ...
237 OS(S-1,0) OS(S-1,1) ... OS(S-1,I-1)
238 </programlisting>
239 </para>
240 <para>
241 For instance, the file containing the information for the example mentioned above is of the form:
242 <programlisting>
243 2 4 4
244
245 0 2
246 0 2
247 1 3
248 1 3
249
250 0 3
251 3 0
252 1 2
253 2 1
254 </programlisting>
255 </para>
256 </listitem>
257
258
259 <listitem>
260 <para>
261 The third way is specific to FSMs representing binary (n,k) conolutional codes. These FSMs are specified by the number of input bits k, the number of output bits n, and the generator matrix, which is a k x n matrix of integers 
262 G = [g<subscript>i,j</subscript>]<subscript>i=1:k, j=1:n</subscript>, given as an one-dimensional STL vector.
263 Each integer g<subscript>i,j</subscript> is the decimal representation of the 
264 polynomial g<subscript>i,j</subscript>(D) (e.g., g<subscript>i,j</subscript>= 6 = 110<subscript>2</subscript> is interpreted as g<subscript>i,j</subscript>(D)=1+D) describing the connections from  the sequence x<subscript>i</subscript> to 
265 y<subscript>j</subscript> (e.g., in the above example y<subscript>j</subscript>(k) = x<subscript>i</subscript>(k) + x<subscript>i</subscript>(k-1)).
266 </para>
267 <programlisting>
268   fsm(int k, int n, const std::vector&lt;int&gt; &amp;G);
269 </programlisting>
270 </listitem>
271
272
273 <listitem>
274 <para>
275 The fourth 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 input symbols x(k-1), x(k-2),...,x(k-L), where L = ch_length-1 and each x(i) belongs to an alphabet of size mod_size. The output is taken to be x(k), x(k-1), x(k-2),...,x(k-L) (in decimal format)
276 </para>
277 <programlisting>
278   fsm(const int mod_size, const int ch_length);
279 </programlisting>
280 </listitem>
281
282 </itemizedlist>
283
284
285 <para>
286 As can be seen from the above description, there are
287 two more variables included in the FSM class implementation, 
288 the PS and the PI matrices. These are computed internally 
289 when an FSM is instantiated and their meaning is as follows.
290 Sometimes (eg, in the traceback operation of the VA) we need
291 to trace the history of the state or the input sequence. 
292 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>
293 and what input symbols x<subscript>k-1</subscript> will get us from 
294 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 
295 convenient format. 
296 In the following we assume that for any state, 
297 the number of incoming transitions is the same as the number of 
298 outgoing transitions, ie, equal to I. All applications of interest 
299 have FSMs satisfying this requirement.
300
301 If we arbitrarily index the incoming transitions to the current state 
302 by "i", then  as i goes from 0 to I-1, PS(s<subscript>k</subscript>,i)
303 gives all previous states s<subscript>k-1</subscript>,
304 and PI(s<subscript>k</subscript>,i) gives all previous inputs x<subscript>k-1</subscript>.
305 In other words, for any given s<subscript>k</subscript> and any index i=0,1,...I-1, starting from 
306 s<subscript>k-1</subscript>=PS(s<subscript>k</subscript>,i)
307 with input 
308 x<subscript>k-1</subscript>=PI(s<subscript>k</subscript>,i)
309 will get us to the state s<subscript>k</subscript>. 
310 More formally, for any i=0,1,...I-1 we have
311 s<subscript>k</subscript> = NS(PS(s<subscript>k</subscript>,i),PI(s<subscript>k</subscript>,i)).
312
313 </para>
314
315
316 <para>
317 Finally, there are
318 two more variables included in the FSM class implementation,
319 the TMl and the TMi matrices. These are both S x S matrices (represented as STL vectors) computed internally
320 when an FSM is instantiated and their meaning is as follows.
321 TMl(i,j) is the minimum number of trellis steps required to go from state i to state j.
322 Similarly, TMi(i,j) is the initial input required to get you from state i to state j in the minimum number of steps. As an example, if TMl(1,4)=2, it means that you need 2 steps in the trellis to get from state 1 to state 4. Further,
323 if TMi(1,4)=0 it means that the first such step will be followed if when at state 1 the input is 0. Furthermore, suppose that NS(1,0)=2. Then, TMl(2,4) should be 1 (ie, one more step is needed when starting from state 2 and having state 4 as the final destination). Finally, TMi(2,4) will give us the second input required to complete the path from 1 to 4.
324 These matrices are useful when we want to implement an encoder with proper state termination. For instance, based on these matrices we can evaluate how many 
325 additional input symbols (and which particular inputs) are required to be appended at the end of an input sequence so that the final state is always 0.
326
327 </para>
328
329
330 </sect1>
331
332
333
334 <!--=====================================================-->
335 <sect1 id="blocks"><title>Blocks Using the FSM structure</title>
336
337 <para>
338 In this section we give a brief description of the basic blocks implemented that make use of the previously described FSM structure.
339 </para>
340
341
342 <!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -->
343 <sect2 id="encoder"><title>Trellis Encoder</title>
344 <para>
345 The <methodname>trellis.encoder_XX(FSM, ST)</methodname> block instantiates an FSM encoder corresponding to the fsm FSM and having initial state ST. The input and output is a sequence of bytes, shorts or integers. 
346 </para>
347 </sect2>
348
349 <!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -->
350 <sect2 id="decoder"><title>Viterbi Decoder</title>
351 <para>
352 The <methodname>trellis.viterbi_X(FSM, K, S0, SK)</methodname> block instantiates a Viterbi decoder 
353 for a sequence of K trellis steps generated by the given FSM and with initial and final states set to S0 and SK, respectively (S0 and/or SK are set to -1
354 if the corresponding states are not fixed/known at the receiver side).
355 The output of this block is a sequence of K bytes, shorts or integers representing the estimated input (i.e., uncoded) sequence.
356 The input is a sequence of K x FSM.O( ) floats, where the k x K + i 
357 float represents the cost associated with the k-th
358 step in the trellis and the i-th FSM output.
359 Observe that these inputs are generated externally and thus the Viterbi block is not informed of their meaning (they can be genarated as soft or hard inputs, etc); the only requirement is that they represent additive costs.
360 </para>
361 </sect2>
362
363 <!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -->
364 <sect2 id="metrics"><title>Metrics Calculator</title>
365 <para>
366 The <methodname>trellis.metrics_X(O,D,TABLE,TYPE)</methodname> block is responsible for 
367 transforming the channel output to the stream of metrics appropriate as 
368 inputs to the Viterbi block described above. For each D input bytes/shorts/integers/floats/complexes it produces O output floats 
369
370 </para>
371
372 <para>
373
374 The parameter TYPE dictates how these metrics are generated:
375
376 <itemizedlist>
377 <listitem><para>
378 TRELLIS_EUCLIDEAN: for each D-dimensional vector 
379 r<subscript>k</subscript>=
380 (r<subscript>k,1</subscript>,r<subscript>k,2</subscript>,...,r<subscript>k,D</subscript>) 
381 evaluates
382 </para>
383 <para>
384 ||r<subscript>k</subscript>-c<subscript>i</subscript>||<superscript>2</superscript> = sum<subscript>j=1</subscript><superscript>D</superscript> |r<subscript>k,j</subscript>-c<subscript>i,j</subscript>|<superscript>2</superscript>
385 </para>
386 <para>
387 for each of the O hypothesized ouput
388 symbols c<subscript>i</subscript> = (c<subscript>i,1</subscript>,c<subscript>i,2</subscript>,...,c<subscript>i,D</subscript>) defined in the vector TABLE,
389 where TABLE[i * D + j] = c<subscript>i,j</subscript>.
390 </para></listitem>
391
392
393 <listitem><para>
394 TRELLIS_HARD_SYMBOL: for each D-dimensional vector 
395 r<subscript>k</subscript>=
396 (r<subscript>k,1</subscript>,r<subscript>k,2</subscript>,...,r<subscript>k,D</subscript>) 
397 evaluates
398 </para>
399 <para>
400 i<subscript>0</subscript>= argmax<subscript>i</subscript> ||r<subscript>k</subscript>-c<subscript>i</subscript>||<superscript>2</superscript> = 
401 argmax<subscript>i</subscript> sum<subscript>j=1</subscript><superscript>D</superscript> |r<subscript>k,j</subscript>-c<subscript>i,j</subscript>|<superscript>2</superscript>
402 </para>
403 <para>
404 and outputs a sequence of O floats of the form (0,...,0,1,0,...,0), where the 
405 i<subscript>0</subscript> position is set to 1. This corresponds to generating hard inputs (based on the symbol-wise Hamming distance) to the Viterbi algorithm.
406 </para></listitem>
407
408
409 <listitem><para>
410 TRELLIS_HARD_BIT (not yet implemented): for each D-dimensional vector 
411 r<subscript>k</subscript>=
412 (r<subscript>k,1</subscript>,r<subscript>k,2</subscript>,...,r<subscript>k,D</subscript>) 
413 evaluates
414 </para>
415 <para>
416 i<subscript>0</subscript>= argmax<subscript>i</subscript> ||r<subscript>k</subscript>-c<subscript>i</subscript>||<superscript>2</superscript> = 
417 argmax<subscript>i</subscript> sum<subscript>j=1</subscript><superscript>D</superscript> (r<subscript>k,j</subscript>-c<subscript>i,j</subscript>)<superscript>2</superscript>
418 </para>
419 <para>
420 and outputs a sequence of O floats of the form (d<subscript>1</subscript>,d<subscript>2</subscript>,...,d<subscript>O</subscript>), where the 
421 d<subscript>i</subscript> is the bitwise Hamming distance between i and  i<subscript>0</subscript>. This corresponds to generating hard inputs (based on the bit-wise Hamming distance) to the Viterbi algorithm.
422 </para></listitem>
423
424
425 </itemizedlist>
426
427
428 </para>
429 </sect2>
430
431
432
433
434 <!-- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -->
435 <sect2 id="viterbi_combined"><title>Combined Metrics Calculator and Viterbi Decoder</title>
436 <para>
437 Although the separation of metric calculation and Viterbi algorithm blocks
438 is consistent with our goal of providing general blocks that can be easily 
439 reused, this separation might result in large input/output buffer sizes
440 betwen blocks. Indeed for an FSM with a large output alphabet, the 
441 output of the metric block/input of the Viterbi block is FSM.O( ) floats for
442 each trellis step. Sometimes this results in buffer overflow even for
443 moderate sequence lengths.
444 To overcome this problem we provide a block that incorporates the metric calculation and Viterbi algorithm into a single GNU Radio block, namely
445 <methodname>trellis.viterbi_combined_X( FSM, K, S0, SK, D, TABLE, TYPE)</methodname> where the arguments are exactly those used in the aforementioned two blocks.
446 </para>
447 </sect2>
448
449
450
451
452
453 </sect1>
454
455
456 <!--=====================================================-->
457 <sect1 id="tcm"><title>A Complete Example: Trellis Coded Modulation (TCM)</title>
458
459 <para>
460 We now discuss through a concrete example how
461 the above FSM model can be used in GNU Radio.
462
463 The communication system that we want to simulate
464 consists of a source generating the
465 input information in packets, a CC encoding each packet separately, 
466 a memoryless modulator,
467 an additive white Gaussian noise (AWGN) channel, and
468 the VA performing MLSD.
469 The program source is as follows.
470 </para>
471
472 &test_tcm_listing;
473
474 <para>
475 The program is called by
476 <literallayout>
477 ./test_tcm.py fsm_fname Es/No_db repetitions
478 </literallayout>
479 where "fsm_fname" is the file containing the FSM specification of the
480 tested TCM code, "Es/No_db" is the SNR in dB, and "repetitions" 
481 are the number of packets to be transmitted and received in order to
482 collect sufficient number of errors for an accurate estimate of the
483 error rate.
484 </para>
485
486 <para>
487 The FSM is first intantiated in "main" by 
488 </para>
489 <programlisting>
490  62      f=trellis.fsm(fname) # get the FSM specification from a file
491 </programlisting>
492
493
494
495
496
497
498
499 <para>
500 Each packet has size Kb bits
501 (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).
502 Assuming that the FSM input has cardinality I, each input symbol consists 
503 of bitspersymbol=log<subscript>2</subscript>( I ). The Kb/16 shorts are now 
504 unpacked to K=Kb/bitspersymbol input
505 symbols that will drive the FSM encoder.
506 </para>
507 <programlisting>
508  63      Kb=1024*16  # packet size in bits (make it multiple of 16 so it can be packed in a short)
509  64      bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
510  65      K=Kb/bitspersymbol # packet size in trellis steps
511 </programlisting>
512
513
514
515 <para>
516 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, D-dimensional
517 modulation is completely specified by a set of M, D-dimensional real vectors. In "fsm_utils.py" file we give a number of useful modulations with the following format: modulation = (D,constellation), where
518 constellation=[c11,c12,...,c1D,c21,c22,...,c2D,...,cM1,cM2,...cMD].
519 The meaning of the above is that every constellation point c_i
520 is an D-dimnsional vector c_i=(ci1,ci2,...,ciD)
521 For instance, 4-ary PAM is represented as
522 (1,[-3, -1, 1, 3]), while QPSK is represented as
523 (2,[1, 0, 0, 1, 0, -1, -1, 0]). In our example we choose QPSK modulation.
524 Clearly, M should be equal to the cardinality of the FSM output, O.
525 Finally the average symbol energy and noise variance are calculated.
526 </para>
527 <programlisting>
528  66      modulation = fsm_utils.psk4 # see fsm_utlis.py for available predefined modulations
529  67      dimensionality = modulation[0]
530  68      constellation = modulation[1]
531  69      if len(constellation)/dimensionality != f.O():
532  70          sys.stderr.write (&apos;Incompatible FSM output cardinality and modulation size.\n&apos;)
533  71          sys.exit (1)
534  72      # calculate average symbol energy
535  73      Es = 0
536  74      for i in range(len(constellation)):
537  75          Es = Es + constellation[i]**2
538  76      Es = Es / (len(constellation)/dimensionality)
539  77      N0=Es/pow(10.0,esn0_db/10.0); # noise variance
540 </programlisting>
541
542
543
544 <para>
545 Then, "run_test" is called with a different "seed" so that we get
546 different noise realizations.
547 </para>
548 <programlisting>
549  82          (s,e)=run_test(f,Kb,bitspersymbol,K,dimensionality,constellation,N0,-long(666+i)) # run experiment with different seed to get different noise realizations
550 </programlisting>
551
552
553
554 <para>
555 Let us examine now the "run_test" function. 
556 First we set up the transmitter blocks.
557 The Kb/16 shorts are first unpacked to 
558 symbols consistent with the FSM input alphabet.
559 The FSm encoder requires the FSM specification,
560 and an initial state (which is set to 0 in this example).
561 </para>
562 <programlisting>
563  15      # TX
564  16      src = gr.lfsr_32k_source_s()
565  17      src_head = gr.head (gr.sizeof_short,Kb/16) # packet size in shorts
566  18      s2fsmi = gr.packed_to_unpacked_ss(bitspersymbol,gr.GR_MSB_FIRST) # unpack shorts to symbols compatible with the FSM input cardinality
567  19      enc = trellis.encoder_ss(f,0) # initial state = 0
568 </programlisting>
569
570
571
572
573 <para>
574 We now need to modulate the FSM output symbols.
575 The "chunks_to_symbols_sf" is essentially a memoryless mapper which 
576 for each input symbol y_k 
577 outputs a sequence of D numbers ci1,ci2,...,ciD  representing the 
578 coordianates of the constellation symbol c_i with i=y_k.
579 </para>
580 <programlisting>
581  20      mod = gr.chunks_to_symbols_sf(constellation,dimensionality)
582 </programlisting>
583
584 <para>
585 The channel is AWGN with appropriate noise variance.
586 For each transmitted symbol c_k=(ck1,ck2,...,ckD) we receive a noisy version
587 r_k=(rk1,rk2,...,rkD).
588 </para>
589 <programlisting>
590  22      # CHANNEL
591  23      add = gr.add_ff()
592  24      noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
593 </programlisting>
594
595
596
597 <para>
598 Part of the design methodology was to decouple the FSM and VA from
599 the details of the modulation, channel, receiver front-end etc.
600 In order for the VA to run, we only need to provide it with
601 a number representing a cost associated with each transition 
602 in the trellis. Then the VA will find the sequence with 
603 the smallest total cost through the trellis. 
604 The cost associated with a transition (s_k,x_k) is only a function
605 of the output y_k = OS(s_k,x_k), and the observation
606 vector r_k. Thus, for each time period, k,
607 we need to label each of the SxI transitions with such a cost.
608 This means that for each time period we need to evaluate 
609 O such numbers (one for each possible output symbol y_k). 
610 This is done 
611 in "metrics_f". In particular, metrics_f is a memoryless device
612 taking D inputs at a time and producing O outputs. The D inputs are
613 rk1,rk2,...,rkD.
614 The O outputs
615 are the costs associated with observations rk1,rk2,...,rkD and
616 hypothesized output symbols c_1,c_2,...,c_M. For instance,
617 if we choose to perform soft-input VA, we need to evaluate
618 the Euclidean distance between r_k and each of c_1,c_2,...,c_M,
619 for each of the K transmitted symbols.
620 Other options are available as well; for instance, we can
621 do hard decision demodulation and feed the VA with 
622 symbol Hamming distances, or even bit Hamming distances, etc.
623 These are all implemented in "metrics_f".
624 </para>
625 <programlisting>
626  26      # RX
627  27      metrics = trellis.metrics_f(f.O(),dimensionality,constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
628 </programlisting>
629
630 <para>
631 Now the VA can run once it is supplied by the initial and final states.
632 The initial state is known to be 0; the final state is usually 
633 forced to some value by padding the information sequence appropriately.
634 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
635 or final state). The VA outputs the estimates of the symbols x_k which
636 are then packed to shorts and compared with the transmitted sequence.
637 </para>
638 <programlisting>
639  28      va = trellis.viterbi_s(f,K,0,-1) # Put -1 if the Initial/Final states are not set.
640  29      fsmi2s = gr.unpacked_to_packed_ss(bitspersymbol,gr.GR_MSB_FIRST) # pack FSM input symbols to shorts
641  30      dst = gr.check_lfsr_32k_s();
642 </programlisting>
643
644
645
646
647 <para>
648 The function returns the number of shorts and the number of shorts in error. Observe that this way the estimated error rate refers to 
649 16-bit-symbol error rate.
650 </para>
651 <programlisting>
652  48      return (ntotal,ntotal-nright)
653 </programlisting>
654
655 </sect1>
656
657
658
659
660
661
662
663
664
665
666
667
668
669 <!--=====================================================-->
670 <sect1 id="isi"><title>Another Complete Example: Viterbi Equalization</title>
671
672 <para>
673 We now discuss through another concrete example how
674 the above FSM model can be used in GNU Radio.
675
676 The communication system that we want to simulate
677 consists of a source generating the
678 input information in packets, an ISI channel with
679 additive white Gaussian noise (AWGN), and
680 the VA performing MLSD.
681 The program source is as follows.
682 </para>
683
684 &test_viterbi_equalization1_listing;
685
686 <para>
687 The program is called by
688 <literallayout>
689 ./test_viterbi_equalization1.py Es/No_db repetitions
690 </literallayout>
691 where
692 "Es/No_db" is the SNR in dB, and "repetitions" 
693 are the number of packets to be transmitted and received in order to
694 collect sufficient number of errors for an accurate estimate of the
695 error rate.
696 </para>
697
698
699 <para>
700 Each packet has size Kb bits.
701 The modulation is chosen to be 4-PAM in this example and the channel is chosen 
702 to be one of the test channels  defined in fsm_utils.py
703 </para>
704 <programlisting>
705  71      Kb=2048  # packet size in bits
706  72      modulation = fsm_utils.pam4 # see fsm_utlis.py for available predefined modulations
707  73      channel = fsm_utils.c_channel # see fsm_utlis.py for available predefined test channels
708 </programlisting>
709
710 <para>
711 The FSM is instantiated in
712 </para>
713 <programlisting>
714  74      f=trellis.fsm(len(modulation[1]),len(channel)) # generate the FSM automatically
715 </programlisting>
716 <para>
717 and generated automatically given the channel length and the modulation size.
718 Since in this example the channel has length 5 and the modulation is 4-ary, the corresponding FSM has 4<superscript>5-1</superscript>=256 states and
719 4<superscript>5</superscript>=1024 outputs (see the documentation on FSM for more explanation).
720 </para>
721
722 <para>
723 Assuming that the FSM input has cardinality I, each input symbol consists 
724 of bitspersymbol=log<subscript>2</subscript>( I ) bits, and thus correspond to K = Kb/bitspersymbol symbols. 
725 </para>
726 <programlisting>
727  75      bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per FSM input symbol
728  76      K=Kb/bitspersymbol # packet size in trellis steps
729 </programlisting>
730
731
732
733 <para>
734 The overall system with input the 4-ary input symbols 
735 x<subscript>k</subscript>, modulated to the
736 4-PAM symbols s<subscript>k</subscript> and passed through the ISI channel to produce the
737 noise-free observations 
738 z<subscript>k</subscript> = 
739 sum<subscript>j=0</subscript><superscript>L-1</superscript> c<subscript>j</subscript> s<subscript>k-j</subscript> (where L is the channel length)
740 can be modeled as a FSM followed by a memoryless modulation. 
741 In particular, the FSM input is the sequence 
742 x<subscript>k</subscript>
743 and its output is the "combined" symbol 
744 y<subscript>k</subscript>=
745 (x<subscript>k</subscript>,x<subscript>k-1</subscript>,...,x<subscript>k-L+1</subscript>) (actually its decimal representation). 
746 The memoryless modulator maps every "combined" symbol
747 y<subscript>k</subscript> to 
748 z<subscript>k</subscript> = 
749 sum<subscript>j=0</subscript><superscript>L-1</superscript> c<subscript>j</subscript> s<subscript>k-j</subscript> 
750 Clearly this modulation is memoryless since 
751 each z<subscript>k</subscript> depends only on y<subscript>k</subscript>; the memory inherent in the ISI is hidden in the FSM structure.
752 This memoryless modulator is automatically generated by a helper function in 
753 fsm_utils.py given the channel and modulation as in line 78, and has the 
754 familiar format tot_channel=(dimensionality,tot_constellation) as described in the TCM example.
755 This is exactly what the metrics block (or the viterbi_combined block) require in order to evaluate the VA metrics from the noisy observations.
756 </para>
757 <programlisting>
758  78      tot_channel = fsm_utils.make_isi_lookup(modulation,channel,True) # generate the lookup table (normalize energy to 1)
759  79      dimensionality = tot_channel[0]
760  80      tot_constellation = tot_channel[1]
761  81      N0=pow(10.0,-esn0_db/10.0); # noise variance
762  82      if len(tot_constellation)/dimensionality != f.O():
763  83          sys.stderr.write (&apos;Incompatible FSM output cardinality and lookup table size.\n&apos;)
764  84          sys.exit (1)
765 </programlisting>
766
767
768
769 <para>
770 Then, "run_test" is called with a different "seed" so that we get
771 different data and noise realizations.
772 </para>
773 <programlisting>
774  91          (s,e)=run_test(f,Kb,bitspersymbol,K,channel,modulation,dimensionality,tot_constellation,N0,-long(666+i)) # run experiment with different seed to get different data and noise realizations
775 </programlisting>
776
777
778
779 <para>
780 Let us examine now the "run_test" function. 
781 First we set up the transmitter blocks.
782 We generate a packet of K random symbols and add a head and a tail of L zeros, 
783 L being the channel length. This is sufficient to drive the initial and final states to 0.
784 </para>
785 <programlisting>
786  14      L = len(channel)
787  15  
788  16      # TX
789  17      # this for loop is TOO slow in python!!!
790  18      packet = [0]*(K+2*L)
791  19      random.seed(seed)
792  20      for i in range(len(packet)):
793  21          packet[i] = random.randint(0, 2**bitspersymbol - 1) # random symbols
794  22      for i in range(L): # first/last L symbols set to 0
795  23          packet[i] = 0
796  24          packet[len(packet)-i-1] = 0
797  25      src = gr.vector_source_s(packet,False)
798  26      mod = gr.chunks_to_symbols_sf(modulation[1],modulation[0])
799 </programlisting>
800
801
802 <para>
803 The modulated symbols are filtered by the ISI channel and AWGN with appropriate noise variance is added.
804 </para>
805 <programlisting>
806  28      # CHANNEL
807  29      isi = gr.fir_filter_fff(1,channel)
808  30      add = gr.add_ff()
809  31      noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
810 </programlisting>
811
812
813
814 <para>
815 Since the output alphabet of the equivalent FSM is quite large (1024) we chose to utilize the combined metrics calculator and Viterbi algorithm block.
816 also note that the first L observations are irrelevant and tus can be skipped.
817 </para>
818 <programlisting>
819  33      # RX
820  34      skip = gr.skiphead(gr.sizeof_float, L) # skip the first L samples since you know they are coming from the L zero symbols
821  35      #metrics = trellis.metrics_f(f.O(),dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # data preprocessing to generate metrics for Viterbi
822  36      #va = trellis.viterbi_s(f,K+L,0,0) # Put -1 if the Initial/Final states are not set.
823  37      va = trellis.viterbi_combined_s(f,K+L,0,0,dimensionality,tot_constellation,trellis.TRELLIS_EUCLIDEAN) # using viterbi_combined_s instead of metrics_f/viterbi_s allows larger packet lengths because metrics_f is complaining for not being able to allocate large buffers. This is due to the large f.O() in this application...                                    
824  38      dst = gr.vector_sink_s()
825 </programlisting>
826
827 <para>
828 Now the VA can run once it is supplied by the initial and final states.
829 In this example both the initial and final states are set to 0.
830 The VA outputs the estimates of the input symbols which
831 are then compared with the transmitted sequence.
832 </para>
833 <programlisting>
834  49      data = dst.data()
835  50      ntotal = len(data) - L
836  51      nright=0
837  52      for i in range(ntotal):
838  53          if packet[i+L]==data[i]:
839  54              nright=nright+1
840  55          #else:
841  56              #print &quot;Error in &quot;, i
842 </programlisting>
843
844
845 <para>
846 The function returns the number of symbols and the number of symbols in error. Observe that this way the estimated error rate refers to 
847 2-bit-symbol error rate.
848 </para>
849 <programlisting>
850  58      return (ntotal,ntotal-nright)
851 </programlisting>
852
853
854
855 </sect1>
856
857
858
859
860
861
862 <!--=====================================================-->
863 <sect1 id="turbo"><title>Support for Concatenated Coding and Turbo Decoding</title>
864
865 <para>
866 To do...
867 </para>
868
869
870
871
872 </sect1>
873
874
875
876
877
878
879
880
881 <!--====================n================================-->
882 <sect1 id="future"><title>Future Work</title>
883
884
885
886 <itemizedlist>
887
888 <listitem>
889 <para>
890 Improve the documentation :-)
891 </para>
892 </listitem>
893
894 <listitem>
895 <para>
896 automate fsm generation from rational functions
897 (feedback form).
898 </para>
899 </listitem>
900
901 <listitem>
902 <para>
903 Optimize the VA code if possible.
904 </para>
905 </listitem>
906
907 <listitem>
908 <para>
909 A host of suboptimal
910 decoders, eg, sphere decoding, M- and T- algorithms, sequential decoding, etc.
911 can be implemented.
912 </para>
913 </listitem>
914
915
916 <listitem>
917 <para>
918 Although turbo decoding is rediculously slow in software, 
919 we can design it in principle. One question is, whether we should 
920 use the encoder, and SISO blocks and connect them
921 through GNU radio or we should implement turbo-decoding
922 as a single block (issues with buffering between blocks).
923 So far the former has been implemented.
924 </para>
925 </listitem>
926
927 </itemizedlist>
928
929 </sect1>
930
931
932 </article>