viterbi_lte.c 14.5 KB
Newer Older
1
/*******************************************************************************
2
    OpenAirInterface
ghaddab's avatar
ghaddab committed
3
    Copyright(c) 1999 - 2014 Eurecom
4

ghaddab's avatar
ghaddab committed
5 6 7 8
    OpenAirInterface is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
9 10


ghaddab's avatar
ghaddab committed
11 12 13 14
    OpenAirInterface is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
15

ghaddab's avatar
ghaddab committed
16
    You should have received a copy of the GNU General Public License
17 18
    along with OpenAirInterface.The full GNU General Public License is
   included in this distribution in the file called "COPYING". If not,
ghaddab's avatar
ghaddab committed
19
   see <http://www.gnu.org/licenses/>.
20 21

  Contact Information
ghaddab's avatar
ghaddab committed
22 23 24
  OpenAirInterface Admin: openair_admin@eurecom.fr
  OpenAirInterface Tech : openair_tech@eurecom.fr
  OpenAirInterface Dev  : openair4g-devel@eurecom.fr
25

ghaddab's avatar
ghaddab committed
26
  Address      : Eurecom, Campus SophiaTech, 450 Route des Chappes, CS 50193 - 06904 Biot Sophia Antipolis cedex, FRANCE
27

ghaddab's avatar
ghaddab committed
28
 *******************************************************************************/
29 30 31 32
/* file: viterbit_lte.c
   purpose: SIMD optimized LTE Viterbi Decoder for rate 1/3 Tail-biting convolutional code.  Performs two iterations
            of code.  First pass does Viterbi with all initial partial metrics set to zero.  Second pass does Viterbi
            with initial partial metrics set to values from final state values after first pass. Max is selected at
33
      end to do trace-back.
34
   author: raymond.knopp@eurecom.fr
35
   date: 21.10.2009
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
*/

#ifdef USER_MODE
#include <stdio.h>
#endif

#ifndef TEST_DEBUG
#include "PHY/defs.h"
#include "PHY/extern.h"
#else
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define msg printf
#endif

52

53
#include "PHY/sse_intrin.h"
54

55
extern uint8_t ccodelte_table[128],ccodelte_table_rev[128];
56 57 58 59




60 61
static int8_t m0_table[64*16*16*16] __attribute__ ((aligned(16)));
static int8_t m1_table[64*16*16*16] __attribute__ ((aligned(16)));
62 63 64


// Set up Viterbi tables for SSE2 implementation
65 66
void phy_generate_viterbi_tables_lte()
{
67

68 69
  int8_t w[8],in0,in1,in2;
  uint8_t state,index0,index1;
70 71

  for (in0 = -8; in0 <8 ; in0++) {   // use 4-bit quantization
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
    for (in1 = -8; in1 <8 ; in1++) {
      for (in2 = -8; in2 <8 ; in2++) {


        w[0] = 24 - in0 - in1 - in2;           // -1,-1,-1
        w[1] = 24 + in0 - in1 - in2;           // -1, 1,-1
        w[2] = 24 - in0 + in1 - in2;           //  1,-1,-1
        w[3] = 24 + in0 + in1 - in2;           //  1, 1,-1
        w[4] = 24 - in0 - in1 + in2;           // -1,-1, 1
        w[5] = 24 + in0 - in1 + in2;           // -1, 1, 1
        w[6] = 24 - in0 + in1 + in2;           //  1,-1, 1
        w[7] = 24 + in0 + in1 + in2;           //  1, 1, 1

        //    printf("w: %d %d %d %d\n",w[0],w[1],w[2],w[3]);
        for (state=0; state<64 ; state++) {

          // input 0
          index0 = (state<<1);

          m0_table[(in0+8 + (16*(in1+8)) + (256*(in2+8)))*64 +state] = w[ccodelte_table_rev[index0]];


          //    if (position < 8)
          //    printf("%d,%d : prev_state0 = %d,m0 = %d,w=%d (%d)\n",position,state,prev_state0%64,m0,w[ccodelte_table[prev_state0]],partial_metrics[prev_state0%64]);

          // input 1
          index1 = (1+ (state<<1));
          m1_table[(in0+8 + (16*(in1+8)) + (256*(in2+8)))*64 +state] = w[ccodelte_table_rev[index1]];

        }
102 103 104 105 106 107 108 109 110 111 112 113
      }
    }
  }
}


#define INIT0 0x00000080
#define RESCALE 0x00000040

//#define DEBUG_VITERBI

#ifdef DEBUG_VITERBI
114 115
void print_bytes(char *s,__m128i *x)
{
116

117
  uint8_t *tempb = (uint8_t *)x;
118 119

  printf("%s  : %d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d\n",s,
120 121
         tempb[0],tempb[1],tempb[2],tempb[3],tempb[4],tempb[5],tempb[6],tempb[7],
         tempb[8],tempb[9],tempb[10],tempb[11],tempb[12],tempb[13],tempb[14],tempb[15]);
122 123 124 125 126 127

}

/*
void print_shorts(__m128i x,char *s) {

128
  int16_t *tempb = (int16_t *)&x;
129 130

  printf("%s  : %d,%d,%d,%d,%d,%d,%d,%d\n",s,
131 132
   tempb[0],tempb[1],tempb[2],tempb[3],tempb[4],tempb[5],tempb[6],tempb[7]
   );
133 134 135 136 137 138 139

}
*/
#endif // USER_MODE



140 141 142
void phy_viterbi_lte_sse2(int8_t *y,uint8_t *decoded_bytes,uint16_t n)
{

143

144 145 146 147 148 149 150 151
#if defined(__x86_64__) || defined(__i386__)
  __m128i  TB[4*8192];
  __m128i *m0_ptr,*m1_ptr,*TB_ptr = &TB[0];
  
  __m128i metrics0_15,metrics16_31,metrics32_47,metrics48_63,even0_30a,even0_30b,even32_62a,even32_62b,odd1_31a,odd1_31b,odd33_63a,odd33_63b,TBeven0_30,TBeven32_62,TBodd1_31,
    TBodd33_63;
  
  __m128i min_state,min_state2;
152

153 154
#elif defined(__arm__)
  uint8x16x2_t TB[2*8192];  // 2 int8x16_t per input bit, 8 bits / byte, 8192 is largest packet size in bits
155

156 157 158 159 160 161 162 163 164 165
  uint8x16_t even0_30a,even0_30b,even32_62a,even32_62b,odd1_31a,odd1_31b,odd33_63a,odd33_63b,TBeven0_30,TBeven32_62,TBodd1_31,TBodd33_63;
  uint8x16x2_t metrics0_31,metrics32_63;

  uint8x16_t min_state;

  uint8x16_t *m0_ptr,*m1_ptr;
  uint8x16x2_t *TB_ptr = &TB[0];


#endif
166 167 168 169 170 171
  int8_t *in = y;
  uint8_t prev_state0,maxm,s;
  static uint8_t *TB_ptr2;
  uint32_t table_offset;
  uint8_t iter;
  int16_t position;
172 173 174 175

  // set initial metrics
  //debug_msg("Doing viterbi\n");

176 177 178
#if defined(__x86_64__) || defined(__i386__)

  metrics0_15  = _mm_setzero_si128();
179 180 181
  metrics16_31 = _mm_setzero_si128();
  metrics32_47 = _mm_setzero_si128();
  metrics48_63 = _mm_setzero_si128();
182 183 184 185 186
#elif defined(__arm__)
    metrics0_31.val[0]  = vdupq_n_u8(0); 
    metrics0_31.val[1]  = vdupq_n_u8(0);
    metrics32_63.val[0] = vdupq_n_u8(0);
    metrics32_63.val[1] = vdupq_n_u8(0);
187 188
#endif

189
  for (iter=0; iter<2; iter++) {
190 191 192
    in = y;
    TB_ptr=&TB[0];

193 194 195
    for (position=0; position<n; position++) {


196 197
      // get branch metric offsets for the 64 states
      table_offset = (in[0]+8 + ((in[1]+8)<<4) + ((in[2]+8)<<8))<<6;
198 199

#if defined(__x86_64__) || defined(__i386__)
200 201
      m0_ptr = (__m128i *)&m0_table[table_offset];
      m1_ptr = (__m128i *)&m1_table[table_offset];
202

203 204 205 206 207
      // even states
      even0_30a  = _mm_adds_epu8(metrics0_15,m0_ptr[0]);
      even32_62a = _mm_adds_epu8(metrics16_31,m0_ptr[1]);
      even0_30b  = _mm_adds_epu8(metrics32_47,m0_ptr[2]);
      even32_62b = _mm_adds_epu8(metrics48_63,m0_ptr[3]);
208 209


210 211 212 213 214
      // odd states
      odd1_31a   = _mm_adds_epu8(metrics0_15,m1_ptr[0]);
      odd33_63a  = _mm_adds_epu8(metrics16_31,m1_ptr[1]);
      odd1_31b   = _mm_adds_epu8(metrics32_47,m1_ptr[2]);
      odd33_63b  = _mm_adds_epu8(metrics48_63,m1_ptr[3]);
215

216
      // select maxima
217

218 219 220 221
      even0_30a  = _mm_max_epu8(even0_30a,even0_30b);
      even32_62a = _mm_max_epu8(even32_62a,even32_62b);
      odd1_31a   = _mm_max_epu8(odd1_31a,odd1_31b);
      odd33_63a  = _mm_max_epu8(odd33_63a,odd33_63b);
222

223
      // Traceback information
224

225 226 227 228
      TBeven0_30  = _mm_cmpeq_epi8(even0_30a,even0_30b);
      TBeven32_62 = _mm_cmpeq_epi8(even32_62a,even32_62b);
      TBodd1_31   = _mm_cmpeq_epi8(odd1_31a,odd1_31b);
      TBodd33_63  = _mm_cmpeq_epi8(odd33_63a,odd33_63b);
229

230 231 232 233
      metrics0_15        = _mm_unpacklo_epi8(even0_30a ,odd1_31a);
      metrics16_31       = _mm_unpackhi_epi8(even0_30a ,odd1_31a);
      metrics32_47       = _mm_unpacklo_epi8(even32_62a,odd33_63a);
      metrics48_63       = _mm_unpackhi_epi8(even32_62a,odd33_63a);
234

235 236 237 238
      TB_ptr[0]  = _mm_unpacklo_epi8(TBeven0_30,TBodd1_31);
      TB_ptr[1] = _mm_unpackhi_epi8(TBeven0_30,TBodd1_31);
      TB_ptr[2] = _mm_unpacklo_epi8(TBeven32_62,TBodd33_63);
      TB_ptr[3] = _mm_unpackhi_epi8(TBeven32_62,TBodd33_63);
239 240


241 242
      in+=3;
      TB_ptr += 4;
243

244 245
      // rescale by subtracting minimum
      /****************************************************
246
      USE SSSE instruction phminpos!!!!!!!
247 248 249 250
      ****************************************************/
      min_state =_mm_min_epu8(metrics0_15,metrics16_31);
      min_state =_mm_min_epu8(min_state,metrics32_47);
      min_state =_mm_min_epu8(min_state,metrics48_63);
251

252 253 254 255
      min_state2 = min_state;
      min_state  = _mm_unpacklo_epi8(min_state,min_state);
      min_state2 = _mm_unpackhi_epi8(min_state2,min_state2);
      min_state  = _mm_min_epu8(min_state,min_state2);
256

257 258 259 260
      min_state2 = min_state;
      min_state  = _mm_unpacklo_epi8(min_state,min_state);
      min_state2 = _mm_unpackhi_epi8(min_state2,min_state2);
      min_state  = _mm_min_epu8(min_state,min_state2);
261

262 263 264 265
      min_state2 = min_state;
      min_state  = _mm_unpacklo_epi8(min_state,min_state);
      min_state2 = _mm_unpackhi_epi8(min_state2,min_state2);
      min_state  = _mm_min_epu8(min_state,min_state2);
266

267 268 269 270
      min_state2 = min_state;
      min_state  = _mm_unpacklo_epi8(min_state,min_state);
      min_state2 = _mm_unpackhi_epi8(min_state2,min_state2);
      min_state  = _mm_min_epu8(min_state,min_state2);
271

272 273 274 275
      metrics0_15  = _mm_subs_epu8(metrics0_15,min_state);
      metrics16_31 = _mm_subs_epu8(metrics16_31,min_state);
      metrics32_47 = _mm_subs_epu8(metrics32_47,min_state);
      metrics48_63 = _mm_subs_epu8(metrics48_63,min_state);
276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
#elif defined(__arm__)
    m0_ptr = (uint8x16_t *)&m0_table[table_offset];
    m1_ptr = (uint8x16_t *)&m1_table[table_offset];


    // even states
    even0_30a  = vqaddq_u8(metrics0_31.val[0],m0_ptr[0]);
    even32_62a = vqaddq_u8(metrics0_31.val[1],m0_ptr[1]);
    even0_30b  = vqaddq_u8(metrics32_63.val[0],m0_ptr[2]);
    even32_62b = vqaddq_u8(metrics32_63.val[1],m0_ptr[3]);

    // odd states
    odd1_31a   = vqaddq_u8(metrics0_31.val[0],m1_ptr[0]);
    odd33_63a  = vqaddq_u8(metrics0_31.val[1],m1_ptr[1]);
    odd1_31b   = vqaddq_u8(metrics32_63.val[0],m1_ptr[2]);
    odd33_63b  = vqaddq_u8(metrics32_63.val[1],m1_ptr[3]);
    // select maxima
    even0_30a  = vmaxq_u8(even0_30a,even0_30b);
    even32_62a = vmaxq_u8(even32_62a,even32_62b);
    odd1_31a   = vmaxq_u8(odd1_31a,odd1_31b);
    odd33_63a  = vmaxq_u8(odd33_63a,odd33_63b);

    // Traceback information
    TBeven0_30  = vceqq_u8(even0_30a,even0_30b);
    TBeven32_62 = vceqq_u8(even32_62a,even32_62b);
    TBodd1_31   = vceqq_u8(odd1_31a,odd1_31b);
    TBodd33_63  = vceqq_u8(odd33_63a,odd33_63b);

    metrics0_31  = vzipq_u8(even0_30a,odd1_31a);
    metrics32_63 = vzipq_u8(even32_62a,odd33_63a);

    TB_ptr[0] = vzipq_u8(TBeven0_30,TBodd1_31);
    TB_ptr[1] = vzipq_u8(TBeven32_62,TBodd33_63);

    in+=2;
    TB_ptr += 2;

    // rescale by subtracting minimum
    /****************************************************
    USE SSSE instruction phminpos!!!!!!!
    ****************************************************/
    min_state =vminq_u8(metrics0_31.val[0],metrics0_31.val[1]);
    min_state =vminq_u8(min_state,metrics32_63.val[0]);
    min_state =vminq_u8(min_state,metrics32_63.val[1]);
    // here we have 16 maximum metrics from the 64 states
    uint8x8_t min_state2 = vpmin_u8(((uint8x8_t*)&min_state)[0],((uint8x8_t*)&min_state)[0]);
    // now the 8 maximum in min_state2
    min_state2 = vpmin_u8(min_state2,min_state2);
    // now the 4 maximum in min_state2, repeated twice
    min_state2 = vpmin_u8(min_state2,min_state2);
    // now the 2 maximum in min_state2, repeated 4 times
    min_state2 = vpmin_u8(min_state2,min_state2);
    // now the 1 maximum in min_state2, repeated 8 times
    min_state  = vcombine_u8(min_state2,min_state2);
    // now the 1 maximum in min_state, repeated 16 times
    metrics0_31.val[0]  = vqsubq_u8(metrics0_31.val[0],min_state);
    metrics0_31.val[1]  = vqsubq_u8(metrics0_31.val[1],min_state);
    metrics32_63.val[0] = vqsubq_u8(metrics32_63.val[0],min_state);
    metrics32_63.val[1] = vqsubq_u8(metrics32_63.val[1],min_state);
#endif
336 337 338 339 340 341 342 343
    }

  } // iteration

  // Traceback
  prev_state0 = 0;
  maxm = 0;

344
#if defined(__x86_64__) || defined(__i386__)
345
  for (s=0; s<16; s++)
346 347
    if (((uint8_t *)&metrics0_15)[s] > maxm) {
      maxm = ((uint8_t *)&metrics0_15)[s];
348 349 350
      prev_state0 = s;
    }

351
  for (s=0; s<16; s++)
352 353
    if (((uint8_t *)&metrics16_31)[s] > maxm) {
      maxm = ((uint8_t *)&metrics16_31)[s];
354 355 356
      prev_state0 = s+16;
    }

357
  for (s=0; s<16; s++)
358 359
    if (((uint8_t *)&metrics32_47)[s] > maxm) {
      maxm = ((uint8_t *)&metrics32_47)[s];
360 361 362
      prev_state0 = s+32;
    }

363
  for (s=0; s<16; s++)
364 365
    if (((uint8_t *)&metrics48_63)[s] > maxm) {
      maxm = ((uint8_t *)&metrics48_63)[s];
366 367 368
      prev_state0 = s+48;
    }

369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395

#elif defined(__arm__)
  for (s=0; s<16; s++)
    if (((uint8_t *)&metrics0_31.val[0])[s] > maxm) {
      maxm = ((uint8_t *)&metrics0_31.val[0])[s];
      prev_state0 = s;
    }

  for (s=0; s<16; s++)
    if (((uint8_t *)&metrics0_31.val[1])[s] > maxm) {
      maxm = ((uint8_t *)&metrics0_31.val[1])[s];
      prev_state0 = s+16;
    }

  for (s=0; s<16; s++)
    if (((uint8_t *)&metrics32_63.val[0])[s] > maxm) {
      maxm = ((uint8_t *)&metrics32_63.val[0])[s];
      prev_state0 = s+32;
    }

  for (s=0; s<16; s++)
    if (((uint8_t *)&metrics32_63.val[1])[s] > maxm) {
      maxm = ((uint8_t *)&metrics32_63.val[1])[s];
      prev_state0 = s+48;
    }
#endif

396
  TB_ptr2 = (uint8_t *)&TB[(n-1)*4];
397

398 399 400
  for (position = n-1 ; position>-1; position--) {

    decoded_bytes[(position)>>3] += (prev_state0 & 0x1)<<(7-(position & 0x7));
401 402 403


    if (TB_ptr2[prev_state0] == 0)
404 405 406 407 408 409
      prev_state0 = (prev_state0 >> 1);
    else
      prev_state0 = 32 + (prev_state0>>1);

    TB_ptr2-=64;
  }
410

411 412

#if defined(__x86_64__) || defined(__i386__)
413 414
  _mm_empty();
  _m_empty();
415
#endif
416 417 418
}

#ifdef TEST_DEBUG
419
int test_viterbi(uint8_t dabflag)
420
{
421 422 423
  uint8_t test[8];
  //_declspec(align(16))  int8_t channel_output[512];
  //_declspec(align(16))  uint8_t output[512],decoded_output[16], *inPtr, *outPtr;
424

425 426 427
  int8_t channel_output[512];
  uint8_t output[512],decoded_output[16], *inPtr, *outPtr;
  uint32_t i;
428

429

430 431 432 433 434 435 436 437
  test[0] = 7;
  test[1] = 0xa5;
  test[2] = 0;
  test[3] = 0xfe;
  test[4] = 0x1a;
  test[5] = 0x33;
  test[6] = 0x99;
  test[7] = 0x14;
438

439 440 441
  if (dabflag==0) {
    ccodelte_init();
    ccodelte_init_inv();
442
  } else {
443 444 445 446
    ccodedab_init();
    ccodedab_init_inv();
    printf("Running with DAB polynomials\n");
  }
447

448 449 450
  inPtr = test;
  outPtr = output;
  phy_generate_viterbi_tables_lte();
451 452 453 454 455
  ccodelte_encode(64,     //input length in bits
                  0,      // add 16-bit crc with rnti scrambling
                  inPtr,  // input pointer
                  outPtr, // output pointer
                  0);     // rnti (optional)
456

457
  for (i = 0; i < 64*3; i++) {
458 459 460 461 462 463
    channel_output[i] = 7*(2*output[i] - 1);
  }

  memset(decoded_output,0,16);
  phy_viterbi_lte_sse2(channel_output,decoded_output,64);
  printf("Optimized Viterbi :");
464

465 466 467 468 469 470 471
  for (i =0 ; i<8 ; i++)
    printf("input %d : %x => %x\n",i,inPtr[i],decoded_output[i]);
}




472 473
int main(int argc, char **argv)
{
474 475

  char c;
476
  uint8_t dabflag=0;
477 478 479 480 481 482 483 484 485 486 487 488 489

  while ((c = getopt (argc, argv, "d")) != -1) {
    if (c=='d')
      dabflag=1;
  }

  test_viterbi(dabflag);
  return(0);
}

#endif // TEST_DEBUG