From ecf140a4e9f25469c8ef533dc2ba6984f4c586cd Mon Sep 17 00:00:00 2001
From: ismagom <gomezi@tcd.ie>
Date: Tue, 22 Jul 2014 18:13:24 +0200
Subject: [PATCH] Added Test and Benchmark for soft_demodulation approximation
 algorithms

---
 lte/phy/include/liblte/phy/channel/ch_awgn.h |   8 +-
 lte/phy/include/liblte/phy/utils/vector.h    |   3 +
 lte/phy/lib/channel/src/ch_awgn.c            |  20 +-
 lte/phy/lib/modem/src/demod_soft.c           |   1 +
 lte/phy/lib/modem/src/soft_algs.c            | 112 ++++-----
 lte/phy/lib/modem/test/CMakeLists.txt        |   8 +
 lte/phy/lib/modem/test/soft_demod_test.c     | 229 +++++++++++++++++++
 lte/phy/lib/phch/src/ue_sync.c               |   2 -
 8 files changed, 318 insertions(+), 65 deletions(-)
 create mode 100644 lte/phy/lib/modem/test/soft_demod_test.c

diff --git a/lte/phy/include/liblte/phy/channel/ch_awgn.h b/lte/phy/include/liblte/phy/channel/ch_awgn.h
index fa209dfe3..940ad332d 100644
--- a/lte/phy/include/liblte/phy/channel/ch_awgn.h
+++ b/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 */
 
diff --git a/lte/phy/include/liblte/phy/utils/vector.h b/lte/phy/include/liblte/phy/utils/vector.h
index 369428fb9..0e6b5a041 100644
--- a/lte/phy/include/liblte/phy/utils/vector.h
+++ b/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);
diff --git a/lte/phy/lib/channel/src/ch_awgn.c b/lte/phy/lib/channel/src/ch_awgn.c
index 54a70b482..05dadf975 100644
--- a/lte/phy/lib/channel/src/ch_awgn.c
+++ b/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) {
-  _Complex float tmp;
-  int i;
+float ch_awgn_get_variance(float ebno_db, float rate) {
+  float esno_db = ebno_db + 10 * log10f(rate);
+  return sqrtf(1 / (powf(10, esno_db / 10)));
+}
 
-  for (i=0;i<buff_sz;i++) {
+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<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) {
-  int i;
+void ch_awgn_f(const float* x, float* y, float variance, uint32_t len) {
+  uint32_t i;
 
-  for (i=0;i<buff_sz;i++) {
+  for (i=0;i<len;i++) {
     y[i] = x[i] + variance * rand_gauss();
   }
 }
diff --git a/lte/phy/lib/modem/src/demod_soft.c b/lte/phy/lib/modem/src/demod_soft.c
index f79863199..d030a3733 100644
--- a/lte/phy/lib/modem/src/demod_soft.c
+++ b/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) {
diff --git a/lte/phy/lib/modem/src/soft_algs.c b/lte/phy/lib/modem/src/soft_algs.c
index 12e12591c..cb8d4e40b 100644
--- a/lte/phy/lib/modem/src/soft_algs.c
+++ b/lte/phy/lib/modem/src/soft_algs.c
@@ -50,42 +50,44 @@
  * \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,
-    _Complex float *symbols, uint32_t (*S)[6][32], float sigma2) {
+void
+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;
   float x, y, d[64];
 
-  for (s=0; s<N; s++) {  	/* recevied symbols */
-	  /* Compute the distances squared d[i] between the received symbol and all constellation points */
-	  for (i=0; i<M; i++) {
-		x = __real__ in[s] - __real__ symbols[i];
-		y = __imag__ in[s] - __imag__ symbols[i];
-		d[i] = x*x + y*y;
-	  }
+  for (s = 0; s < N; s++) {     /* recevied symbols */
+    /* Compute the distances squared d[i] between the received symbol and all constellation points */
+    for (i = 0; i < M; i++) {
+      x = __real__ in[s] - __real__ symbols[i];
+      y = __imag__ in[s] - __imag__ symbols[i];
+      d[i] = x * x + y * y;
+    }
 
-	  for (b=0; b<B; b++) {/* bits per symbol*/
-	  	/* initiate num[b] and den[b] */
-		num = d[S[0][b][0]];
-		den = d[S[1][b][0]];
+    for (b = 0; b < B; b++) {   /* bits per symbol */
+      /* initiate num[b] and den[b] */
+      num = d[S[0][b][0]];
+      den = d[S[1][b][0]];
 
-		/* Minimum distance squared search between recevied symbol and a constellation point with a
-		 '1' and a '0' for each bit position */
-	        for (i=1; i<M/2; i++) { /* half the constellation points have '1'|'0' at any given bit position */
-	        	if (d[S[0][b][i]] < num) {
-	        		num = d[S[0][b][i]];
-	        	}
-	        	if (d[S[1][b][i]] < den) {
-	        		den = d[S[1][b][i]];
-	        	}
-	        }
-	        /* Theoretical LLR and approximate LLR values are positive if
-	         * symbol(s) with '0' is/are closer and negative if symbol(s)
-	         * with '1' are closer.
-	         * Change sign if mapping negative to '0' and positive to '1' */
-	        out[s*B+b] = change_sign*(den-num)/sigma2;
-	  }
+      /* Minimum distance squared search between recevied symbol and a constellation point with a
+         '1' and a '0' for each bit position */
+      for (i = 1; i < M / 2; i++) {     /* half the constellation points have '1'|'0' at any given bit position */
+        if (d[S[0][b][i]] < num) {
+          num = d[S[0][b][i]];
+        }
+        if (d[S[1][b][i]] < den) {
+          den = d[S[1][b][i]];
+        }
+      }
+      /* Theoretical LLR and approximate LLR values are positive if
+       * symbol(s) with '0' is/are closer and negative if symbol(s)
+       * with '1' are closer.
+       * Change sign if mapping negative to '0' and positive to '1' */
+      out[s * B + b] = change_sign * (den - num) / sigma2;
+    }
   }
 
 }
@@ -106,36 +108,38 @@ 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,
-    _Complex float *symbols, uint32_t (*S)[6][32], float sigma2) {
+void
+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;
   float x, y, d[64];
 
-  for (s=0; s<N; s++) {  	/* recevied symbols */
-    	/* Compute exp{·} of the distances squared d[i] between the received symbol and all constellation points */
-    	for (i=0; i<M; i++) {
-    		x = __real__ in[s] - __real__ symbols[i];
-    		y = __imag__ in[s] - __imag__ symbols[i];
-    		d[i] = exp(-1*(x*x + y*y)/sigma2);
-    	}
+  for (s = 0; s < N; s++) {     /* recevied symbols */
+    /* Compute exp{·} of the distances squared d[i] between the received symbol and all constellation points */
+    for (i = 0; i < M; i++) {
+      x = __real__ in[s] - __real__ symbols[i];
+      y = __imag__ in[s] - __imag__ symbols[i];
+      d[i] = exp(-1 * (x * x + y * y) / sigma2);
+    }
 
-	/* Sum up the corresponding d[i]'s for each bit position */
-    	for (b=0; b<B; b++) {/* bits per symbol*/
-    		/* initiate num[b] and den[b] */
-    		num = 0;
-    		den = 0;
-    	    	
-    	        for (i=0; i<M/2; i++) { /* half the constellation points have '1'|'0' at any given bit position */
-    	        	num += d[S[0][b][i]];
-    	        	den += d[S[1][b][i]];
-    	        }
-    	        /* Theoretical LLR and approximate LLR values are positive if
-    	         * symbol(s) with '0' is/are closer and negative if symbol(s)
-    	         * with '1' are closer.
-    	         * Change sign if mapping negative to '0' and positive to '1' */
-    	        out[s*B+b] = change_sign*log(num/den);
-	}
+    /* Sum up the corresponding d[i]'s for each bit position */
+    for (b = 0; b < B; b++) {   /* bits per symbol */
+      /* initiate num[b] and den[b] */
+      num = 0;
+      den = 0;
+
+      for (i = 0; i < M / 2; i++) {     /* half the constellation points have '1'|'0' at any given bit position */
+        num += d[S[0][b][i]];
+        den += d[S[1][b][i]];
+      }
+      /* Theoretical LLR and approximate LLR values are positive if
+       * symbol(s) with '0' is/are closer and negative if symbol(s)
+       * with '1' are closer.
+       * Change sign if mapping negative to '0' and positive to '1' */
+      out[s * B + b] = change_sign * log(num / den);
+    }
   }
 }
diff --git a/lte/phy/lib/modem/test/CMakeLists.txt b/lte/phy/lib/modem/test/CMakeLists.txt
index 70308be83..24ebc1986 100644
--- a/lte/phy/lib/modem/test/CMakeLists.txt
+++ b/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)
+ 
 
 
diff --git a/lte/phy/lib/modem/test/soft_demod_test.c b/lte/phy/lib/modem/test/soft_demod_test.c
new file mode 100644
index 000000000..a5621e305
--- /dev/null
+++ b/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);
+}
diff --git a/lte/phy/lib/phch/src/ue_sync.c b/lte/phy/lib/phch/src/ue_sync.c
index b877cf3b8..c62739da9 100644
--- a/lte/phy/lib/phch/src/ue_sync.c
+++ b/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)