tcc_dvb-rcs_pc.c

MAP 算法 源程序 第三代 数字通信 turbo码相关

    }
  }

  for (i=0; i<2; i++){                /* turbo code permutation Level two */
    for (j=0; j<N; j++)
      permu[i][j] = *(pinfo_array + 2*j + i);
  }

  for (i=0; i<2; i++){
    for (j=0; j<N; j++)
      *(pinfo_array + 2*j + i) = permu[i][*(alpha + j)];
  }

  #if R == 1                          /* code rate = 1/3 */
  for (j=0; j<N; j++){
    for (i=0; i<2; i++){
      yy[i][j] = *(in_array + 2*N + 2*j + i);
      ww[i][j] = *(in_array + 4*N + 2*j + i);
    }
  }
  #endif

  #if R == 2                          /* code rate = 2/5 */
  for (j=0; j<N; j++){
    for (i=0; i<2; i++)
      yy[i][j] = *(in_array + 2*N + 2*j + i);
  }
  for (j=0; j<M; j++){
    for (i=0; i<2; i++){
     ww[i][2*j] = *(in_array + 4*N + 2*j + i);
     ww[i][2*j+1] = 0;
    }
  }
  #endif

  #if R == 3                          /* code rate = 1/2 */
  for (j=0; j<N; j++){
    for (i=0; i<2; i++)
      yy[i][j] = *(in_array + 2*N + 2*j + i);
  }
  #endif

  #if R == 4                          /* code rate = 2/3 */
  for (j=0; j<M; j++){
    for (i=0; i<2; i++){
     yy[i][2*j] = *(in_array + 2*N + 2*j + i);
     yy[i][2*j+1] = 0;
    }
  }
  #endif

  #if R == 5                          /* code rate = 3/4 */
  if (N%3 == 0){         /* N modulo 3 = 1, M=N/3 */
    for (j=0; j<M; j++){
      for (i=0; i<2; i++)
	yy[i][3*j] = *(in_array + 2*N + 2*j + i);
      for (i=1; i<3; i++)
	yy[i-1][3*j+i] = 0;
    }
  }
  else if (N%3 == 1){    /* N modulo 3 = 1, M=(N-4)/3+2 */
    for (j=0; j<(M-1); j++){
      for (i=0; i<2; i++)
	yy[i][3*j] = *(in_array + 2*N + 2*j + i);
      for (i=1; i<3; i++) 
	yy[i-1][3*j+i] = 0;
    }
    yy[0][N-1] = *(in_array + 2*N + 2*M - 2);
    yy[1][N-1] = *(in_array + 2*N + 2*M - 1);
  }
  else if (N%3 == 2){    /* N modulo 3 = 2, M=(N-8)/3+3 */
    for (j=0; j<(M-1); j++){
      for (i=0; i<2; i++)
	yy[i][3*j] = *(in_array + 2*N + 2*j + i);
      for (i=1; i<3; i++)
	yy[i-1][3*j+i] = 0;
    }
    yy[0][N-2] = *(in_array + 2*N + 2*M - 2);
    yy[0][N-1] = 0;
    yy[1][N-2] = *(in_array + 2*N + 2*M - 1);
    yy[1][N-1] = 0;
  }
  #endif

  #if R == 6                         /* code rate = 4/5 */
  for (j=0; j<M; j++){
    for (i=0; i<2; i++)
      yy[i][4*j] = *(in_array + 2*N + 2*j + i);
    for (i=1; i<4; i++){
      yy[0][4*j+i] = 0;
      yy[1][4*j+i] = 0;
    }
  }
  #endif

  #if R == 7                        /* code rate = 6/7 */
  if (N%3 == 0){         /* N modulo 3 = 1, M=N/6 */
    for (j=0; j<M; j++){
      for (i=0; i<2; i++)
	yy[i][6*j] = *(in_array + 2*N + 2*j + i);
      for (i=1; i<6; i++){
	yy[0][6*j+i] = 0;
	yy[1][6*j+i] = 0;
      }
    }
  }
  else if (N%3 == 1){    /* N modulo 3 = 1, M=(N-4)/6+1 */
    for (j=0; j<(M-1); j++){
      for (i=0; i<2; i++)
	yy[i][6*j] = *(in_array + 2*N + 2*j + i);
      for (i=1; i<6; i++){
	yy[0][6*j + i] = 0;
	yy[1][6*j + i] = 0;
      }
    }
    yy[0][N-4] = *(in_array + 2*N + 2*M - 2);
    for (j=3; j>0; j--)
      yy[0][N-j] = 0;
    yy[1][N-4] = *(in_array + 2*N + 2*M - 1);
    for (j=3; j>0; j--)
      yy[1][N-j] = 0;
  }
  else if (N%3 == 2){    /* N modulo 3 = 2, M=(N-8)/6+2 */
    for (j=0; j<(M-1); j++){
      for (i=0; i<2; i++)
	yy[i][6*j] = *(in_array + 2*N + 2*j + i);
      for (i=1; i<6; i++){
	yy[0][6*j + i] = 0;
	yy[1][6*j + i] = 0;
      }
    }
    yy[0][N-2] = *(in_array + 2*N + 2*M - 2);
    yy[0][N-1] = 0;
    yy[1][N-2] = *(in_array + 2*N + 2*M - 1);
    yy[1][N-1] = 0;
  }
  #endif

  if (R >= 3){
    for (j=0; j<N; j++){
      for (i=0; i<2; i++)
        ww[i][j] = 0;
    }
  }

  for (j=0; j<N; j++){
    *(parity1_array + 2*j) = ww[0][j];
    *(parity1_array + 2*j + 1) = yy[0][j];
    *(parity2_array + 2*j) = ww[1][j];
    *(parity2_array + 2*j + 1) = yy[1][j];
  }
}

/*******************************************/
/* Non-binary turbo code permutation table */
/*******************************************/
void permutation(int *alpha)
{
  int j;         /* loop variables */

  #if N == 48
  int p[4] = {11,24,0,24};
  #endif

  #if N == 64
  int p[4] = {7,34,32,2};
  #endif
 
  #if N == 212
  int p[4] = {13,106,108,2};
  #endif

  #if N == 220
  int p[4] = {23,112,4,116};
  #endif

  #if N == 228
  int p[4] = {17,116,72,188};
  #endif

  #if N == 424
  int p[4] = {11,6,8,2};
  #endif

  #if N == 432
  int p[4] = {13,0,4,8};
  #endif

  #if N == 440
  int p[4] = {13,10,4,2};
  #endif

  #if N == 752
  int p[4] = {19,376,224,600};
  #endif

  #if N == 848
  int p[4] = {19,2,16,6};
  #endif

  #if N == 856
  int p[4] = {19,428,224,652};
  #endif

  #if N == 864
  int p[4] = {19,2,16,6};
  #endif

  for (j=0; j<N; j++){
    if (j%4 == 0) *(alpha + j) = (p[0]*j + 1) % N;
    else if (j%4 == 1) *(alpha + j) = (p[0]*j + N/2 + p[1] + 1) % N;
    else if (j%4 == 2) *(alpha + j) = (p[0]*j + p[2] + 1) % N;
    else if (j%4 == 3) *(alpha + j) = (p[0]*j + N/2 + p[3] + 1) % N;
  }
}

/****************************************/
/* Turbo code decoding with Max_log_map */
/****************************************/
void trellis(int g[3][K], int nextstate[NUMSTATES][K], int nextoutput[NUMSTATES][NUMSTATES],
	     int prestate[NUMSTATES][K], int preoutput[NUMSTATES][NUMSTATES]);
float max4(float x, float y, float z, float w);
float max8(float x, float y, float z, float w, float A, float B, float C, float D);

void maxlogmap(int g[3][K], float *rec_s, float *rec_p, float *L_a, float *L_u, float *Alpha0, float *BetaN)
{
  int i,j,l,ll,t;                         /* loop variables */
  int nextoutput[NUMSTATES][NUMSTATES];   /* maps current/next states to next output */
  int nextstate[NUMSTATES][K];            /* for current state, gives next for given input */
  int preoutput[NUMSTATES][NUMSTATES];    /* maps previous/current states to previous output */
  int prestate[NUMSTATES][K];             /* gives conv. encoder output */
  int binary_input[2],binary_output[2];   /* branch input and output */

  float gamma[NUMSTATES] = {0.0};         /* the branch transition probabilities */
  float Alpha[N+1][NUMSTATES];
  float Beta[N+1][NUMSTATES];

  float Alpha_temp[4] = {0.0};            /* temporary values of Alpha or Beta */
  float Beta_temp[4] = {0.0};
  float tempmax[N] = {0.0};


  float temp0[NUMSTATES] = {0.0};
  float temp1[NUMSTATES] = {0.0};
  float temp2[NUMSTATES] = {0.0};
  float temp3[NUMSTATES] = {0.0};
  float max8temp0 = {0.0};
  float gamma0[4] = {0.0};
  float p[N][4] = {0.0};

  for (i = 0; i < NUMSTATES; i++){        /* initialize Alpha and Beta data array */
    Alpha[0][i] = *(Alpha0 + i);
    Beta[N][i] = *(BetaN + i);
  }

  for (j=0; j<N; j++){
    p[j][0] = -max4(0,*(L_a+j),*(L_a+j+N),*(L_a+j+2*N));
    p[j][1] = *(L_a+j) - max4(0,*(L_a+j),*(L_a+j+N),*(L_a+j+2*N));
    p[j][2] = *(L_a+j+N) - max4(0,*(L_a+j),*(L_a+j+N),*(L_a+j+2*N));
    p[j][3] = *(L_a+j+2*N) - max4(0,*(L_a+j),*(L_a+j+N),*(L_a+j+2*N));
  }

  trellis(g, nextstate, nextoutput, prestate, preoutput);

  for (t=1; t<(N+1); t++) {        /* trace forward, compute Alpha */
    for (j=0; j<NUMSTATES; j++) {
      for (i=0; i<4; i++){
	deci2bin(i, 2, binary_input);
	deci2bin(preoutput[j][prestate[j][i]], 2, binary_output);
	for (ll=0; ll<2; ll++){
	  binary_input[ll] = 1 - 2*binary_input[ll];
	  binary_output[ll] = 1 - 2*binary_output[ll];
	}
        gamma[prestate[j][i]] = *(rec_s+2*t-2)*binary_input[1] + *(rec_s+2*t-1)*binary_input[0]
                                + *(rec_p+2*t-2)*binary_output[1] + *(rec_p+2*t-1)*binary_output[0] + p[t-1][i];
	Alpha_temp[i] = gamma[prestate[j][i]] + Alpha[t-1][prestate[j][i]];
      }
      Alpha[t][j] = max4(Alpha_temp[0], Alpha_temp[1], Alpha_temp[2], Alpha_temp[3]);
    }
    tempmax[t-1] = max8(Alpha[t][0],Alpha[t][1],Alpha[t][2],Alpha[t][3],Alpha[t][4],Alpha[t][5],Alpha[t][6],Alpha[t][7]);

    for (j=0; j<NUMSTATES; j++)    /* to be sure that Alpha do not overflow */
      Alpha[t][j] = Alpha[t][j] - tempmax[t-1];
  }

  for (t=N-1; t>=0; t--) {         /* trace forward, compute Beta */
    for (j=0; j<NUMSTATES; j++) {
      for (i=0; i<4; i++){
	deci2bin(i, 2, binary_input);
	deci2bin(nextoutput[j][nextstate[j][i]], 2, binary_output);
	for (ll=0; ll<2; ll++){
	  binary_input[ll] = 1 - 2*binary_input[ll];
	  binary_output[ll] = 1 - 2*binary_output[ll];
	}
        gamma[nextstate[j][i]] = *(rec_s+2*t)*binary_input[1] + *(rec_s+2*t+1)*binary_input[0]
                                + *(rec_p+2*t)*binary_output[1] + *(rec_p+2*t+1)*binary_output[0] + p[t][i];
	Beta_temp[i] = gamma[nextstate[j][i]] + Beta[t+1][nextstate[j][i]];
      }
      Beta[t][j] = max4(Beta_temp[0], Beta_temp[1], Beta_temp[2], Beta_temp[3]);
    }

    for (j=0; j<NUMSTATES; j++)    /* to be sure that Beta do not overflow */
      Beta[t][j] = Beta[t][j] - tempmax[t];
  }

  for (t=0; t<N; t++) {            /* compute the soft output, log_likelihood ratio */
    for (j=0; j<NUMSTATES; j++) {
      for (i=0; i<4; i++){
	deci2bin(i, 2, binary_input);
	deci2bin(preoutput[j][prestate[j][i]], 2, binary_output);
	for (ll=0; ll<2; ll++){
	  binary_input[ll] = 1 - 2*binary_input[ll];
	  binary_output[ll] = 1 - 2*binary_output[ll];
	}
        gamma0[i] = *(rec_s+2*t)*binary_input[1] + *(rec_s+2*t+1)*binary_input[0]
                   + *(rec_p+2*t)*binary_output[1] + *(rec_p+2*t+1)*binary_output[0] + p[t][i];
      }

      temp0[j] = gamma0[0] + Alpha[t][prestate[j][0]] + Beta[t+1][j]; 

      temp1[j] = gamma0[1] + Alpha[t][prestate[j][1]] + Beta[t+1][j];

      temp2[j] = gamma0[2] + Alpha[t][prestate[j][2]] + Beta[t+1][j];

      temp3[j] = gamma0[3] + Alpha[t][prestate[j][3]] + Beta[t+1][j];
    }
    max8temp0 = max8(temp0[0],temp0[1],temp0[2],temp0[3],temp0[4],temp0[5],temp0[6],temp0[7]);

    *(L_u+t) = max8(temp1[0],temp1[1],temp1[2],temp1[3],temp1[4],temp1[5],temp1[6],temp1[7]) - max8temp0;

    *(L_u+t+N) = max8(temp2[0],temp2[1],temp2[2],temp2[3],temp2[4],temp2[5],temp2[6],temp2[7]) - max8temp0;

    *(L_u+t+2*N) = max8(temp3[0],temp3[1],temp3[2],temp3[3],temp3[4],temp3[5],temp3[6],temp3[7]) - max8temp0;
  }

  for (i = 0; i < NUMSTATES; i++){
    *(Alpha0 + i) = Alpha[N][i];
    *(BetaN + i) = Beta[0][i];
  }
}

float max4(float x, float y, float z, float w)
{
  float maxmum;

  maxmum = x>y?x:y;
  maxmum = maxmum>z?maxmum:z;
  maxmum = maxmum>w?maxmum:w;

  return(maxmum);
}

float max8(float x, float y, float z, float w, float A, float B, float C, float D)
{
  float maxmum;

  maxmum = x>y?x:y;
  maxmum = maxmum>z?maxmum:z;
  maxmum = maxmum>w?maxmum:w;
  maxmum = maxmum>A?maxmum:A;
  maxmum = maxmum>B?maxmum:B;
  maxmum = maxmum>C?maxmum:C;
  maxmum = maxmum>D?maxmum:D;

  return(maxmum);
}

/*************************************************/
/* Set up the trellis for a given code generator */
/*************************************************/
void deci2bin(int d, int size, int *b);
int bin2deci(int *b, int size);
int nxt_stat(int g[3][K], int current_state, int input, int *memory_contents);

void trellis(int g[3][K], int nextstate[NUMSTATES][K], int nextoutput[NUMSTATES][NUMSTATES],
	     int prestate[NUMSTATES][K], int preoutput[NUMSTATES][NUMSTATES])
{
  int i,j,l;                   /* loop variables */
  int memory_contents[K];      /* feedback bit + conv. encoder sr */
  int encd_output[2];          /* store trial encoder output */
  int next_state;

  for (i = 0; i < NUMSTATES; i++) {     /* initialize data structures */
    for (j = 0; j < NUMSTATES; j++){
      preoutput[i][j] = 0;
      nextoutput[i][j] = 0;
    }

    for (j = 0; j < K; j++){
      prestate[i][j] = 0;
      nextstate[i][j] = 0;
    }
  }

  for (j = 0; j < NUMSTATES; j++) {     /* generate the matrices */
    for (l = 0; l < K; l++) {
      next_state = nxt_stat(g, j, l, memory_contents);

      encd_output[0] = 0;
      encd_output[1] = 0;
      for (i = 0; i < K; i++){
        encd_output[0] ^= memory_contents[i] & g[1][i];
        encd_output[1] ^= memory_contents[i] & g[2][i];
      }

      nextstate[j][l] = next_state;                        /* next state, given current state and input */
      nextoutput[j][next_state] = bin2deci(encd_output,2); /* generate the next output */

      prestate[next_state][l] = j;
      preoutput[next_state][j] = nextoutput[j][next_state];
    }  /* end of l for loop */
  }  /* end of j for loop */
}

/* ************************************************************************** */
int nxt_stat(int g[3][K], int current_state, int input, int *memory_contents) {

  int binary_state[K - 1];              /* binary value of current state */
  int binary_input[2];                  /* binary value of input */
  int next_state_binary[K - 1];         /* binary value of next state */
  int next_state;                       /* decimal value of next state */
  int feedback;                         /* feedback bit */
  int i;                                /* loop variable */

  deci2bin(current_state, K - 1, binary_state);  /* convert the decimal value of the current state number to binary */
  deci2bin(input, 2, binary_input);              /* convert the decimal value of the input to binary */

  feedback = binary_input[0] ^ binary_input[1];
  for (i=0; i<K-1; i++){
    feedback ^= binary_state[i] & g[0][i+1];     /* get the feedback bit */
    next_state_binary[i] = 0;
  }

  next_state_binary[0] = feedback;       /* given the input and current state number, compute the next state number */
  for (i = 1; i < K - 1; i++)
    next_state_binary[i] = binary_state[i - 1] ^ binary_input[0];
  next_state = bin2deci(next_state_binary, K - 1);  /* convert the binary value of the next state number to decimal */

  memory_contents[0] = feedback;         /* memory_contents are feedback bits in the encoder */
  for (i = 1; i < K; i++)
    memory_contents[i] = binary_state[i - 1];

  return(next_state);
}

void deci2bin(int d, int size, int *b) {  /* converts a decimal number to a binary number */
   int i;

   for(i = 0; i < size; i++)
     b[i] = 0;

   b[size - 1] = d & 0x01;

   for (i = size - 2; i >= 0; i--) {
     d = d >> 1;
     b[i] = d & 0x01;
   }
}

int bin2deci(int *b, int size) {         /* converts a binary number to decimal number */
   int i, d;

   d = 0;

   for (i = 0; i < size; i++)
     d = d + (int)pow(2, size - i - 1) * b[i];

   return(d);
}