Added Test and Benchmark for soft_demodulation approximation algorithms

lte/phy/include/liblte/phy/channel/ch_awgn.h

@@ -27,6 +27,8 @@
 
 
 #include <complex.h>
+#include <stdint.h>
+
 #include "liblte/config.h"
 
 #ifndef CH_AWGN_
@@ -37,13 +39,15 @@ typedef _Complex float cf_t;
 LIBLTE_API void ch_awgn_c(const cf_t* input, 
                           cf_t* output, 
                           float variance, 
-                          int buff_sz);
+                          uint32_t len);
 
 LIBLTE_API void ch_awgn_f(const float* x, 
                           float* y, 
                           float variance, 
-                          int buff_sz);
+                          uint32_t len);
 
+LIBLTE_API float ch_awgn_get_variance(float ebno_db, 
+                                      float rate);
 
 /* High-level API */
 

lte/phy/include/liblte/phy/utils/vector.h

@@ -35,6 +35,9 @@
 
 typedef _Complex float cf_t;
 
+#define EXPAVERAGE(data, average, nframes) ((data + average * nframes) / (nframes + 1))  
+
+
 /** Return the sum of all the elements */
 LIBLTE_API int vec_acc_ii(int *x, uint32_t len);
 LIBLTE_API float vec_acc_ff(float *x, uint32_t len);

lte/phy/lib/channel/src/ch_awgn.c

@@ -29,25 +29,31 @@
 #include <complex.h>
 #include <stdlib.h>
 #include <strings.h>
+#include <math.h>
 
 #include "gauss.h"
 #include "liblte/phy/channel/ch_awgn.h"
 
-void ch_awgn_c(const cf_t* x, cf_t* y, float variance, int buff_sz) {
+float ch_awgn_get_variance(float ebno_db, float rate) {
-  _Complex float tmp;
+  float esno_db = ebno_db + 10 * log10f(rate);
-  int i;
+  return sqrtf(1 / (powf(10, esno_db / 10)));
+}
+
+void ch_awgn_c(const cf_t* x, cf_t* y, float variance, uint32_t len) {
+  cf_t tmp;
+  uint32_t i;
 
-  for (i=0;i<buff_sz;i++) {
+  for (i=0;i<len;i++) {
     __real__ tmp = rand_gauss();
     __imag__ tmp = rand_gauss();
     tmp *= variance;
     y[i] = tmp + x[i];
   }
 }
-void ch_awgn_f(const float* x, float* y, float variance, int buff_sz) {
+void ch_awgn_f(const float* x, float* y, float variance, uint32_t len) {
-  int i;
+  uint32_t i;
 
-  for (i=0;i<buff_sz;i++) {
+  for (i=0;i<len;i++) {
     y[i] = x[i] + variance * rand_gauss();
   }
 }

lte/phy/lib/modem/src/demod_soft.c

@@ -36,6 +36,7 @@
 
 void demod_soft_init(demod_soft_t *q) {
   bzero((void*)q,sizeof(demod_soft_t));
+  q->sigma = 1.0; 
 }
 
 void demod_soft_table_set(demod_soft_t *q, modem_table_t *table) {

lte/phy/lib/modem/src/soft_algs.c

@@ -50,8 +50,10 @@
  * \param S Soft demapping auxiliary matrix
  * \param sigma2 Noise vatiance
  */
-void llr_approx(const _Complex float *in, float *out, int N, int M, int B,
+void
-    _Complex float *symbols, uint32_t (*S)[6][32], float sigma2) {
+llr_approx(const _Complex float *in, float *out, int N, int M, int B,
+           _Complex float *symbols, uint32_t(*S)[6][32], float sigma2)
+{
   int i, s, b;
   float num, den;
   int change_sign = -1;
@@ -106,8 +108,10 @@ void llr_approx(const _Complex float *in, float *out, int N, int M, int B,
  * \param S Soft demapping auxiliary matrix
  * \param sigma2 Noise vatiance
  */
-void llr_exact(const _Complex float *in, float *out, int N, int M, int B,
+void
-    _Complex float *symbols, uint32_t (*S)[6][32], float sigma2) {
+llr_exact(const _Complex float *in, float *out, int N, int M, int B,
+          _Complex float *symbols, uint32_t(*S)[6][32], float sigma2)
+{
   int i, s, b;
   float num, den;
   int change_sign = -1;

lte/phy/lib/modem/test/CMakeLists.txt

@@ -36,5 +36,13 @@ ADD_TEST(modem_qpsk_soft modem_test -n 1020 -m 2 -s)
 ADD_TEST(modem_qam16_soft modem_test -n 1020 -m 4 -s)
 ADD_TEST(modem_qam64_soft modem_test -n 1020 -m 6 -s)
  
+ADD_EXECUTABLE(soft_demod_test soft_demod_test.c)
+TARGET_LINK_LIBRARIES(soft_demod_test lte_phy)
+
+ADD_TEST(modem_bpsk_soft modem_test -n 1020 -m 1) 
+ADD_TEST(modem_qpsk_soft modem_test -n 1020 -m 2)
+ADD_TEST(modem_qam16_soft modem_test -n 1020 -m 4)
+ADD_TEST(modem_qam64_soft modem_test -n 1020 -m 6)
+ 
 
 

lte/phy/lib/modem/test/soft_demod_test.c

@@ -0,0 +1,229 @@
+/**
+ *
+ * \section COPYRIGHT
+ *
+ * Copyright 2013-2014 The libLTE Developers. See the
+ * COPYRIGHT file at the top-level directory of this distribution.
+ *
+ * \section LICENSE
+ *
+ * This file is part of the libLTE library.
+ *
+ * libLTE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of
+ * the License, or (at your option) any later version.
+ *
+ * libLTE 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 Lesser General Public License for more details.
+ *
+ * A copy of the GNU Lesser General Public License can be found in
+ * the LICENSE file in the top-level directory of this distribution
+ * and at http://www.gnu.org/licenses/.
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <strings.h>
+#include <unistd.h>
+#include <math.h>
+#include <time.h>
+#include <stdbool.h>
+
+#include "liblte/phy/phy.h"
+
+int nof_frames = 10; 
+int num_bits = 1000;
+lte_mod_t modulation = 0;
+
+void usage(char *prog) {
+  printf("Usage: %s [nfv] -m modulation (1: BPSK, 2: QPSK, 4: QAM16, 6: QAM64)\n", prog);
+  printf("\t-n num_bits [Default %d]\n", num_bits);
+  printf("\t-f nof_frames [Default %d]\n", nof_frames);
+  printf("\t-v verbose [Default None]\n");
+}
+
+void parse_args(int argc, char **argv) {
+  int opt;
+  while ((opt = getopt(argc, argv, "nmvf")) != -1) {
+    switch (opt) {
+    case 'n':
+      num_bits = atoi(argv[optind]);
+      break;
+    case 'f':
+      nof_frames = atoi(argv[optind]);
+      break;
+    case 'v':
+      verbose++;
+      break;
+    case 'm':
+      switch(atoi(argv[optind])) {
+      case 1:
+        modulation = LTE_BPSK;
+        break;
+      case 2:
+        modulation = LTE_QPSK;
+        break;
+      case 4:
+        modulation = LTE_QAM16;
+        break;
+      case 6:
+        modulation = LTE_QAM64;
+        break;
+      default:
+        fprintf(stderr, "Invalid modulation %d. Possible values: "
+            "(1: BPSK, 2: QPSK, 4: QAM16, 6: QAM64)\n", atoi(argv[optind]));
+        break;
+      }
+      break;
+    default:
+      usage(argv[0]);
+      exit(-1);
+    }
+  }
+  if (modulation == 0) {
+    usage(argv[0]);
+    exit(-1);
+  }
+}
+
+float mse_threshold() {
+  switch(modulation) {
+    case LTE_BPSK: 
+      return 1.0e-6;
+    case LTE_QPSK:
+      return 1.0e-6; 
+    case LTE_QAM16: 
+      return 0.11; 
+    case LTE_QAM64:
+      return 0.18;
+    default:
+      return -1.0;
+  }
+}
+
+int main(int argc, char **argv) {
+  int i;
+  modem_table_t mod;
+  demod_soft_t demod_soft;
+  char *input, *output;
+  cf_t *symbols;
+  float *llr_exact, *llr_approx;
+
+  parse_args(argc, argv);
+
+  /* initialize objects */
+  if (modem_table_lte(&mod, modulation, true)) {
+    fprintf(stderr, "Error initializing modem table\n");
+    exit(-1);
+  }
+
+  /* check that num_bits is multiple of num_bits x symbol */
+  num_bits = mod.nbits_x_symbol * (num_bits / mod.nbits_x_symbol);
+
+  demod_soft_init(&demod_soft);
+  demod_soft_table_set(&demod_soft, &mod);
+  demod_soft_sigma_set(&demod_soft, 2.0 / mod.nbits_x_symbol);
+
+  /* allocate buffers */
+  input = malloc(sizeof(char) * num_bits);
+  if (!input) {
+    perror("malloc");
+    exit(-1);
+  }
+  output = malloc(sizeof(char) * num_bits);
+  if (!output) {
+    perror("malloc");
+    exit(-1);
+  }
+  symbols = malloc(sizeof(cf_t) * num_bits / mod.nbits_x_symbol);
+  if (!symbols) {
+    perror("malloc");
+    exit(-1);
+  }
+
+  llr_exact = malloc(sizeof(float) * num_bits);
+  if (!llr_exact) {
+    perror("malloc");
+    exit(-1);
+  }
+
+  llr_approx = malloc(sizeof(float) * num_bits);
+  if (!llr_approx) {
+    perror("malloc");
+    exit(-1);
+  }
+
+  /* generate random data */
+  srand(time(NULL));
+  
+  int ret = -1;
+  double mse;
+  struct timeval t[3]; 
+  float mean_texec = 0.0; 
+  for (int n=0;n<nof_frames;n++) {
+    for (i=0;i<num_bits;i++) {
+      input[i] = rand()%2;
+    }
+
+    /* modulate */
+    mod_modulate(&mod, input, symbols, num_bits);
+
+    /* add noise */
+    ch_awgn_c(symbols, symbols, ch_awgn_get_variance(5.0, mod.nbits_x_symbol), num_bits / mod.nbits_x_symbol);
+    
+    /* Compare exact with approximation algorithms */
+    demod_soft_alg_set(&demod_soft, EXACT);
+    demod_soft_demodulate(&demod_soft, symbols, llr_exact, num_bits / mod.nbits_x_symbol);
+    
+    demod_soft_alg_set(&demod_soft, APPROX);
+    gettimeofday(&t[1], NULL);
+    demod_soft_demodulate(&demod_soft, symbols, llr_approx, num_bits / mod.nbits_x_symbol);
+    gettimeofday(&t[2], NULL);
+    get_time_interval(t);
+    
+    /* compute exponentially averaged execution time */
+    mean_texec = EXPAVERAGE((float) t[0].tv_usec, mean_texec, n);
+    
+    /* check MSE */
+    mse = 0.0;
+    for (i=0;i<num_bits;i++) {
+      float e = llr_exact[i] - llr_approx[i];
+      mse += e*e;
+    }
+    mse/=num_bits;
+
+    if (VERBOSE_ISDEBUG()) {
+      printf("exact=");
+      vec_fprint_f(stdout, llr_exact, num_bits);
+
+      printf("approx=");
+      vec_fprint_f(stdout, llr_approx, num_bits);
+    }
+    
+    if (mse > mse_threshold()) {
+        goto clean_exit; 
+    }
+  }
+  ret = 0; 
+
+clean_exit:  
+  free(llr_exact);
+  free(llr_approx);
+  free(symbols);
+  free(output);
+  free(input);
+
+  modem_table_free(&mod);
+
+  if (ret == 0) {
+    printf("Ok Mean Throughput: %.2f. Mbps ExTime: %.2f us\n", num_bits/mean_texec, mean_texec);    
+  } else {
+    printf("Error: MSE too large (%f > %f)\n", mse, mse_threshold());
+  }
+  exit(ret);
+}

lte/phy/lib/phch/src/ue_sync.c

@@ -40,8 +40,6 @@
 #define MAX_TIME_OFFSET 128
 cf_t dummy[MAX_TIME_OFFSET];
 
-#define EXPAVERAGE(data, average, nframes) ((data + average * nframes) / (nframes + 1))  
-
 #define CURRENT_FFTSIZE   lte_symbol_sz(q->cell.nof_prb)
 #define CURRENT_SFLEN     SF_LEN(CURRENT_FFTSIZE, q->cell.cp)