src/conv.c

   1 /*
   2  * conv.c
   3  *
   4  * Generic convolutional encoding / decoding
   5  *
   6  * Copyright (C) 2011  Sylvain Munaut <tnt@246tNt.com>
   7  *
   8  * All Rights Reserved
   9  *
  10  * This program is free software; you can redistribute it and/or modify
  11  * it under the terms of the GNU General Public License as published by
  12  * the Free Software Foundation; either version 2 of the License, or
  13  * (at your option) any later version.
  14  *
  15  * This program is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  * GNU General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU General Public License along
  21  * with this program; if not, write to the Free Software Foundation, Inc.,
  22  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  23  */
  24
  25 #include <alloca.h>
  26 #include <stdint.h>
  27 #include <stdlib.h>
  28 #include <string.h>
  29
  30 #include <osmocom/core/bits.h>
  31 #include <osmocom/core/conv.h>
  32
  33
  34 /* ------------------------------------------------------------------------ */
  35 /* Encoding                                                                 */
  36 /* ------------------------------------------------------------------------ */
  37
  38 void
  39 osmo_conv_encode_init(struct osmo_conv_encoder *encoder,
  40                       const struct osmo_conv_code *code)
  41 {
  42         memset(encoder, 0x00, sizeof(struct osmo_conv_encoder));
  43         encoder->code = code;
  44 }
  45
  46 static inline int
  47 _conv_encode_do_output(struct osmo_conv_encoder *encoder,
  48                        uint8_t out, ubit_t *output)
  49 {
  50         const struct osmo_conv_code *code = encoder->code;
  51         int o_idx = 0;
  52         int j;
  53
  54         if (code->puncture) {
  55                 for (j=0; j<code->N; j++)
  56                 {
  57                         int bit_no = code->N - j - 1;
  58                         int r_idx = encoder->i_idx * code->N + j;
  59
  60                         if (code->puncture[encoder->p_idx] == r_idx)
  61                                 encoder->p_idx++;
  62                         else
  63                                 output[o_idx++] = (out >> bit_no) & 1;
  64                 }
  65         } else {
  66                 for (j=0; j<code->N; j++)
  67                 {
  68                         int bit_no = code->N - j - 1;
  69                         output[o_idx++] = (out >> bit_no) & 1;
  70                 }
  71         }
  72
  73         return o_idx;
  74 }
  75
  76 int
  77 osmo_conv_encode_raw(struct osmo_conv_encoder *encoder,
  78                      const ubit_t *input, ubit_t *output, int n)
  79 {
  80         const struct osmo_conv_code *code = encoder->code;
  81         uint8_t state;
  82         int i;
  83         int o_idx;
  84
  85         o_idx = 0;
  86         state = encoder->state;
  87
  88         for (i=0; i<n; i++) {
  89                 int bit = input[i];
  90                 uint8_t out;
  91
  92                 out   = code->next_output[state][bit];
  93                 state = code->next_state[state][bit];
  94
  95                 o_idx += _conv_encode_do_output(encoder, out, &output[o_idx]);
  96
  97                 encoder->i_idx++;
  98         }
  99
 100         encoder->state = state;
 101
 102         return o_idx;
 103 }
 104
 105 int
 106 osmo_conv_encode_finish(struct osmo_conv_encoder *encoder,
 107                         ubit_t *output)
 108 {
 109         const struct osmo_conv_code *code = encoder->code;
 110         uint8_t state;
 111         int n;
 112         int i;
 113         int o_idx;
 114
 115         n = code->K - 1;
 116
 117         o_idx = 0;
 118         state = encoder->state;
 119
 120         for (i=0; i<n; i++) {
 121                 uint8_t out;
 122
 123                 if (code->next_term_output) {
 124                         out   = code->next_term_output[state];
 125                         state = code->next_term_state[state];
 126                 } else {
 127                         out   = code->next_output[state][0];
 128                         state = code->next_state[state][0];
 129                 }
 130
 131                 o_idx += _conv_encode_do_output(encoder, out, &output[o_idx]);
 132
 133                 encoder->i_idx++;
 134         }
 135
 136         encoder->state = state;
 137
 138         return o_idx;
 139 }
 140
 141 int
 142 osmo_conv_encode(const struct osmo_conv_code *code,
 143                  const ubit_t *input, ubit_t *output)
 144 {
 145         struct osmo_conv_encoder encoder;
 146         int l;
 147
 148         osmo_conv_encode_init(&encoder, code);
 149         l  = osmo_conv_encode_raw(&encoder, input, output, code->len);
 150         l += osmo_conv_encode_finish(&encoder, &output[l]);
 151
 152         return l;
 153 }
 154
 155
 156 /* ------------------------------------------------------------------------ */
 157 /* Decoding (viterbi)                                                       */
 158 /* ------------------------------------------------------------------------ */
 159
 160 #define MAX_AE 0x00ffffff
 161
 162 void
 163 osmo_conv_decode_init(struct osmo_conv_decoder *decoder,
 164                       const struct osmo_conv_code *code, int len)
 165 {
 166         int n_states;
 167
 168         /* Init */
 169         if (len <= 0)
 170                 len =  code->len;
 171
 172         n_states = 1 << (code->K - 1);
 173
 174         memset(decoder, 0x00, sizeof(struct osmo_conv_decoder));
 175
 176         decoder->code = code;
 177         decoder->n_states = n_states;
 178         decoder->len = len;
 179
 180         /* Allocate arrays */
 181         decoder->ae      = malloc(sizeof(unsigned int) * n_states);
 182         decoder->ae_next = malloc(sizeof(unsigned int) * n_states);
 183
 184         decoder->state_history = malloc(sizeof(uint8_t) * n_states * (len + decoder->code->K - 1));
 185
 186         /* Classic reset */
 187         osmo_conv_decode_reset(decoder);
 188 }
 189
 190 void
 191 osmo_conv_decode_reset(struct osmo_conv_decoder *decoder)
 192 {
 193         int i;
 194
 195         /* Reset indexes */
 196         decoder->o_idx = 0;
 197         decoder->p_idx = 0;
 198
 199         /* Initial error (only state 0 is valid) */
 200         decoder->ae[0] = 0;
 201         for (i=1; i<decoder->n_states; i++) {
 202                 decoder->ae[i] = MAX_AE;
 203         }
 204 }
 205
 206 void
 207 osmo_conv_decode_deinit(struct osmo_conv_decoder *decoder)
 208 {
 209         free(decoder->ae);
 210         free(decoder->ae_next);
 211         free(decoder->state_history);
 212
 213         memset(decoder, 0x00, sizeof(struct osmo_conv_decoder));
 214 }
 215
 216 int
 217 osmo_conv_decode_scan(struct osmo_conv_decoder *decoder,
 218                       const sbit_t *input, int n)
 219 {
 220         const struct osmo_conv_code *code = decoder->code;
 221
 222         int i, s, b, j;
 223
 224         int n_states;
 225         unsigned int *ae;
 226         unsigned int *ae_next;
 227         uint8_t *state_history;
 228         sbit_t *in_sym;
 229
 230         int i_idx, p_idx;
 231
 232         /* Prepare */
 233         n_states = decoder->n_states;
 234
 235         ae      = decoder->ae;
 236         ae_next = decoder->ae_next;
 237         state_history = &decoder->state_history[n_states * decoder->o_idx];
 238
 239         in_sym  = alloca(sizeof(sbit_t) * code->N);
 240
 241         i_idx = 0;
 242         p_idx = decoder->p_idx;
 243
 244         /* Scan the treillis */
 245         for (i=0; i<n; i++)
 246         {
 247                 /* Reset next accumulated error */
 248                 for (s=0; s<n_states; s++) {
 249                         ae_next[s] = MAX_AE;
 250                 }
 251
 252                 /* Get input */
 253                 if (code->puncture) {
 254                         /* Hard way ... */
 255                         for (j=0; j<code->N; j++) {
 256                                 int idx = ((decoder->o_idx + i) * code->N) + j;
 257                                 if (idx == code->puncture[p_idx]) {
 258                                         in_sym[j] = 0;  /* Undefined */
 259                                         p_idx++;
 260                                 } else {
 261                                         in_sym[j] = input[i_idx];
 262                                         i_idx++;
 263                                 }
 264                         }
 265                 } else {
 266                         /* Easy, just copy N bits */
 267                         memcpy(in_sym, &input[i_idx], code->N);
 268                         i_idx += code->N;
 269                 }
 270
 271                 /* Scan all state */
 272                 for (s=0; s<n_states; s++)
 273                 {
 274                         /* Scan possible input bits */
 275                         for (b=0; b<2; b++)
 276                         {
 277                                 int nae, ov, e;
 278                                 uint8_t m;
 279
 280                                 /* Next output and state */
 281                                 uint8_t out   = code->next_output[s][b];
 282                                 uint8_t state = code->next_state[s][b];
 283
 284                                 /* New error for this path */
 285                                 nae = ae[s];                    /* start from last error */
 286                                 m = 1 << (code->N - 1);         /* mask for 'out' bit selection */
 287
 288                                 for (j=0; j<code->N; j++) {
 289                                         ov = (out & m) ? -127 : 127; /* sbit_t value for it */
 290                                         e = ((int)in_sym[j]) - ov;   /* raw error for this bit */
 291                                         nae += (e * e) >> 9;         /* acc the squared/scaled value */
 292                                         m >>= 1;                     /* next mask bit */
 293                                 }
 294
 295                                 /* Is it survivor ? */
 296                                 if (ae_next[state] > nae) {
 297                                         ae_next[state] = nae;
 298                                         state_history[(n_states * i) + state] = s;
 299                                 }
 300                         }
 301                 }
 302
 303                 /* Copy accumulated error */
 304                 memcpy(ae, ae_next, sizeof(unsigned int) * n_states);
 305         }
 306
 307         /* Update decoder state */
 308         decoder->p_idx = p_idx;
 309         decoder->o_idx += n;
 310
 311         return i_idx;
 312 }
 313
 314 int
 315 osmo_conv_decode_finish(struct osmo_conv_decoder *decoder,
 316                         const sbit_t *input)
 317 {
 318         const struct osmo_conv_code *code = decoder->code;
 319
 320         int i, s, j;
 321
 322         int n_states;
 323         unsigned int *ae;
 324         unsigned int *ae_next;
 325         uint8_t *state_history;
 326         sbit_t *in_sym;
 327
 328         int i_idx, p_idx;
 329
 330         /* Prepare */
 331         n_states = decoder->n_states;
 332
 333         ae      = decoder->ae;
 334         ae_next = decoder->ae_next;
 335         state_history = &decoder->state_history[n_states * decoder->o_idx];
 336
 337         in_sym  = alloca(sizeof(sbit_t) * code->N);
 338
 339         i_idx = 0;
 340         p_idx = decoder->p_idx;
 341
 342         /* Scan the treillis */
 343         for (i=0; i<code->K-1; i++)
 344         {
 345                 /* Reset next accumulated error */
 346                 for (s=0; s<n_states; s++) {
 347                         ae_next[s] = MAX_AE;
 348                 }
 349
 350                 /* Get input */
 351                 if (code->puncture) {
 352                         /* Hard way ... */
 353                         for (j=0; j<code->N; j++) {
 354                                 int idx = ((decoder->o_idx + i) * code->N) + j;
 355                                 if (idx == code->puncture[p_idx]) {
 356                                         in_sym[j] = 0;  /* Undefined */
 357                                         p_idx++;
 358                                 } else {
 359                                         in_sym[j] = input[i_idx];
 360                                         i_idx++;
 361                                 }
 362                         }
 363                 } else {
 364                         /* Easy, just copy N bits */
 365                         memcpy(in_sym, &input[i_idx], code->N);
 366                         i_idx += code->N;
 367                 }
 368
 369                 /* Scan all state */
 370                 for (s=0; s<n_states; s++)
 371                 {
 372                         int nae, ov, e;
 373                         uint8_t m;
 374
 375                         /* Next output and state */
 376                         uint8_t out;
 377                         uint8_t state;
 378
 379                         if (code->next_term_output) {
 380                                 out   = code->next_term_output[s];
 381                                 state = code->next_term_state[s];
 382                         } else {
 383                                 out   = code->next_output[s][0];
 384                                 state = code->next_state[s][0];
 385                         }
 386
 387                         /* New error for this path */
 388                         nae = ae[s];                    /* start from last error */
 389                         m = 1 << (code->N - 1);         /* mask for 'out' bit selection */
 390
 391                         for (j=0; j<code->N; j++) {
 392                                 int is = (int)in_sym[j];
 393                                 if (is) {
 394                                         ov = (out & m) ? -127 : 127; /* sbit_t value for it */
 395                                         e = is - ov;                 /* raw error for this bit */
 396                                         nae += (e * e) >> 9;         /* acc the squared/scaled value */
 397                                 }
 398                                 m >>= 1;                     /* next mask bit */
 399                         }
 400
 401                         /* Is it survivor ? */
 402                         if (ae_next[state] > nae) {
 403                                 ae_next[state] = nae;
 404                                 state_history[(n_states * i) + state] = s;
 405                         }
 406                 }
 407
 408                 /* Copy accumulated error */
 409                 memcpy(ae, ae_next, sizeof(unsigned int) * n_states);
 410         }
 411
 412         /* Update decoder state */
 413         decoder->p_idx = p_idx;
 414         decoder->o_idx += code->K - 1;
 415
 416         return i_idx;
 417 }
 418
 419 int
 420 osmo_conv_decode_get_output(struct osmo_conv_decoder *decoder,
 421                             ubit_t *output, int has_finish)
 422 {
 423         const struct osmo_conv_code *code = decoder->code;
 424
 425         int min_ae;
 426         uint8_t min_state, cur_state;
 427         int i, s, n;
 428
 429         uint8_t *sh_ptr;
 430
 431         /* Find state with least error */
 432         min_ae = MAX_AE;
 433         min_state = 0xff;
 434
 435         for (s=0; s<decoder->n_states; s++)
 436         {
 437                 if (decoder->ae[s] < min_ae) {
 438                         min_ae = decoder->ae[s];
 439                         min_state = s;
 440                 }
 441         }
 442
 443         if (min_state == 0xff)
 444                 return -1;
 445
 446         /* Traceback */
 447         cur_state = min_state;
 448
 449         n = decoder->o_idx;
 450
 451         sh_ptr = &decoder->state_history[decoder->n_states * (n-1)];
 452
 453                 /* No output for the K-1 termination input bits */
 454         if (has_finish) {
 455                 for (i=0; i<code->K-1; i++) {
 456                         cur_state = sh_ptr[cur_state];
 457                         sh_ptr -= decoder->n_states;
 458                 }
 459                 n -= code->K - 1;
 460         }
 461
 462                 /* Generate output backward */
 463         for (i=n-1; i>=0; i--)
 464         {
 465                 min_state = cur_state;
 466                 cur_state = sh_ptr[cur_state];
 467
 468                 sh_ptr -= decoder->n_states;
 469
 470                 if (code->next_state[cur_state][0] == min_state)
 471                         output[i] = 0;
 472                 else
 473                         output[i] = 1;
 474         }
 475
 476         return min_ae;
 477 }
 478
 479 int
 480 osmo_conv_decode(const struct osmo_conv_code *code,
 481                  const sbit_t *input, ubit_t *output)
 482 {
 483         struct osmo_conv_decoder decoder;
 484         int rv, l;
 485
 486         osmo_conv_decode_init(&decoder, code, 0);
 487
 488         l = osmo_conv_decode_scan(&decoder, input, code->len);
 489         l = osmo_conv_decode_finish(&decoder, &input[l]);
 490
 491         rv = osmo_conv_decode_get_output(&decoder, output, 1);
 492
 493         osmo_conv_decode_deinit(&decoder);
 494
 495         return rv;
 496 }