git.gag.com Git - debian/gnuradio/blob - gr-trellis/src/lib/trellis_siso_f.cc

   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 2, 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 #ifdef HAVE_CONFIG_H
  24 #include "config.h"
  25 #endif
  26
  27 #include <trellis_siso_f.h>
  28 #include <gr_io_signature.h>
  29 #include <stdexcept>
  30 #include <assert.h>
  31 #include <iostream>
  32
  33 static const float INF = 1.0e9;
  34
  35 trellis_siso_f_sptr
  36 trellis_make_siso_f (
  37     const fsm &FSM,
  38     int K,
  39     int S0,
  40     int SK,
  41     bool POSTI,
  42     bool POSTO,
  43     trellis_siso_type_t SISO_TYPE)
  44 {
  45   return trellis_siso_f_sptr (new trellis_siso_f (FSM,K,S0,SK,POSTI,POSTO,SISO_TYPE));
  46 }
  47
  48 trellis_siso_f::trellis_siso_f (
  49     const fsm &FSM,
  50     int K,
  51     int S0,
  52     int SK,
  53     bool POSTI,
  54     bool POSTO,
  55     trellis_siso_type_t SISO_TYPE)
  56   : gr_block ("siso_f",
  57                           gr_make_io_signature (1, -1, sizeof (float)),
  58                           gr_make_io_signature (1, -1, sizeof (float))),
  59   d_FSM (FSM),
  60   d_K (K),
  61   d_S0 (S0),
  62   d_SK (SK),
  63   d_POSTI (POSTI),
  64   d_POSTO (POSTO),
  65   d_SISO_TYPE (SISO_TYPE)//,
  66   //d_alpha(FSM.S()*(K+1)),
  67   //d_beta(FSM.S()*(K+1))
  68 {
  69     int multiple;
  70     if (d_POSTI && d_POSTO)
  71         multiple = d_FSM.I()+d_FSM.O();
  72     else if(d_POSTI)
  73         multiple = d_FSM.I();
  74     else if(d_POSTO)
  75         multiple = d_FSM.O();
  76     else
  77         throw std::runtime_error ("Not both POSTI and POSTO can be false.");
  78     //printf("constructor: Multiple = %d\n",multiple);
  79     set_output_multiple (d_K*multiple);
  80     //what is the meaning of relative rate for a block with 2 inputs?
  81     //set_relative_rate ( multiple / ((double) d_FSM.I()) );
  82     // it turns out that the above gives problems in the scheduler, so
  83     // let's try (assumption O>I)
  84     //set_relative_rate ( multiple / ((double) d_FSM.O()) );
  85     // I am tempted to automate like this
  86     if(d_FSM.I() <= d_FSM.O())
  87       set_relative_rate ( multiple / ((double) d_FSM.O()) );
  88     else
  89       set_relative_rate ( multiple / ((double) d_FSM.I()) );
  90 }
  91
  92
  93 void
  94 trellis_siso_f::forecast (int noutput_items, gr_vector_int &ninput_items_required)
  95 {
  96   int multiple;
  97   if (d_POSTI && d_POSTO)
  98       multiple = d_FSM.I()+d_FSM.O();
  99   else if(d_POSTI)
 100       multiple = d_FSM.I();
 101   else if(d_POSTO)
 102       multiple = d_FSM.O();
 103   else
 104       throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 105   //printf("forecast: Multiple = %d\n",multiple);
 106   assert (noutput_items % (d_K*multiple) == 0);
 107   int input_required1 =  d_FSM.I() * (noutput_items/multiple) ;
 108   int input_required2 =  d_FSM.O() * (noutput_items/multiple) ;
 109   //printf("forecast: Output requirements:  %d\n",noutput_items);
 110   //printf("forecast: Input requirements:  %d   %d\n",input_required1,input_required2);
 111   unsigned ninputs = ninput_items_required.size();
 112   assert(ninputs % 2 == 0);
 113   for (unsigned int i = 0; i < ninputs/2; i++) {
 114     ninput_items_required[2*i] = input_required1;
 115     ninput_items_required[2*i+1] = input_required2;
 116   }
 117 }
 118
 119 inline float min(float a, float b)
 120 {
 121   return a <= b ? a : b;
 122 }
 123
 124 inline float min_star(float a, float b)
 125 {
 126   return (a <= b ? a : b)-log(1+exp(a <= b ? a-b : b-a));
 127 }
 128
 129 void siso_algorithm(int I, int S, int O,
 130              const std::vector<int> &NS,
 131              const std::vector<int> &OS,
 132              const std::vector<int> &PS,
 133              const std::vector<int> &PI,
 134              int K,
 135              int S0,int SK,
 136              bool POSTI, bool POSTO,
 137              float (*p2mymin)(float,float),
 138              const float *priori, const float *prioro, float *post//,
 139              //std::vector<float> &alpha,
 140              //std::vector<float> &beta
 141              )
 142 {
 143   float norm,mm,minm;
 144   std::vector<float> alpha(S*(K+1));
 145   std::vector<float> beta(S*(K+1));
 146
 147
 148   if(S0<0) { // initial state not specified
 149       for(int i=0;i<S;i++) alpha[0*S+i]=0;
 150   }
 151   else {
 152       for(int i=0;i<S;i++) alpha[0*S+i]=INF;
 153       alpha[0*S+S0]=0.0;
 154   }
 155
 156   for(int k=0;k<K;k++) { // forward recursion
 157       norm=INF;
 158       for(int j=0;j<S;j++) {
 159           minm=INF;
 160           for(int i=0;i<I;i++) {
 161               int i0 = j*I+i;
 162               mm=alpha[k*S+PS[i0]]+priori[k*I+PI[i0]]+prioro[k*O+OS[PS[i0]*I+PI[i0]]];
 163               minm=(*p2mymin)(minm,mm);
 164           }
 165           alpha[(k+1)*S+j]=minm;
 166           if(minm<norm) norm=minm;
 167       }
 168       for(int j=0;j<S;j++)
 169           alpha[(k+1)*S+j]-=norm; // normalize total metrics so they do not explode
 170   }
 171
 172   if(SK<0) { // final state not specified
 173       for(int i=0;i<S;i++) beta[K*S+i]=0;
 174   }
 175   else {
 176       for(int i=0;i<S;i++) beta[K*S+i]=INF;
 177       beta[K*S+SK]=0.0;
 178   }
 179
 180   for(int k=K-1;k>=0;k--) { // backward recursion
 181       norm=INF;
 182       for(int j=0;j<S;j++) {
 183           minm=INF;
 184           for(int i=0;i<I;i++) {
 185               int i0 = j*I+i;
 186               mm=beta[(k+1)*S+NS[i0]]+priori[k*I+i]+prioro[k*O+OS[i0]];
 187               minm=(*p2mymin)(minm,mm);
 188           }
 189           beta[k*S+j]=minm;
 190           if(minm<norm) norm=minm;
 191       }
 192       for(int j=0;j<S;j++)
 193           beta[k*S+j]-=norm; // normalize total metrics so they do not explode
 194   }
 195
 196
 197 if (POSTI && POSTO)
 198 {
 199   for(int k=0;k<K;k++) { // input combining
 200       norm=INF;
 201       for(int i=0;i<I;i++) {
 202           minm=INF;
 203           for(int j=0;j<S;j++) {
 204               mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 205               minm=(*p2mymin)(minm,mm);
 206           }
 207           post[k*(I+O)+i]=minm;
 208           if(minm<norm) norm=minm;
 209       }
 210       for(int i=0;i<I;i++)
 211           post[k*(I+O)+i]-=norm; // normalize metrics
 212   }
 213
 214
 215   for(int k=0;k<K;k++) { // output combining
 216       norm=INF;
 217       for(int n=0;n<O;n++) {
 218           minm=INF;
 219           for(int j=0;j<S;j++) {
 220               for(int i=0;i<I;i++) {
 221                   mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 222                   minm=(*p2mymin)(minm,mm);
 223               }
 224           }
 225           post[k*(I+O)+I+n]=minm;
 226           if(minm<norm) norm=minm;
 227       }
 228       for(int n=0;n<O;n++)
 229           post[k*(I+O)+I+n]-=norm; // normalize metrics
 230   }
 231 }
 232 else if(POSTI)
 233 {
 234   for(int k=0;k<K;k++) { // input combining
 235       norm=INF;
 236       for(int i=0;i<I;i++) {
 237           minm=INF;
 238           for(int j=0;j<S;j++) {
 239               mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 240               minm=(*p2mymin)(minm,mm);
 241           }
 242           post[k*I+i]=minm;
 243           if(minm<norm) norm=minm;
 244       }
 245       for(int i=0;i<I;i++)
 246           post[k*I+i]-=norm; // normalize metrics
 247   }
 248 }
 249 else if(POSTO)
 250 {
 251   for(int k=0;k<K;k++) { // output combining
 252       norm=INF;
 253       for(int n=0;n<O;n++) {
 254           minm=INF;
 255           for(int j=0;j<S;j++) {
 256               for(int i=0;i<I;i++) {
 257                   mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 258                   minm=(*p2mymin)(minm,mm);
 259               }
 260           }
 261           post[k*O+n]=minm;
 262           if(minm<norm) norm=minm;
 263       }
 264       for(int n=0;n<O;n++)
 265           post[k*O+n]-=norm; // normalize metrics
 266   }
 267 }
 268 else
 269     throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 270
 271 }
 272
 273
 274
 275
 276
 277
 278 int
 279 trellis_siso_f::general_work (int noutput_items,
 280                         gr_vector_int &ninput_items,
 281                         gr_vector_const_void_star &input_items,
 282                         gr_vector_void_star &output_items)
 283 {
 284   assert (input_items.size() == 2*output_items.size());
 285   int nstreams = output_items.size();
 286   //printf("general_work:Streams:  %d\n",nstreams);
 287   int multiple;
 288   if (d_POSTI && d_POSTO)
 289       multiple = d_FSM.I()+d_FSM.O();
 290   else if(d_POSTI)
 291       multiple = d_FSM.I();
 292   else if(d_POSTO)
 293       multiple = d_FSM.O();
 294   else
 295       throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 296
 297   assert (noutput_items % (d_K*multiple) == 0);
 298   int nblocks = noutput_items / (d_K*multiple);
 299   //printf("general_work:Blocks:  %d\n",nblocks);
 300   //for(int i=0;i<ninput_items.size();i++)
 301       //printf("general_work:Input items available:  %d\n",ninput_items[i]);
 302
 303   float (*p2min)(float, float) = NULL;
 304   if(d_SISO_TYPE == TRELLIS_MIN_SUM)
 305     p2min = &min;
 306   else if(d_SISO_TYPE == TRELLIS_SUM_PRODUCT)
 307     p2min = &min_star;
 308
 309
 310   for (int m=0;m<nstreams;m++) {
 311     const float *in1 = (const float *) input_items[2*m];
 312     const float *in2 = (const float *) input_items[2*m+1];
 313     float *out = (float *) output_items[m];
 314     for (int n=0;n<nblocks;n++) {
 315       siso_algorithm(d_FSM.I(),d_FSM.S(),d_FSM.O(),
 316         d_FSM.NS(),d_FSM.OS(),d_FSM.PS(),d_FSM.PI(),
 317         d_K,d_S0,d_SK,
 318         d_POSTI,d_POSTO,
 319         p2min,
 320         &(in1[n*d_K*d_FSM.I()]),&(in2[n*d_K*d_FSM.O()]),
 321         &(out[n*d_K*multiple])//,
 322         //d_alpha,d_beta
 323         );
 324     }
 325   }
 326
 327   for (unsigned int i = 0; i < input_items.size()/2; i++) {
 328     consume(2*i,d_FSM.I() * noutput_items / multiple );
 329     consume(2*i+1,d_FSM.O() * noutput_items / multiple );
 330   }
 331
 332   return noutput_items;
 333 }