gfq_ldpc_ntt.c

LDPC编解码器程序

      sum += dfna;
    }
    //normalize
    sum = log2(sum);
    for (a = 0; a < Q; a++) {
      logfna[i][a] = float2fix(_logfna[a] - sum);
    }
  }
  sum = 0;
  for (i = 0; i < Q; i++) {
    if (count[i]) sum += count[i] * log2(count[i]);
  }
  printf("m/n=%g, ", (double)m/n);
  printf("noise entropy = %g bits, rate = %g\n", -sum/n + log2(n),
         (-sum/n + log2(n)) / log2(Q));
  printf("PSNR = %g\n", 10*log((Q-1)*(Q-1)/mse(x,y,n)));
}

void enc(int x[], int s[])
{
  int i, j;
  for (j = 0; j < m; j++) {
    int sum = 0;
    for (i = 0; i < row_weight[j]; i++) {
      int xHmn = GF_mul(x[row_col[j][i]], H[j][i]);
      sum = GF_add(sum, xHmn);
    }
    s[j] = sum;
  }
}

NTT **malloc2DNTT(int a, int b) // allocates array[a][b]
{
  int i;
  NTT **pp = malloc(sizeof(NTT *) * a);
  NTT *p = malloc(sizeof(NTT) * a * b);
  if (pp == NULL || p == NULL) exit(-1);
  for (i = 0; i < a; i++) {
    pp[i] = p + b*i;
  }
  return pp;
}

int ***malloc2Dintp(int a, int b) // allocates array[a][b] (int pointer)
{
  int j;
  int ***ppp= malloc(sizeof(int **) * a);
  int **pp  = malloc(sizeof(int *) * a * b);
  if (ppp == NULL || pp == NULL) exit(-1);
  for (j = 0; j < a; j++) {
    ppp[j] = pp + b*j;
  }
  return ppp;
}

int **malloc2Dint(int a, int b) // allocates array[a][b]
{
  int i;
  int **pp = malloc(sizeof(int *) * a);
  int *p = malloc(sizeof(int) * a * b);
  if (pp == NULL || p == NULL) exit(-1);
  for (i = 0; i < a; i++) {
    pp[i] = p + b*i;
  }
  return pp;
}

int ***malloc3Dint(int a, int b, int c) // allocates array[a][b][c]
{
  int i, j;
  int ***ppp= malloc(sizeof(int **) * a);
  int **pp  = malloc(sizeof(int *) * a * b);
  int *p    = malloc(sizeof(int) * a * b * c);
  if (ppp == NULL || pp == NULL || p == NULL) exit(-1);
  for (j = 0; j < a; j++) {
    for (i = 0; i < b; i++) {
      pp[i+b*j] = p + c*i + c*b*j;
    }
    ppp[j] = pp + b*j;
  }
  return ppp;
}

// Sum Product Decoder
// logfna: prior
// z: syndrome
// loop_max: max iteration
// x[]: original signal (just for reference)
int dec(int **logfna, int z[], int loop_max, int x[])
{
  int a, i, j, k, loop;
  int iir, prev = 999999, nodecr = 0;

  for (i = 0; i < n; i++) {
    for (k = 0; k < col_weight[i]; k++) {
      memcpy(logqacol[i][k], logfna[i], sizeof(int)*Q);
    }
  }

  for (loop = 0; loop < loop_max; loop++) {
    for (j = 0; j < m; j++) {
      int row_weightj = row_weight[j];
      int logprodqa[Q], sgnsum[Q];
      memset(sgnsum, 0, sizeof(sgnsum));
      memset(logprodqa, 0, sizeof(logprodqa));
      memset(*isnegative, 0, rmax * Q * sizeof(int));
      for (k = 0; k < row_weightj; k++) {
        for (a = 0; a < Q; a++) {
          fQa[k][GF_mul(a, H[j][k])] = Fexp(logqa[j][k][a]);
        }
        ntt(fQa[k]);
        //Flog(|fQa[k][a]|)傪偳偙偐偵曐懚偟屻偵巊偆偺傕傾儕
        for (a = 0; a < Q; a++) {
          if (fQa[k][a] < 0) {
            fQa[k][a] = -fQa[k][a];
            isnegative[k][a] = 1;
            sgnsum[a] ^= 1;     // 晧偺偲偒偩偗僇僂儞僩傾僢僾
          }
          fQa[k][a] = (fQa[k][a] + (1<<(Log2Q/2-1))) >> (Log2Q/2);
          logprodqa[a] += Flog(fQa[k][a]);
          assert(logprodqa[a] > -(1<<29));
        }
      }
      for (k = 0; k < row_weightj; k++) {
        NTT fRa[Q];
        for (a = 0; a < Q; a++) {
          fRa[a] = Fexp(logprodqa[a] - Flog(fQa[k][a]) + Log2Q*FMUL); //prod. except k
          if (isnegative[k][a] != sgnsum[a])
            fRa[a] = -fRa[a];
        }
        ntt(fRa);          // intt
        for (a = 0; a < Q; a++) {
          logra[j][k][a] = Flog(fRa[GF_sub(z[j], GF_mul(a, H[j][k]))]);
        }
      }
    }
    // update qa & tentative decode
    for (i = 0; i < n; i++) {
      int col = col_weight[i];
      int argmaxa = -1;
      int maxp = -(1<<30);
      for (a = 0; a < Q; a++) {
        int logprod = logfna[i][a];
        for (k = 0; k < col; k++) {
          logprod += logracol[i][k][a];
        }
        for (k = 0; k < col; k++) {
          logqacol[i][k][a] = logprod - logracol[i][k][a]; //taka product except k
        }
        if (logprod > maxp) {
          maxp = logprod;
          argmaxa = a;
        }
      }
      tmp_x[i] = argmaxa;
    }
    //normalize qa
    for (j = 0; j < m; j++) {
      for (k = 0; k < row_weight[j]; k++) {
        int *qajk = logqa[j][k];
        int sum = 0;
        for (a = 0; a < Q; a++) {
          sum += Fexp(*qajk++);
        }
        if (sum <= 0) {//all zero happenes
          //printf("%d,%d:sum%d ", j, k, sum);
          sum = Q;//taka051207
        }
        sum = Flog(sum);
        qajk = logqa[j][k];
        for (a = 0; a < Q; a++) {
          *qajk++ -= sum;
        }
      }
    }
    printf("%2d: HamDist(x,x^)=%d ", loop+1, HamDist(x, tmp_x, n));

    enc(tmp_x, tmp_z);
    {
      int dist = HamDist(z, tmp_z, m);
      printf("HamDist(s,synd(x^))=%d\n", dist);
      if (dist == 0)           // nothing more can be done
        return 0;
      // nonconvergence detection
      if (loop == 0) iir = dist;
      else iir = (int)(iir * 0.85 + dist * 0.15 + 0.5);

      if (prev <= dist) nodecr++;
      else nodecr = 0;
      if (dist > iir * 1.1 || nodecr > 10) break; // no conversion
      prev = dist;
    }
  }

  return -1;
}

//http://www.inference.phy.cam.ac.uk/mackay/codes/alist.html
void initdec(char *s)
{
  int i, j, q, *count;
  int isbinary = 0;
  char buf[BUFSIZ];
  FILE *fp = fopen(s, "rt");
  if (fp == NULL) {
    fprintf(stderr, "cannot open %s\n", s);
    exit(-2);
  }
  fgets(buf, sizeof(buf), fp);
  if (sscanf(buf, "%d%d%d", &n, &m, &q) == 2) {
    fprintf(stderr, "Warning! A binary matrix seems to be specified.\n");
    fprintf(stderr, "I'll use it anyways.\n");
    isbinary = 1;
    q = 2;
  }
  if (q > Q) exit(-1);
  if (q != Q) {
    fprintf(stderr, "Warning! GF order(%d) does not match with matrix file(%d)\n", Q, q);
    fprintf(stderr, "I'll use it anyways.\n");
  }
  fscanf(fp, "%d%d", &cmax, &rmax);
  col_weight = malloc(sizeof(int) * n);
  for (i = 0; i < n; i++) {
    fscanf(fp, "%d", &col_weight[i]);
  }
  row_weight = malloc(sizeof(int) * m);
  for (j = 0; j < m; j++)
    fscanf(fp, "%d", &row_weight[j]);

  count = malloc(sizeof(int) * n);
  memset(count, 0, sizeof(int) * n);

  row_col = malloc2Dint(m, rmax);
  row_N   = malloc2Dint(m, rmax);
  col_row = malloc2Dint(n, cmax);
  col_N   = malloc2Dint(n, cmax);
  H            = malloc2Dint(m, rmax); // Hmn
  fQa          = malloc2DNTT(rmax, Q);
  isnegative   = malloc2Dint(rmax, Q);

  {//skip n lines
    for (i = 0; i < n; i++) {
      for (j = 0; j < cmax; j++)
        fscanf(fp, isbinary ? "%*d" : "%*d%*d");
    }
  }

  for (j = 0; j < m; j++) {
    for (i = 0; i < row_weight[j]; i++) {
      int v;
      if (isbinary) {
        fscanf(fp, "%d", &v);
        H[j][i] = 1;
      } else
        fscanf(fp, "%d%d", &v, &H[j][i]);
      v--;
      row_col[j][i] = v;
      row_N[j][i] = count[v];
      col_row[v][count[v]] = j;
      col_N[v][count[v]] = i;
      count[v]++;
    }

    for (; i < rmax; i++) {//skip fillers
      int a,b;
      fscanf(fp, "%d%d", &a, &b);
      if (a!=0 || b!=0) {
        printf("error at row %d, %d %d\n", i, a, b);
        exit(-1);
      }
    }
  }
  free(count);
  fclose(fp);
  tmp_z = malloc(sizeof(int) * m);
  tmp_x = malloc(sizeof(int) * n);

  logra  = malloc3Dint(m, rmax, Q);
  logqa  = malloc3Dint(m, rmax, Q);
  logracol=malloc2Dintp(n, cmax);
  logqacol=malloc2Dintp(n, cmax);
  for (i = 0; i < n; i++) {
    int k;
    int *cri = col_row[i], *cNi = col_N[i];
    int col = col_weight[i];
    for (k = 0; k < col; k++) {
      logracol[i][k] = logra[cri[k]][cNi[k]];
      logqacol[i][k] = logqa[cri[k]][cNi[k]];
    }
  }
}

int main(int argc, char **argv)
{
  int i, j, *s, *s0, *x, *y;
  int **logfna;
  double error;
  void (*channel)(int x[], int y[], double p, int **logfna) = lap;
  int iteration = 1000;

  initlogexptab2();

  if (argc == 6 && strcmp(argv[1], "-iter") == 0) {
    iteration = atoi(argv[2]);
    argc -= 2;
    argv += 2;
  }
  if (argc != 4) {
  usage:
    fprintf(stderr, "GFq_LDPC_NTT - GF(%d) simulator\n", Q);
    fprintf(stderr, "usage: GFq_LDPC_NTT [-iter num] <lap|bsc> NoiseLevel ParityMatrix\n");
    fprintf(stderr, "e.g.: GFq_LDPC_NTT lap 0.5 q8.sp.6000.4000.3000.1\n");
    fprintf(stderr, "      GFq_LDPC_NTT bsc 0.08 q4.sp.9000.6000.4500.1\n");
    return -1;
  }
  if (strcmp(argv[1], "bsc") == 0) channel = bsc;
  else if (strcmp(argv[1], "lap") != 0) {
    fprintf(stderr, "please specify 'lap' or 'bsc'\n");
    goto usage;
  }
  error = atof(argv[2]);
  initdec(argv[3]);
  s = malloc(sizeof(int) * m);  // syndrome
  s0= malloc(sizeof(int) * m);  // syndrome
  x = malloc(sizeof(int) * n);  // source
  y = malloc(sizeof(int) * n);  // side information
  logfna=malloc2Dint(n, Q);
  // iterate experiments
  for (j = 1; j <= 3; j++) {
    printf("\nGF(%d) experiment %d:\n", Q, j);
    SRand(j);
    for (i = 0; i < n; i++)
      x[i] = Rand() % Q;
    enc(x, s);
    channel(x, y, error, logfna);
    enc(y, s0);
    printf("start: HamDist(x,y)=%d HamDist(synd(x),synd(y))=%d\n",
           HamDist(x,y,n), HamDist(s,s0,m));
    puts(
      dec(logfna, s, iteration, x) ?
      "failed." : "converged."
      );
  }
  return 0;
}