~ubuntu-branches/ubuntu/trusty/gnuradio/trusty-updates

« back to all changes in this revision

Viewing changes to gr-trellis/src/lib/core_algorithms.cc

  • Committer: Package Import Robot
  • Author(s): A. Maitland Bottoms
  • Date: 2012-02-26 21:26:16 UTC
  • mfrom: (1.1.4)
  • Revision ID: package-import@ubuntu.com-20120226212616-vsfkbi1158xshdql
Tags: 3.5.1-1
* new upstream version, re-packaged from scratch with modern tools
    closes: #642716, #645332, #394849, #616832, #590048, #642580,
    #647018, #557050, #559640, #631863
* CMake build

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* -*- c++ -*- */
 
2
/*
 
3
 * Copyright 2004 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 <cstring>
 
24
#include <stdexcept>
 
25
//#include <cstdio>
 
26
#include <iostream>
 
27
#include "core_algorithms.h"
 
28
#include "calc_metric.h"
 
29
 
 
30
static const float INF = 1.0e9;
 
31
 
 
32
float min(float a, float b)
 
33
{
 
34
  return a <= b ? a : b;
 
35
}
 
36
 
 
37
float min_star(float a, float b)
 
38
{
 
39
  return (a <= b ? a : b)-log(1+exp(a <= b ? a-b : b-a));
 
40
}
 
41
 
 
42
 
 
43
 
 
44
 
 
45
template <class T> 
 
46
void viterbi_algorithm(int I, int S, int O, 
 
47
             const std::vector<int> &NS,
 
48
             const std::vector<int> &OS,
 
49
             const std::vector< std::vector<int> > &PS,
 
50
             const std::vector< std::vector<int> > &PI,
 
51
             int K,
 
52
             int S0,int SK,
 
53
             const float *in, T *out)//,
 
54
             //std::vector<int> &trace) 
 
55
{
 
56
  std::vector<int> trace(S*K);
 
57
  std::vector<float> alpha(S*2);
 
58
  int alphai;
 
59
  float norm,mm,minm;
 
60
  int minmi;
 
61
  int st;
 
62
 
 
63
 
 
64
  if(S0<0) { // initial state not specified
 
65
      for(int i=0;i<S;i++) alpha[0*S+i]=0;
 
66
  }
 
67
  else {
 
68
      for(int i=0;i<S;i++) alpha[0*S+i]=INF;
 
69
      alpha[0*S+S0]=0.0;
 
70
  }
 
71
 
 
72
  alphai=0;
 
73
  for(int k=0;k<K;k++) {
 
74
      norm=INF;
 
75
      for(int j=0;j<S;j++) { // for each next state do ACS
 
76
          minm=INF;
 
77
          minmi=0;
 
78
          for(unsigned int i=0;i<PS[j].size();i++) {
 
79
              //int i0 = j*I+i;
 
80
              if((mm=alpha[alphai*S+PS[j][i]]+in[k*O+OS[PS[j][i]*I+PI[j][i]]])<minm)
 
81
                  minm=mm,minmi=i;
 
82
          }
 
83
          trace[k*S+j]=minmi;
 
84
          alpha[((alphai+1)%2)*S+j]=minm;
 
85
          if(minm<norm) norm=minm;
 
86
      }
 
87
      for(int j=0;j<S;j++) 
 
88
          alpha[((alphai+1)%2)*S+j]-=norm; // normalize total metrics so they do not explode
 
89
      alphai=(alphai+1)%2;
 
90
  }
 
91
 
 
92
  if(SK<0) { // final state not specified
 
93
      minm=INF;
 
94
      minmi=0;
 
95
      for(int i=0;i<S;i++)
 
96
          if((mm=alpha[alphai*S+i])<minm) minm=mm,minmi=i;
 
97
      st=minmi;
 
98
  }
 
99
  else {
 
100
      st=SK;
 
101
  }
 
102
 
 
103
  for(int k=K-1;k>=0;k--) { // traceback
 
104
      int i0=trace[k*S+st];
 
105
      out[k]= (T) PI[st][i0];
 
106
      st=PS[st][i0];
 
107
  }
 
108
 
 
109
}
 
110
 
 
111
 
 
112
template
 
113
void viterbi_algorithm<unsigned char>(int I, int S, int O,
 
114
             const std::vector<int> &NS,
 
115
             const std::vector<int> &OS,
 
116
             const std::vector< std::vector<int> > &PS,
 
117
             const std::vector< std::vector<int> > &PI,
 
118
             int K,
 
119
             int S0,int SK,
 
120
             const float *in, unsigned char *out);
 
121
 
 
122
 
 
123
template
 
124
void viterbi_algorithm<short>(int I, int S, int O,
 
125
             const std::vector<int> &NS,
 
126
             const std::vector<int> &OS,
 
127
             const std::vector< std::vector<int> > &PS,
 
128
             const std::vector< std::vector<int> > &PI,
 
129
             int K,
 
130
             int S0,int SK,
 
131
             const float *in, short *out);
 
132
 
 
133
template
 
134
void viterbi_algorithm<int>(int I, int S, int O,
 
135
             const std::vector<int> &NS,
 
136
             const std::vector<int> &OS,
 
137
             const std::vector< std::vector<int> > &PS,
 
138
             const std::vector< std::vector<int> > &PI,
 
139
             int K,
 
140
             int S0,int SK,
 
141
             const float *in, int *out);
 
142
 
 
143
 
 
144
 
 
145
//==============================================
 
146
 
 
147
template <class Ti, class To>
 
148
void viterbi_algorithm_combined(int I, int S, int O,
 
149
             const std::vector<int> &NS,
 
150
             const std::vector<int> &OS,
 
151
             const std::vector< std::vector<int> > &PS,
 
152
             const std::vector< std::vector<int> > &PI,
 
153
             int K,
 
154
             int S0,int SK,
 
155
             int D,
 
156
             const std::vector<Ti> &TABLE,
 
157
             trellis_metric_type_t TYPE,
 
158
             const Ti *in, To *out
 
159
)
 
160
{
 
161
  std::vector<int> trace(S*K);
 
162
  std::vector<float> alpha(S*2);
 
163
  float *metric = new float[O];
 
164
  int alphai;
 
165
  float norm,mm,minm;
 
166
  int minmi;
 
167
  int st;
 
168
 
 
169
  if(S0<0) { // initial state not specified
 
170
      for(int i=0;i<S;i++) alpha[0*S+i]=0;
 
171
  }
 
172
  else {
 
173
      for(int i=0;i<S;i++) alpha[0*S+i]=INF;
 
174
      alpha[0*S+S0]=0.0;
 
175
  }
 
176
 
 
177
  alphai=0;
 
178
  for(int k=0;k<K;k++) {
 
179
      calc_metric(O, D, TABLE, &(in[k*D]), metric,TYPE); // calc metrics
 
180
      norm=INF;
 
181
      for(int j=0;j<S;j++) { // for each next state do ACS
 
182
          minm=INF;
 
183
          minmi=0;
 
184
          for(unsigned int i=0;i<PS[j].size();i++) {
 
185
              //int i0 = j*I+i;
 
186
              if((mm=alpha[alphai*S+PS[j][i]]+metric[OS[PS[j][i]*I+PI[j][i]]])<minm)
 
187
                  minm=mm,minmi=i;
 
188
          }
 
189
          trace[k*S+j]=minmi;
 
190
          alpha[((alphai+1)%2)*S+j]=minm;
 
191
          if(minm<norm) norm=minm;
 
192
      }
 
193
      for(int j=0;j<S;j++) 
 
194
          alpha[((alphai+1)%2)*S+j]-=norm; // normalize total metrics so they do not explode
 
195
      alphai=(alphai+1)%2;
 
196
  }
 
197
 
 
198
  if(SK<0) { // final state not specified
 
199
      minm=INF;
 
200
      minmi=0;
 
201
      for(int i=0;i<S;i++)
 
202
          if((mm=alpha[alphai*S+i])<minm) minm=mm,minmi=i;
 
203
      st=minmi;
 
204
  }
 
205
  else {
 
206
      st=SK;
 
207
  }
 
208
 
 
209
  for(int k=K-1;k>=0;k--) { // traceback
 
210
      int i0=trace[k*S+st];
 
211
      out[k]= (To) PI[st][i0];
 
212
      st=PS[st][i0];
 
213
  }
 
214
  
 
215
  delete [] metric;
 
216
 
 
217
}
 
218
 
 
219
// Ti = s i f c
 
220
// To = b s i
 
221
 
 
222
//---------------
 
223
 
 
224
template
 
225
void viterbi_algorithm_combined<short,unsigned char>(int I, int S, int O,
 
226
             const std::vector<int> &NS,
 
227
             const std::vector<int> &OS,
 
228
             const std::vector< std::vector<int> > &PS,
 
229
             const std::vector< std::vector<int> > &PI,
 
230
             int K,
 
231
             int S0,int SK,
 
232
             int D,
 
233
             const std::vector<short> &TABLE,
 
234
             trellis_metric_type_t TYPE,
 
235
             const short *in, unsigned char *out
 
236
);
 
237
 
 
238
template
 
239
void viterbi_algorithm_combined<int,unsigned char>(int I, int S, int O,
 
240
             const std::vector<int> &NS,
 
241
             const std::vector<int> &OS,
 
242
             const std::vector< std::vector<int> > &PS,
 
243
             const std::vector< std::vector<int> > &PI,
 
244
             int K,
 
245
             int S0,int SK,
 
246
             int D,
 
247
             const std::vector<int> &TABLE,
 
248
             trellis_metric_type_t TYPE,
 
249
             const int *in, unsigned char *out
 
250
);
 
251
 
 
252
template
 
253
void viterbi_algorithm_combined<float,unsigned char>(int I, int S, int O,
 
254
             const std::vector<int> &NS,
 
255
             const std::vector<int> &OS,
 
256
             const std::vector< std::vector<int> > &PS,
 
257
             const std::vector< std::vector<int> > &PI,
 
258
             int K,
 
259
             int S0,int SK,
 
260
             int D,
 
261
             const std::vector<float> &TABLE,
 
262
             trellis_metric_type_t TYPE,
 
263
             const float *in, unsigned char *out
 
264
);
 
265
 
 
266
template
 
267
void viterbi_algorithm_combined<gr_complex,unsigned char>(int I, int S, int O,
 
268
             const std::vector<int> &NS,
 
269
             const std::vector<int> &OS,
 
270
             const std::vector< std::vector<int> > &PS,
 
271
             const std::vector< std::vector<int> > &PI,
 
272
             int K,
 
273
             int S0,int SK,
 
274
             int D,
 
275
             const std::vector<gr_complex> &TABLE,
 
276
             trellis_metric_type_t TYPE,
 
277
             const gr_complex *in, unsigned char *out
 
278
);
 
279
 
 
280
//---------------
 
281
 
 
282
template
 
283
void viterbi_algorithm_combined<short,short>(int I, int S, int O,
 
284
             const std::vector<int> &NS,
 
285
             const std::vector<int> &OS,
 
286
             const std::vector< std::vector<int> > &PS,
 
287
             const std::vector< std::vector<int> > &PI,
 
288
             int K,
 
289
             int S0,int SK,
 
290
             int D,
 
291
             const std::vector<short> &TABLE,
 
292
             trellis_metric_type_t TYPE,
 
293
             const short *in, short *out
 
294
);
 
295
 
 
296
template
 
297
void viterbi_algorithm_combined<int,short>(int I, int S, int O,
 
298
             const std::vector<int> &NS,
 
299
             const std::vector<int> &OS,
 
300
             const std::vector< std::vector<int> > &PS,
 
301
             const std::vector< std::vector<int> > &PI,
 
302
             int K,
 
303
             int S0,int SK,
 
304
             int D,
 
305
             const std::vector<int> &TABLE,
 
306
             trellis_metric_type_t TYPE,
 
307
             const int *in, short *out
 
308
);
 
309
 
 
310
template
 
311
void viterbi_algorithm_combined<float,short>(int I, int S, int O,
 
312
             const std::vector<int> &NS,
 
313
             const std::vector<int> &OS,
 
314
             const std::vector< std::vector<int> > &PS,
 
315
             const std::vector< std::vector<int> > &PI,
 
316
             int K,
 
317
             int S0,int SK,
 
318
             int D,
 
319
             const std::vector<float> &TABLE,
 
320
             trellis_metric_type_t TYPE,
 
321
             const float *in, short *out
 
322
);
 
323
 
 
324
template
 
325
void viterbi_algorithm_combined<gr_complex,short>(int I, int S, int O,
 
326
             const std::vector<int> &NS,
 
327
             const std::vector<int> &OS,
 
328
             const std::vector< std::vector<int> > &PS,
 
329
             const std::vector< std::vector<int> > &PI,
 
330
             int K,
 
331
             int S0,int SK,
 
332
             int D,
 
333
             const std::vector<gr_complex> &TABLE,
 
334
             trellis_metric_type_t TYPE,
 
335
             const gr_complex *in, short *out
 
336
);
 
337
 
 
338
//--------------
 
339
 
 
340
template
 
341
void viterbi_algorithm_combined<short,int>(int I, int S, int O,
 
342
             const std::vector<int> &NS,
 
343
             const std::vector<int> &OS,
 
344
             const std::vector< std::vector<int> > &PS,
 
345
             const std::vector< std::vector<int> > &PI,
 
346
             int K,
 
347
             int S0,int SK,
 
348
             int D,
 
349
             const std::vector<short> &TABLE,
 
350
             trellis_metric_type_t TYPE,
 
351
             const short *in, int *out
 
352
);
 
353
 
 
354
template
 
355
void viterbi_algorithm_combined<int,int>(int I, int S, int O,
 
356
             const std::vector<int> &NS,
 
357
             const std::vector<int> &OS,
 
358
             const std::vector< std::vector<int> > &PS,
 
359
             const std::vector< std::vector<int> > &PI,
 
360
             int K,
 
361
             int S0,int SK,
 
362
             int D,
 
363
             const std::vector<int> &TABLE,
 
364
             trellis_metric_type_t TYPE,
 
365
             const int *in, int *out
 
366
);
 
367
 
 
368
template
 
369
void viterbi_algorithm_combined<float,int>(int I, int S, int O,
 
370
             const std::vector<int> &NS,
 
371
             const std::vector<int> &OS,
 
372
             const std::vector< std::vector<int> > &PS,
 
373
             const std::vector< std::vector<int> > &PI,
 
374
             int K,
 
375
             int S0,int SK,
 
376
             int D,
 
377
             const std::vector<float> &TABLE,
 
378
             trellis_metric_type_t TYPE,
 
379
             const float *in, int *out
 
380
);
 
381
 
 
382
template
 
383
void viterbi_algorithm_combined<gr_complex,int>(int I, int S, int O,
 
384
             const std::vector<int> &NS,
 
385
             const std::vector<int> &OS,
 
386
             const std::vector< std::vector<int> > &PS,
 
387
             const std::vector< std::vector<int> > &PI,
 
388
             int K,
 
389
             int S0,int SK,
 
390
             int D,
 
391
             const std::vector<gr_complex> &TABLE,
 
392
             trellis_metric_type_t TYPE,
 
393
             const gr_complex *in, int *out
 
394
);
 
395
 
 
396
 
 
397
 
 
398
 
 
399
 
 
400
 
 
401
 
 
402
 
 
403
 
 
404
 
 
405
 
 
406
 
 
407
 
 
408
 
 
409
 
 
410
 
 
411
 
 
412
 
 
413
 
 
414
 
 
415
//===============================================
 
416
 
 
417
 
 
418
void siso_algorithm(int I, int S, int O, 
 
419
             const std::vector<int> &NS,
 
420
             const std::vector<int> &OS,
 
421
             const std::vector< std::vector<int> > &PS,
 
422
             const std::vector< std::vector<int> > &PI,
 
423
             int K,
 
424
             int S0,int SK,
 
425
             bool POSTI, bool POSTO,
 
426
             float (*p2mymin)(float,float),
 
427
             const float *priori, const float *prioro, float *post//,
 
428
             //std::vector<float> &alpha,
 
429
             //std::vector<float> &beta
 
430
             ) 
 
431
{
 
432
  float norm,mm,minm;
 
433
  std::vector<float> alpha(S*(K+1));
 
434
  std::vector<float> beta(S*(K+1));
 
435
 
 
436
 
 
437
  if(S0<0) { // initial state not specified
 
438
      for(int i=0;i<S;i++) alpha[0*S+i]=0;
 
439
  }
 
440
  else {
 
441
      for(int i=0;i<S;i++) alpha[0*S+i]=INF;
 
442
      alpha[0*S+S0]=0.0;
 
443
  }
 
444
 
 
445
  for(int k=0;k<K;k++) { // forward recursion
 
446
      norm=INF;
 
447
      for(int j=0;j<S;j++) {
 
448
          minm=INF;
 
449
          for(unsigned int i=0;i<PS[j].size();i++) {
 
450
              //int i0 = j*I+i;
 
451
              mm=alpha[k*S+PS[j][i]]+priori[k*I+PI[j][i]]+prioro[k*O+OS[PS[j][i]*I+PI[j][i]]];
 
452
              minm=(*p2mymin)(minm,mm);
 
453
          }
 
454
          alpha[(k+1)*S+j]=minm;
 
455
          if(minm<norm) norm=minm;
 
456
      }
 
457
      for(int j=0;j<S;j++) 
 
458
          alpha[(k+1)*S+j]-=norm; // normalize total metrics so they do not explode
 
459
  }
 
460
 
 
461
  if(SK<0) { // final state not specified
 
462
      for(int i=0;i<S;i++) beta[K*S+i]=0;
 
463
  }
 
464
  else {
 
465
      for(int i=0;i<S;i++) beta[K*S+i]=INF;
 
466
      beta[K*S+SK]=0.0;
 
467
  }
 
468
 
 
469
  for(int k=K-1;k>=0;k--) { // backward recursion
 
470
      norm=INF;
 
471
      for(int j=0;j<S;j++) { 
 
472
          minm=INF;
 
473
          for(int i=0;i<I;i++) {
 
474
              int i0 = j*I+i;
 
475
              mm=beta[(k+1)*S+NS[i0]]+priori[k*I+i]+prioro[k*O+OS[i0]];
 
476
              minm=(*p2mymin)(minm,mm);
 
477
          }
 
478
          beta[k*S+j]=minm;
 
479
          if(minm<norm) norm=minm;
 
480
      }
 
481
      for(int j=0;j<S;j++)
 
482
          beta[k*S+j]-=norm; // normalize total metrics so they do not explode
 
483
  }
 
484
 
 
485
 
 
486
if (POSTI && POSTO)
 
487
{
 
488
  for(int k=0;k<K;k++) { // input combining
 
489
      norm=INF;
 
490
      for(int i=0;i<I;i++) {
 
491
          minm=INF;
 
492
          for(int j=0;j<S;j++) {
 
493
              mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 
494
              minm=(*p2mymin)(minm,mm);
 
495
          }
 
496
          post[k*(I+O)+i]=minm;
 
497
          if(minm<norm) norm=minm;
 
498
      }
 
499
      for(int i=0;i<I;i++)
 
500
          post[k*(I+O)+i]-=norm; // normalize metrics
 
501
  }
 
502
 
 
503
 
 
504
  for(int k=0;k<K;k++) { // output combining
 
505
      norm=INF;
 
506
      for(int n=0;n<O;n++) {
 
507
          minm=INF;
 
508
          for(int j=0;j<S;j++) {
 
509
              for(int i=0;i<I;i++) {
 
510
                  mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 
511
                  minm=(*p2mymin)(minm,mm);
 
512
              }
 
513
          }
 
514
          post[k*(I+O)+I+n]=minm;
 
515
          if(minm<norm) norm=minm;
 
516
      }
 
517
      for(int n=0;n<O;n++)
 
518
          post[k*(I+O)+I+n]-=norm; // normalize metrics
 
519
  }
 
520
 
521
else if(POSTI) 
 
522
{
 
523
  for(int k=0;k<K;k++) { // input combining
 
524
      norm=INF;
 
525
      for(int i=0;i<I;i++) {
 
526
          minm=INF;
 
527
          for(int j=0;j<S;j++) {
 
528
              mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 
529
              minm=(*p2mymin)(minm,mm);
 
530
          }
 
531
          post[k*I+i]=minm;
 
532
          if(minm<norm) norm=minm;
 
533
      }
 
534
      for(int i=0;i<I;i++)
 
535
          post[k*I+i]-=norm; // normalize metrics
 
536
  }
 
537
}
 
538
else if(POSTO)
 
539
{
 
540
  for(int k=0;k<K;k++) { // output combining
 
541
      norm=INF;
 
542
      for(int n=0;n<O;n++) {
 
543
          minm=INF;
 
544
          for(int j=0;j<S;j++) {
 
545
              for(int i=0;i<I;i++) {
 
546
                  mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 
547
                  minm=(*p2mymin)(minm,mm);
 
548
              }
 
549
          }
 
550
          post[k*O+n]=minm;
 
551
          if(minm<norm) norm=minm;
 
552
      }
 
553
      for(int n=0;n<O;n++)
 
554
          post[k*O+n]-=norm; // normalize metrics
 
555
  }
 
556
}
 
557
else
 
558
    throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 
559
 
 
560
}
 
561
 
 
562
 
 
563
//===========================================================
 
564
 
 
565
template <class T>
 
566
void siso_algorithm_combined(int I, int S, int O, 
 
567
             const std::vector<int> &NS,
 
568
             const std::vector<int> &OS,
 
569
             const std::vector< std::vector<int> > &PS,
 
570
             const std::vector< std::vector<int> > &PI,
 
571
             int K,
 
572
             int S0,int SK,
 
573
             bool POSTI, bool POSTO,
 
574
             float (*p2mymin)(float,float),
 
575
             int D,
 
576
             const std::vector<T> &TABLE,
 
577
             trellis_metric_type_t TYPE,
 
578
             const float *priori, const T *observations, float *post
 
579
 
580
{
 
581
  float norm,mm,minm;
 
582
  std::vector<float> alpha(S*(K+1));
 
583
  std::vector<float> beta(S*(K+1));
 
584
  float *prioro = new float[O*K];
 
585
 
 
586
 
 
587
  if(S0<0) { // initial state not specified
 
588
      for(int i=0;i<S;i++) alpha[0*S+i]=0;
 
589
  }
 
590
  else {
 
591
      for(int i=0;i<S;i++) alpha[0*S+i]=INF;
 
592
      alpha[0*S+S0]=0.0;
 
593
  }
 
594
 
 
595
  for(int k=0;k<K;k++) { // forward recursion
 
596
      calc_metric(O, D, TABLE, &(observations[k*D]), &(prioro[k*O]),TYPE); // calc metrics
 
597
      norm=INF;
 
598
      for(int j=0;j<S;j++) {
 
599
          minm=INF;
 
600
          for(unsigned int i=0;i<PS[j].size();i++) {
 
601
              //int i0 = j*I+i;
 
602
              mm=alpha[k*S+PS[j][i]]+priori[k*I+PI[j][i]]+prioro[k*O+OS[PS[j][i]*I+PI[j][i]]];
 
603
              minm=(*p2mymin)(minm,mm);
 
604
          }
 
605
          alpha[(k+1)*S+j]=minm;
 
606
          if(minm<norm) norm=minm;
 
607
      }
 
608
      for(int j=0;j<S;j++) 
 
609
          alpha[(k+1)*S+j]-=norm; // normalize total metrics so they do not explode
 
610
  }
 
611
 
 
612
  if(SK<0) { // final state not specified
 
613
      for(int i=0;i<S;i++) beta[K*S+i]=0;
 
614
  }
 
615
  else {
 
616
      for(int i=0;i<S;i++) beta[K*S+i]=INF;
 
617
      beta[K*S+SK]=0.0;
 
618
  }
 
619
 
 
620
  for(int k=K-1;k>=0;k--) { // backward recursion
 
621
      norm=INF;
 
622
      for(int j=0;j<S;j++) { 
 
623
          minm=INF;
 
624
          for(int i=0;i<I;i++) {
 
625
              int i0 = j*I+i;
 
626
              mm=beta[(k+1)*S+NS[i0]]+priori[k*I+i]+prioro[k*O+OS[i0]];
 
627
              minm=(*p2mymin)(minm,mm);
 
628
          }
 
629
          beta[k*S+j]=minm;
 
630
          if(minm<norm) norm=minm;
 
631
      }
 
632
      for(int j=0;j<S;j++)
 
633
          beta[k*S+j]-=norm; // normalize total metrics so they do not explode
 
634
  }
 
635
 
 
636
 
 
637
  if (POSTI && POSTO)
 
638
  {
 
639
    for(int k=0;k<K;k++) { // input combining
 
640
        norm=INF;
 
641
        for(int i=0;i<I;i++) {
 
642
            minm=INF;
 
643
            for(int j=0;j<S;j++) {
 
644
                mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 
645
                minm=(*p2mymin)(minm,mm);
 
646
            }
 
647
            post[k*(I+O)+i]=minm;
 
648
            if(minm<norm) norm=minm;
 
649
        }
 
650
        for(int i=0;i<I;i++)
 
651
            post[k*(I+O)+i]-=norm; // normalize metrics
 
652
    }
 
653
 
 
654
 
 
655
    for(int k=0;k<K;k++) { // output combining
 
656
        norm=INF;
 
657
        for(int n=0;n<O;n++) {
 
658
            minm=INF;
 
659
            for(int j=0;j<S;j++) {
 
660
                for(int i=0;i<I;i++) {
 
661
                    mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 
662
                    minm=(*p2mymin)(minm,mm);
 
663
                }
 
664
            }
 
665
            post[k*(I+O)+I+n]=minm;
 
666
            if(minm<norm) norm=minm;
 
667
        }
 
668
        for(int n=0;n<O;n++)
 
669
            post[k*(I+O)+I+n]-=norm; // normalize metrics
 
670
    }
 
671
  } 
 
672
  else if(POSTI) 
 
673
  {
 
674
    for(int k=0;k<K;k++) { // input combining
 
675
        norm=INF;
 
676
        for(int i=0;i<I;i++) {
 
677
            minm=INF;
 
678
            for(int j=0;j<S;j++) {
 
679
                mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 
680
                minm=(*p2mymin)(minm,mm);
 
681
            }
 
682
            post[k*I+i]=minm;
 
683
            if(minm<norm) norm=minm;
 
684
        }
 
685
        for(int i=0;i<I;i++)
 
686
            post[k*I+i]-=norm; // normalize metrics
 
687
    }
 
688
  }
 
689
  else if(POSTO)
 
690
  {
 
691
    for(int k=0;k<K;k++) { // output combining
 
692
        norm=INF;
 
693
        for(int n=0;n<O;n++) {
 
694
            minm=INF;
 
695
            for(int j=0;j<S;j++) {
 
696
                for(int i=0;i<I;i++) {
 
697
                    mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 
698
                    minm=(*p2mymin)(minm,mm);
 
699
                }
 
700
            }
 
701
            post[k*O+n]=minm;
 
702
            if(minm<norm) norm=minm;
 
703
        }
 
704
        for(int n=0;n<O;n++)
 
705
            post[k*O+n]-=norm; // normalize metrics
 
706
    }
 
707
  }
 
708
  else
 
709
    throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 
710
 
 
711
  delete [] prioro;
 
712
 
 
713
}
 
714
 
 
715
//---------
 
716
 
 
717
template
 
718
void siso_algorithm_combined<short>(int I, int S, int O,
 
719
             const std::vector<int> &NS,
 
720
             const std::vector<int> &OS,
 
721
             const std::vector< std::vector<int> > &PS,
 
722
             const std::vector< std::vector<int> > &PI,
 
723
             int K,
 
724
             int S0,int SK,
 
725
             bool POSTI, bool POSTO,
 
726
             float (*p2mymin)(float,float),
 
727
             int D,
 
728
             const std::vector<short> &TABLE,
 
729
             trellis_metric_type_t TYPE,
 
730
             const float *priori, const short *observations, float *post
 
731
);
 
732
 
 
733
template
 
734
void siso_algorithm_combined<int>(int I, int S, int O,
 
735
             const std::vector<int> &NS,
 
736
             const std::vector<int> &OS,
 
737
             const std::vector< std::vector<int> > &PS,
 
738
             const std::vector< std::vector<int> > &PI,
 
739
             int K,
 
740
             int S0,int SK,
 
741
             bool POSTI, bool POSTO,
 
742
             float (*p2mymin)(float,float),
 
743
             int D,
 
744
             const std::vector<int> &TABLE,
 
745
             trellis_metric_type_t TYPE,
 
746
             const float *priori, const int *observations, float *post
 
747
);
 
748
 
 
749
template
 
750
void siso_algorithm_combined<float>(int I, int S, int O,
 
751
             const std::vector<int> &NS,
 
752
             const std::vector<int> &OS,
 
753
             const std::vector< std::vector<int> > &PS,
 
754
             const std::vector< std::vector<int> > &PI,
 
755
             int K,
 
756
             int S0,int SK,
 
757
             bool POSTI, bool POSTO,
 
758
             float (*p2mymin)(float,float),
 
759
             int D,
 
760
             const std::vector<float> &TABLE,
 
761
             trellis_metric_type_t TYPE,
 
762
             const float *priori, const float *observations, float *post
 
763
);
 
764
 
 
765
template
 
766
void siso_algorithm_combined<gr_complex>(int I, int S, int O,
 
767
             const std::vector<int> &NS,
 
768
             const std::vector<int> &OS,
 
769
             const std::vector< std::vector<int> > &PS,
 
770
             const std::vector< std::vector<int> > &PI,
 
771
             int K,
 
772
             int S0,int SK,
 
773
             bool POSTI, bool POSTO,
 
774
             float (*p2mymin)(float,float),
 
775
             int D,
 
776
             const std::vector<gr_complex> &TABLE,
 
777
             trellis_metric_type_t TYPE,
 
778
             const float *priori, const gr_complex *observations, float *post
 
779
);
 
780
 
 
781
//=========================================================
 
782
 
 
783
template<class Ti, class To>
 
784
void sccc_decoder_combined(
 
785
      const fsm &FSMo, int STo0, int SToK,
 
786
      const fsm &FSMi, int STi0, int STiK,
 
787
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
788
      float (*p2mymin)(float,float),
 
789
      int D, const std::vector<Ti> &TABLE,
 
790
      trellis_metric_type_t METRIC_TYPE,
 
791
      float scaling,
 
792
      const Ti *observations, To *data
 
793
)
 
794
{
 
795
 
 
796
//allocate space for priori, prioro and posti of inner FSM
 
797
std::vector<float> ipriori(blocklength*FSMi.I(),0.0);
 
798
std::vector<float> iprioro(blocklength*FSMi.O());
 
799
std::vector<float> iposti(blocklength*FSMi.I());
 
800
 
 
801
//allocate space for priori, prioro and posto of outer FSM
 
802
std::vector<float> opriori(blocklength*FSMo.I(),0.0);
 
803
std::vector<float> oprioro(blocklength*FSMo.O());
 
804
std::vector<float> oposti(blocklength*FSMo.I());
 
805
std::vector<float> oposto(blocklength*FSMo.O());
 
806
 
 
807
// turn observations to neg-log-priors
 
808
for(int k=0;k<blocklength;k++) {
 
809
  calc_metric(FSMi.O(), D, TABLE, &(observations[k*D]), &(iprioro[k*FSMi.O()]),METRIC_TYPE);
 
810
  iprioro[k*FSMi.O()] *= scaling;
 
811
}
 
812
 
 
813
for(int rep=0;rep<iterations;rep++) {
 
814
  // run inner SISO
 
815
  siso_algorithm(FSMi.I(),FSMi.S(),FSMi.O(),
 
816
             FSMi.NS(), FSMi.OS(), FSMi.PS(), FSMi.PI(),
 
817
             blocklength,
 
818
             STi0,STiK,
 
819
             true, false,
 
820
             p2mymin,
 
821
             &(ipriori[0]),  &(iprioro[0]), &(iposti[0])
 
822
  );
 
823
 
 
824
  //interleave soft info inner -> outer
 
825
  for(int k=0;k<blocklength;k++) {
 
826
    int ki = INTERLEAVER.DEINTER()[k];
 
827
    //for(int i=0;i<FSMi.I();i++) {
 
828
      //oprioro[k*FSMi.I()+i]=iposti[ki*FSMi.I()+i];
 
829
    //}
 
830
    memcpy(&(oprioro[k*FSMi.I()]),&(iposti[ki*FSMi.I()]),FSMi.I()*sizeof(float));
 
831
  } 
 
832
 
 
833
  // run outer SISO
 
834
 
 
835
  if(rep<iterations-1) { // do not produce posti
 
836
    siso_algorithm(FSMo.I(),FSMo.S(),FSMo.O(),
 
837
             FSMo.NS(), FSMo.OS(), FSMo.PS(), FSMo.PI(),
 
838
             blocklength,
 
839
             STo0,SToK,
 
840
             false, true,
 
841
             p2mymin,
 
842
             &(opriori[0]),  &(oprioro[0]), &(oposto[0])
 
843
    );
 
844
 
 
845
    //interleave soft info outer --> inner
 
846
    for(int k=0;k<blocklength;k++) {
 
847
      int ki = INTERLEAVER.DEINTER()[k];
 
848
      //for(int i=0;i<FSMi.I();i++) {
 
849
        //ipriori[ki*FSMi.I()+i]=oposto[k*FSMi.I()+i];
 
850
      //}
 
851
      memcpy(&(ipriori[ki*FSMi.I()]),&(oposto[k*FSMi.I()]),FSMi.I()*sizeof(float));
 
852
    } 
 
853
  }
 
854
  else // produce posti but not posto
 
855
    
 
856
    siso_algorithm(FSMo.I(),FSMo.S(),FSMo.O(),
 
857
             FSMo.NS(), FSMo.OS(), FSMo.PS(), FSMo.PI(),
 
858
             blocklength,
 
859
             STo0,SToK,
 
860
             true, false,
 
861
             p2mymin,
 
862
             &(opriori[0]),  &(oprioro[0]), &(oposti[0])
 
863
    );
 
864
   
 
865
    /*
 
866
    viterbi_algorithm(FSMo.I(),FSMo.S(),FSMo.O(),
 
867
             FSMo.NS(), FSMo.OS(), FSMo.PS(), FSMo.PI(),
 
868
             blocklength,
 
869
             STo0,SToK,
 
870
             &(oprioro[0]), data
 
871
    );
 
872
    */
 
873
 
 
874
}
 
875
 
 
876
 
 
877
// generate hard decisions
 
878
for(int k=0;k<blocklength;k++) {
 
879
  float min=INF;
 
880
  int mini=0;
 
881
  for(int i=0;i<FSMo.I();i++) {
 
882
    if(oposti[k*FSMo.I()+i]<min) {
 
883
      min=oposti[k*FSMo.I()+i];
 
884
      mini=i;
 
885
    }
 
886
  }
 
887
  data[k]=(To)mini;
 
888
}
 
889
 
 
890
 
 
891
 
 
892
}
 
893
 
 
894
//-------
 
895
 
 
896
template
 
897
void sccc_decoder_combined<float,unsigned char>(
 
898
      const fsm &FSMo, int STo0, int SToK,
 
899
      const fsm &FSMi, int STi0, int STiK,
 
900
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
901
      float (*p2mymin)(float,float),
 
902
      int D, const std::vector<float> &TABLE,
 
903
      trellis_metric_type_t METRIC_TYPE,
 
904
      float scaling,
 
905
      const float *observations, unsigned char *data
 
906
);
 
907
 
 
908
template
 
909
void sccc_decoder_combined<float,short>(
 
910
      const fsm &FSMo, int STo0, int SToK,
 
911
      const fsm &FSMi, int STi0, int STiK,
 
912
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
913
      float (*p2mymin)(float,float),
 
914
      int D, const std::vector<float> &TABLE,
 
915
      trellis_metric_type_t METRIC_TYPE,
 
916
      float scaling,
 
917
      const float *observations, short *data
 
918
);
 
919
 
 
920
template
 
921
void sccc_decoder_combined<float,int>(
 
922
      const fsm &FSMo, int STo0, int SToK,
 
923
      const fsm &FSMi, int STi0, int STiK,
 
924
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
925
      float (*p2mymin)(float,float),
 
926
      int D, const std::vector<float> &TABLE,
 
927
      trellis_metric_type_t METRIC_TYPE,
 
928
      float scaling,
 
929
      const float *observations, int *data
 
930
);
 
931
 
 
932
template
 
933
void sccc_decoder_combined<gr_complex,unsigned char>(
 
934
      const fsm &FSMo, int STo0, int SToK,
 
935
      const fsm &FSMi, int STi0, int STiK,
 
936
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
937
      float (*p2mymin)(float,float),
 
938
      int D, const std::vector<gr_complex> &TABLE,
 
939
      trellis_metric_type_t METRIC_TYPE,
 
940
      float scaling,
 
941
      const gr_complex *observations, unsigned char *data
 
942
);
 
943
 
 
944
template
 
945
void sccc_decoder_combined<gr_complex,short>(
 
946
      const fsm &FSMo, int STo0, int SToK,
 
947
      const fsm &FSMi, int STi0, int STiK,
 
948
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
949
      float (*p2mymin)(float,float),
 
950
      int D, const std::vector<gr_complex> &TABLE,
 
951
      trellis_metric_type_t METRIC_TYPE,
 
952
      float scaling,
 
953
      const gr_complex *observations, short *data
 
954
);
 
955
 
 
956
template
 
957
void sccc_decoder_combined<gr_complex,int>(
 
958
      const fsm &FSMo, int STo0, int SToK,
 
959
      const fsm &FSMi, int STi0, int STiK,
 
960
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
961
      float (*p2mymin)(float,float),
 
962
      int D, const std::vector<gr_complex> &TABLE,
 
963
      trellis_metric_type_t METRIC_TYPE,
 
964
      float scaling,
 
965
      const gr_complex *observations, int *data
 
966
);
 
967
 
 
968
 
 
969
 
 
970
//=========================================================
 
971
 
 
972
template<class T>
 
973
void sccc_decoder(
 
974
      const fsm &FSMo, int STo0, int SToK,
 
975
      const fsm &FSMi, int STi0, int STiK,
 
976
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
977
      float (*p2mymin)(float,float),
 
978
      const float *iprioro, T *data
 
979
)
 
980
{
 
981
  //allocate space for priori, and posti of inner FSM
 
982
  std::vector<float> ipriori(blocklength*FSMi.I(),0.0);
 
983
  std::vector<float> iposti(blocklength*FSMi.I());
 
984
 
 
985
  //allocate space for priori, prioro and posto of outer FSM
 
986
  std::vector<float> opriori(blocklength*FSMo.I(),0.0);
 
987
  std::vector<float> oprioro(blocklength*FSMo.O());
 
988
  std::vector<float> oposti(blocklength*FSMo.I());
 
989
  std::vector<float> oposto(blocklength*FSMo.O());
 
990
 
 
991
  for(int rep=0;rep<iterations;rep++) {
 
992
    // run inner SISO
 
993
    siso_algorithm(FSMi.I(),FSMi.S(),FSMi.O(),
 
994
             FSMi.NS(), FSMi.OS(), FSMi.PS(), FSMi.PI(),
 
995
             blocklength,
 
996
             STi0,STiK,
 
997
             true, false,
 
998
             p2mymin,
 
999
             &(ipriori[0]),  &(iprioro[0]), &(iposti[0])
 
1000
    );
 
1001
 
 
1002
    //interleave soft info inner -> outer
 
1003
    for(int k=0;k<blocklength;k++) {
 
1004
      int ki = INTERLEAVER.DEINTER()[k];
 
1005
      //for(int i=0;i<FSMi.I();i++) {
 
1006
        //oprioro[k*FSMi.I()+i]=iposti[ki*FSMi.I()+i];
 
1007
      //}
 
1008
      memcpy(&(oprioro[k*FSMi.I()]),&(iposti[ki*FSMi.I()]),FSMi.I()*sizeof(float));
 
1009
    } 
 
1010
 
 
1011
    // run outer SISO
 
1012
 
 
1013
    if(rep<iterations-1) { // do not produce posti
 
1014
      siso_algorithm(FSMo.I(),FSMo.S(),FSMo.O(),
 
1015
             FSMo.NS(), FSMo.OS(), FSMo.PS(), FSMo.PI(),
 
1016
             blocklength,
 
1017
             STo0,SToK,
 
1018
             false, true,
 
1019
             p2mymin,
 
1020
             &(opriori[0]),  &(oprioro[0]), &(oposto[0])
 
1021
      );
 
1022
 
 
1023
      //interleave soft info outer --> inner
 
1024
      for(int k=0;k<blocklength;k++) {
 
1025
        int ki = INTERLEAVER.DEINTER()[k];
 
1026
        //for(int i=0;i<FSMi.I();i++) {
 
1027
          //ipriori[ki*FSMi.I()+i]=oposto[k*FSMi.I()+i];
 
1028
        //}
 
1029
        memcpy(&(ipriori[ki*FSMi.I()]),&(oposto[k*FSMi.I()]),FSMi.I()*sizeof(float));
 
1030
      } 
 
1031
    }
 
1032
    else {// produce posti but not posto
 
1033
    
 
1034
      siso_algorithm(FSMo.I(),FSMo.S(),FSMo.O(),
 
1035
             FSMo.NS(), FSMo.OS(), FSMo.PS(), FSMo.PI(),
 
1036
             blocklength,
 
1037
             STo0,SToK,
 
1038
             true, false,
 
1039
             p2mymin,
 
1040
             &(opriori[0]),  &(oprioro[0]), &(oposti[0])
 
1041
      );
 
1042
   
 
1043
      /*
 
1044
      viterbi_algorithm(FSMo.I(),FSMo.S(),FSMo.O(),
 
1045
             FSMo.NS(), FSMo.OS(), FSMo.PS(), FSMo.PI(),
 
1046
             blocklength,
 
1047
             STo0,SToK,
 
1048
             &(oprioro[0]), data
 
1049
      );
 
1050
      */
 
1051
    }
 
1052
 
 
1053
  } // end iterations
 
1054
 
 
1055
  // generate hard decisions
 
1056
  for(int k=0;k<blocklength;k++) {
 
1057
    float min=INF;
 
1058
    int mini=0;
 
1059
    for(int i=0;i<FSMo.I();i++) {
 
1060
      if(oposti[k*FSMo.I()+i]<min) {
 
1061
        min=oposti[k*FSMo.I()+i];
 
1062
        mini=i;
 
1063
      }
 
1064
    }
 
1065
    data[k]=(T)mini;
 
1066
  }
 
1067
 
 
1068
 
 
1069
 
 
1070
}
 
1071
 
 
1072
//-------
 
1073
 
 
1074
template
 
1075
void sccc_decoder<unsigned char>(
 
1076
      const fsm &FSMo, int STo0, int SToK,
 
1077
      const fsm &FSMi, int STi0, int STiK,
 
1078
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1079
      float (*p2mymin)(float,float),
 
1080
      const float *iprioro, unsigned char *data
 
1081
);
 
1082
 
 
1083
template
 
1084
void sccc_decoder<short>(
 
1085
      const fsm &FSMo, int STo0, int SToK,
 
1086
      const fsm &FSMi, int STi0, int STiK,
 
1087
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1088
      float (*p2mymin)(float,float),
 
1089
      const float *iprioro, short *data
 
1090
);
 
1091
 
 
1092
template
 
1093
void sccc_decoder<int>(
 
1094
      const fsm &FSMo, int STo0, int SToK,
 
1095
      const fsm &FSMi, int STi0, int STiK,
 
1096
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1097
      float (*p2mymin)(float,float),
 
1098
      const float *iprioro, int *data
 
1099
);
 
1100
 
 
1101
 
 
1102
//====================================================
 
1103
 
 
1104
template<class T>
 
1105
void pccc_decoder(
 
1106
      const fsm &FSM1, int ST10, int ST1K,
 
1107
      const fsm &FSM2, int ST20, int ST2K,
 
1108
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1109
      float (*p2mymin)(float,float),
 
1110
      const float *cprioro, T *data
 
1111
)
 
1112
{
 
1113
 
 
1114
  //allocate space for priori, prioro and posti of FSM1
 
1115
  std::vector<float> priori1(blocklength*FSM1.I(),0.0);
 
1116
  std::vector<float> prioro1(blocklength*FSM1.O());
 
1117
  std::vector<float> posti1(blocklength*FSM1.I());
 
1118
 
 
1119
  //allocate space for priori, prioro and posti of FSM2
 
1120
  std::vector<float> priori2(blocklength*FSM2.I(),0.0);
 
1121
  std::vector<float> prioro2(blocklength*FSM2.O());
 
1122
  std::vector<float> posti2(blocklength*FSM2.I());
 
1123
 
 
1124
  //generate prioro1,2 (metrics are not updated per iteration: this is not the best you can do...)
 
1125
  for (int k=0;k<blocklength;k++) {
 
1126
    //std::cout << k << std::endl;
 
1127
    for(int i=0;i<FSM1.O();i++) {
 
1128
      float x=cprioro[k*FSM1.O()*FSM2.O()+i*FSM1.O()+0];
 
1129
      for(int j=1;j<FSM2.O();j++)
 
1130
        x = (*p2mymin)(x,cprioro[k*FSM1.O()*FSM2.O()+i*FSM1.O()+j]);
 
1131
      prioro1[k*FSM1.O()+i]=x;
 
1132
      //std::cout <<  prioro1[k*FSM1.O()+i] << ", ";
 
1133
    }
 
1134
    //std::cout << std::endl;
 
1135
    for(int i=0;i<FSM2.O();i++) {
 
1136
      float x=cprioro[k*FSM1.O()*FSM2.O()+0*FSM1.O()+i];
 
1137
      for(int j=1;j<FSM1.O();j++)
 
1138
        x = (*p2mymin)(x,cprioro[k*FSM1.O()*FSM2.O()+j*FSM1.O()+i]);
 
1139
      prioro2[k*FSM2.O()+i]=x;
 
1140
    }
 
1141
  }
 
1142
 
 
1143
  for(int rep=0;rep<iterations;rep++) {
 
1144
    // run  SISO 1
 
1145
    siso_algorithm(FSM1.I(),FSM1.S(),FSM1.O(),
 
1146
             FSM1.NS(), FSM1.OS(), FSM1.PS(), FSM1.PI(),
 
1147
             blocklength,
 
1148
             ST10,ST1K,
 
1149
             true, false,
 
1150
             p2mymin,
 
1151
             &(priori1[0]),  &(prioro1[0]), &(posti1[0])
 
1152
    );
 
1153
   
 
1154
    //for(int k=0;k<blocklength;k++){
 
1155
      //for(int i=0;i<FSM1.I();i++)
 
1156
        //std::cout << posti1[k*FSM1.I()+i] << ", ";
 
1157
      //std::cout << std::endl;
 
1158
    //}
 
1159
 
 
1160
    //interleave soft info 1 -> 2
 
1161
    for(int k=0;k<blocklength;k++) {
 
1162
      int ki = INTERLEAVER.INTER()[k];
 
1163
      //for(int i=0;i<FSMi.I();i++) {
 
1164
        //oprioro[k*FSMi.I()+i]=iposti[ki*FSMi.I()+i];
 
1165
      //}
 
1166
      memcpy(&(priori2[k*FSM2.I()]),&(posti1[ki*FSM1.I()]),FSM1.I()*sizeof(float));
 
1167
    } 
 
1168
 
 
1169
    // run SISO 2
 
1170
    siso_algorithm(FSM2.I(),FSM2.S(),FSM2.O(),
 
1171
           FSM2.NS(), FSM2.OS(), FSM2.PS(), FSM2.PI(),
 
1172
           blocklength,
 
1173
           ST20,ST2K,
 
1174
           true, false,
 
1175
           p2mymin,
 
1176
           &(priori2[0]),  &(prioro2[0]), &(posti2[0])
 
1177
    );
 
1178
 
 
1179
    //interleave soft info 2 --> 1
 
1180
    for(int k=0;k<blocklength;k++) {
 
1181
      int ki = INTERLEAVER.INTER()[k];
 
1182
      //for(int i=0;i<FSMi.I();i++) {
 
1183
        //ipriori[ki*FSMi.I()+i]=oposto[k*FSMi.I()+i];
 
1184
      //}
 
1185
      memcpy(&(priori1[ki*FSM1.I()]),&(posti2[k*FSM2.I()]),FSM1.I()*sizeof(float));
 
1186
    }
 
1187
 
 
1188
  } // end iterations
 
1189
 
 
1190
  // generate hard decisions
 
1191
  for(int k=0;k<blocklength;k++) {
 
1192
    for(int i=0;i<FSM1.I();i++)
 
1193
      posti1[k*FSM1.I()+i]  = (*p2mymin)(priori1[k*FSM1.I()+i],posti1[k*FSM1.I()+i]);
 
1194
    float min=INF;
 
1195
    int mini=0;
 
1196
    for(int i=0;i<FSM1.I();i++) {
 
1197
      if(posti1[k*FSM1.I()+i]<min) {
 
1198
        min=posti1[k*FSM1.I()+i];
 
1199
        mini=i;
 
1200
      }
 
1201
    }
 
1202
    data[k]=(T)mini;
 
1203
    //std::cout << data[k] << ", "<< std::endl;
 
1204
  }
 
1205
  //std::cout << std::endl;
 
1206
 
 
1207
}
 
1208
 
 
1209
//----------------
 
1210
 
 
1211
template
 
1212
void pccc_decoder<unsigned char>(
 
1213
      const fsm &FSM1, int ST10, int ST1K,
 
1214
      const fsm &FSM2, int ST20, int ST2K,
 
1215
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1216
      float (*p2mymin)(float,float),
 
1217
      const float *cprioro, unsigned char *data
 
1218
);
 
1219
 
 
1220
template
 
1221
void pccc_decoder<short>(
 
1222
      const fsm &FSM1, int ST10, int ST1K,
 
1223
      const fsm &FSM2, int ST20, int ST2K,
 
1224
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1225
      float (*p2mymin)(float,float),
 
1226
      const float *cprioro, short *data
 
1227
);
 
1228
 
 
1229
template
 
1230
void pccc_decoder<int>(
 
1231
      const fsm &FSM1, int ST10, int ST1K,
 
1232
      const fsm &FSM2, int ST20, int ST2K,
 
1233
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1234
      float (*p2mymin)(float,float),
 
1235
      const float *cprioro, int *data
 
1236
);
 
1237
 
 
1238
 
 
1239
 
 
1240
//----------------
 
1241
 
 
1242
 
 
1243
template<class Ti, class To>
 
1244
void pccc_decoder_combined(
 
1245
      const fsm &FSM1, int ST10, int ST1K,
 
1246
      const fsm &FSM2, int ST20, int ST2K,
 
1247
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1248
      float (*p2mymin)(float,float),
 
1249
      int D, const std::vector<Ti> &TABLE,
 
1250
      trellis_metric_type_t METRIC_TYPE,
 
1251
      float scaling,
 
1252
      const Ti *observations, To *data
 
1253
)
 
1254
{
 
1255
 
 
1256
  //allocate space for cprioro
 
1257
  std::vector<float> cprioro(blocklength*FSM1.O()*FSM2.O(),0.0);
 
1258
 
 
1259
  //allocate space for priori, prioro and posti of FSM1
 
1260
  std::vector<float> priori1(blocklength*FSM1.I(),0.0);
 
1261
  std::vector<float> prioro1(blocklength*FSM1.O());
 
1262
  std::vector<float> posti1(blocklength*FSM1.I());
 
1263
 
 
1264
  //allocate space for priori, prioro and posti of FSM2
 
1265
  std::vector<float> priori2(blocklength*FSM2.I(),0.0);
 
1266
  std::vector<float> prioro2(blocklength*FSM2.O());
 
1267
  std::vector<float> posti2(blocklength*FSM2.I());
 
1268
 
 
1269
  // turn observations to neg-log-priors for cprioiro
 
1270
  int O=FSM1.O()*FSM2.O();
 
1271
  for(int k=0;k<blocklength;k++) {
 
1272
    calc_metric(O, D, TABLE, &(observations[k*D]), &(cprioro[k*O]),METRIC_TYPE);
 
1273
    cprioro[k*O] *= scaling;
 
1274
  }
 
1275
 
 
1276
  //generate prioro1,2 (metrics are not updated per iteration: this is not the best you can do...)
 
1277
  for (int k=0;k<blocklength;k++) {
 
1278
    //std::cout << k << std::endl;
 
1279
    for(int i=0;i<FSM1.O();i++) {
 
1280
      float x=cprioro[k*FSM1.O()*FSM2.O()+i*FSM1.O()+0];
 
1281
      for(int j=1;j<FSM2.O();j++)
 
1282
        x = (*p2mymin)(x,cprioro[k*FSM1.O()*FSM2.O()+i*FSM1.O()+j]);
 
1283
      prioro1[k*FSM1.O()+i]=x;
 
1284
      //std::cout <<  prioro1[k*FSM1.O()+i] << ", ";
 
1285
    }
 
1286
    //std::cout << std::endl;
 
1287
    for(int i=0;i<FSM2.O();i++) {
 
1288
      float x=cprioro[k*FSM1.O()*FSM2.O()+0*FSM1.O()+i];
 
1289
      for(int j=1;j<FSM1.O();j++)
 
1290
        x = (*p2mymin)(x,cprioro[k*FSM1.O()*FSM2.O()+j*FSM1.O()+i]);
 
1291
      prioro2[k*FSM2.O()+i]=x;
 
1292
    }
 
1293
  }
 
1294
 
 
1295
  for(int rep=0;rep<iterations;rep++) {
 
1296
    // run  SISO 1
 
1297
    siso_algorithm(FSM1.I(),FSM1.S(),FSM1.O(),
 
1298
             FSM1.NS(), FSM1.OS(), FSM1.PS(), FSM1.PI(),
 
1299
             blocklength,
 
1300
             ST10,ST1K,
 
1301
             true, false,
 
1302
             p2mymin,
 
1303
             &(priori1[0]),  &(prioro1[0]), &(posti1[0])
 
1304
    );
 
1305
   
 
1306
    //for(int k=0;k<blocklength;k++){
 
1307
      //for(int i=0;i<FSM1.I();i++)
 
1308
        //std::cout << posti1[k*FSM1.I()+i] << ", ";
 
1309
      //std::cout << std::endl;
 
1310
    //}
 
1311
 
 
1312
    //interleave soft info 1 -> 2
 
1313
    for(int k=0;k<blocklength;k++) {
 
1314
      int ki = INTERLEAVER.INTER()[k];
 
1315
      //for(int i=0;i<FSMi.I();i++) {
 
1316
        //oprioro[k*FSMi.I()+i]=iposti[ki*FSMi.I()+i];
 
1317
      //}
 
1318
      memcpy(&(priori2[k*FSM2.I()]),&(posti1[ki*FSM1.I()]),FSM1.I()*sizeof(float));
 
1319
    } 
 
1320
 
 
1321
    // run SISO 2
 
1322
    siso_algorithm(FSM2.I(),FSM2.S(),FSM2.O(),
 
1323
           FSM2.NS(), FSM2.OS(), FSM2.PS(), FSM2.PI(),
 
1324
           blocklength,
 
1325
           ST20,ST2K,
 
1326
           true, false,
 
1327
           p2mymin,
 
1328
           &(priori2[0]),  &(prioro2[0]), &(posti2[0])
 
1329
    );
 
1330
 
 
1331
    //interleave soft info 2 --> 1
 
1332
    for(int k=0;k<blocklength;k++) {
 
1333
      int ki = INTERLEAVER.INTER()[k];
 
1334
      //for(int i=0;i<FSMi.I();i++) {
 
1335
        //ipriori[ki*FSMi.I()+i]=oposto[k*FSMi.I()+i];
 
1336
      //}
 
1337
      memcpy(&(priori1[ki*FSM1.I()]),&(posti2[k*FSM2.I()]),FSM1.I()*sizeof(float));
 
1338
    }
 
1339
 
 
1340
  } // end iterations
 
1341
 
 
1342
  // generate hard decisions
 
1343
  for(int k=0;k<blocklength;k++) {
 
1344
    for(int i=0;i<FSM1.I();i++)
 
1345
      posti1[k*FSM1.I()+i]  = (*p2mymin)(priori1[k*FSM1.I()+i],posti1[k*FSM1.I()+i]);
 
1346
    float min=INF;
 
1347
    int mini=0;
 
1348
    for(int i=0;i<FSM1.I();i++) {
 
1349
      if(posti1[k*FSM1.I()+i]<min) {
 
1350
        min=posti1[k*FSM1.I()+i];
 
1351
        mini=i;
 
1352
      }
 
1353
    }
 
1354
    data[k]=(To)mini;
 
1355
    //std::cout << data[k] << ", "<< std::endl;
 
1356
  }
 
1357
  //std::cout << std::endl;
 
1358
 
 
1359
}
 
1360
 
 
1361
 
 
1362
template
 
1363
void pccc_decoder_combined(
 
1364
      const fsm &FSM1, int ST10, int ST1K,
 
1365
      const fsm &FSM2, int ST20, int ST2K,
 
1366
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1367
      float (*p2mymin)(float,float),
 
1368
      int D, const std::vector<float> &TABLE,
 
1369
      trellis_metric_type_t METRIC_TYPE,
 
1370
      float scaling,
 
1371
      const float *observations, unsigned char *data
 
1372
);
 
1373
 
 
1374
 
 
1375
template
 
1376
void pccc_decoder_combined(
 
1377
      const fsm &FSM1, int ST10, int ST1K,
 
1378
      const fsm &FSM2, int ST20, int ST2K,
 
1379
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1380
      float (*p2mymin)(float,float),
 
1381
      int D, const std::vector<float> &TABLE,
 
1382
      trellis_metric_type_t METRIC_TYPE,
 
1383
      float scaling,
 
1384
      const float *observations, short *data
 
1385
);
 
1386
 
 
1387
 
 
1388
template
 
1389
void pccc_decoder_combined(
 
1390
      const fsm &FSM1, int ST10, int ST1K,
 
1391
      const fsm &FSM2, int ST20, int ST2K,
 
1392
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1393
      float (*p2mymin)(float,float),
 
1394
      int D, const std::vector<float> &TABLE,
 
1395
      trellis_metric_type_t METRIC_TYPE,
 
1396
      float scaling,
 
1397
      const float *observations, int *data
 
1398
);
 
1399
 
 
1400
 
 
1401
template
 
1402
void pccc_decoder_combined(
 
1403
      const fsm &FSM1, int ST10, int ST1K,
 
1404
      const fsm &FSM2, int ST20, int ST2K,
 
1405
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1406
      float (*p2mymin)(float,float),
 
1407
      int D, const std::vector<gr_complex> &TABLE,
 
1408
      trellis_metric_type_t METRIC_TYPE,
 
1409
      float scaling,
 
1410
      const gr_complex *observations, unsigned char *data
 
1411
);
 
1412
 
 
1413
 
 
1414
template
 
1415
void pccc_decoder_combined(
 
1416
      const fsm &FSM1, int ST10, int ST1K,
 
1417
      const fsm &FSM2, int ST20, int ST2K,
 
1418
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1419
      float (*p2mymin)(float,float),
 
1420
      int D, const std::vector<gr_complex> &TABLE,
 
1421
      trellis_metric_type_t METRIC_TYPE,
 
1422
      float scaling,
 
1423
      const gr_complex *observations, short *data
 
1424
);
 
1425
 
 
1426
 
 
1427
template
 
1428
void pccc_decoder_combined(
 
1429
      const fsm &FSM1, int ST10, int ST1K,
 
1430
      const fsm &FSM2, int ST20, int ST2K,
 
1431
      const interleaver &INTERLEAVER, int blocklength, int iterations,
 
1432
      float (*p2mymin)(float,float),
 
1433
      int D, const std::vector<gr_complex> &TABLE,
 
1434
      trellis_metric_type_t METRIC_TYPE,
 
1435
      float scaling,
 
1436
      const gr_complex *observations, int *data
 
1437
);