Added DL channel estimator algorithm options

lib/examples/pdsch_ue.c

@@ -91,7 +91,7 @@ typedef struct {
   uint32_t file_nof_ports;
   uint32_t file_cell_id;
   bool     enable_cfo_ref;
-  bool     average_subframe;
+  char*    estimator_alg;
   char*    rf_dev;
   char*    rf_args;
   uint32_t rf_nof_rx_ant;
@@ -133,7 +133,7 @@ void args_default(prog_args_t* args)
   args->rf_freq                            = -1.0;
   args->rf_nof_rx_ant                      = 1;
   args->enable_cfo_ref                     = false;
-  args->average_subframe                   = false;
+  args->estimator_alg                      = "interpolate";
   args->enable_256qam                      = false;
 #ifdef ENABLE_AGC_DEFAULT
   args->rf_gain = -1.0;
@@ -176,7 +176,7 @@ void usage(prog_args_t* args, char* prog)
   printf("\t-l Force N_id_2 [Default best]\n");
   printf("\t-C Disable CFO correction [Default %s]\n", args->disable_cfo ? "Disabled" : "Enabled");
   printf("\t-F Enable RS-based CFO correction [Default %s]\n", !args->enable_cfo_ref ? "Disabled" : "Enabled");
-  printf("\t-R Average channel estimates on 1 ms [Default %s]\n", !args->average_subframe ? "Disabled" : "Enabled");
+  printf("\t-R Channel estimates algorithm (average, interpolate, wiener) [Default %s]\n", args->estimator_alg);
   printf("\t-t Add time offset [Default %d]\n", args->time_offset);
   printf("\t-T Set TDD special subframe configuration [Default %d]\n", args->tdd_special_sf);
   printf("\t-G Set TDD uplink/downlink configuration [Default %d]\n", args->sf_config);
@@ -242,7 +242,7 @@ void parse_args(prog_args_t* args, int argc, char** argv)
         args->enable_cfo_ref = true;
         break;
       case 'R':
-        args->average_subframe = true;
+        args->estimator_alg = argv[optind];
         break;
       case 't':
         args->time_offset = (uint32_t)strtol(argv[optind], NULL, 10);
@@ -606,13 +606,13 @@ int main(int argc, char** argv)
   srslte_chest_dl_cfg_t chest_pdsch_cfg = {};
   chest_pdsch_cfg.cfo_estimate_enable   = prog_args.enable_cfo_ref;
   chest_pdsch_cfg.cfo_estimate_sf_mask  = 1023;
-  chest_pdsch_cfg.interpolate_subframe  = !prog_args.average_subframe;
+  chest_pdsch_cfg.estimator_alg        = srslte_chest_dl_str2estimator_alg(prog_args.estimator_alg);
 
   // Special configuration for MBSFN channel estimation
   srslte_chest_dl_cfg_t chest_mbsfn_cfg = {};
   chest_mbsfn_cfg.filter_type           = SRSLTE_CHEST_FILTER_TRIANGLE;
   chest_mbsfn_cfg.filter_coef[0]        = 0.1;
-  chest_mbsfn_cfg.interpolate_subframe  = true;
+  chest_mbsfn_cfg.estimator_alg        = SRSLTE_ESTIMATOR_ALG_INTERPOLATE;
   chest_mbsfn_cfg.noise_alg             = SRSLTE_NOISE_ALG_PSS;
 
   // Allocate softbuffer buffers

lib/include/srslte/phy/ch_estimation/chest_dl.h

@@ -67,12 +67,20 @@ typedef struct SRSLTE_API {
   float    sync_error;
 } srslte_chest_dl_res_t;
 
+// Noise estimation algorithm
 typedef enum SRSLTE_API {
   SRSLTE_NOISE_ALG_REFS = 0,
   SRSLTE_NOISE_ALG_PSS,
   SRSLTE_NOISE_ALG_EMPTY,
 } srslte_chest_dl_noise_alg_t;
 
+// Channel estimator algorithm
+typedef enum SRSLTE_API {
+  SRSLTE_ESTIMATOR_ALG_AVERAGE = 0,
+  SRSLTE_ESTIMATOR_ALG_INTERPOLATE,
+  SRSLTE_ESTIMATOR_ALG_WIENER,
+} srslte_chest_dl_estimator_alg_t;
+
 typedef struct SRSLTE_API {
   srslte_cell_t cell;
   uint32_t      nof_rx_antennas;
@@ -114,12 +122,13 @@ typedef struct SRSLTE_API {
 
 typedef struct SRSLTE_API {
 
+  srslte_chest_dl_estimator_alg_t estimator_alg;
   srslte_chest_dl_noise_alg_t     noise_alg;
+
   srslte_chest_filter_t       filter_type;
   float                       filter_coef[2];
 
   uint16_t mbsfn_area_id;
-  bool     interpolate_subframe;
   bool     rsrp_neighbour;
   bool     cfo_estimate_enable;
   uint32_t cfo_estimate_sf_mask;
@@ -158,4 +167,6 @@ SRSLTE_API int srslte_chest_dl_estimate_cfg(srslte_chest_dl_t*     q,
                                             cf_t*                  input[SRSLTE_MAX_PORTS],
                                             srslte_chest_dl_res_t* res);
 
+SRSLTE_API srslte_chest_dl_estimator_alg_t srslte_chest_dl_str2estimator_alg(const char* str);
+
 #endif // SRSLTE_CHEST_DL_H

lib/src/phy/ch_estimation/chest_dl.c

@@ -446,7 +446,7 @@ static void interpolate_pilots(srslte_chest_dl_t*     q,
   /* Interpolate in the frequency domain */
 
   uint32_t freq_nsymbols = nsymbols;
-  if (!cfg->interpolate_subframe) {
+  if (cfg->estimator_alg == SRSLTE_ESTIMATOR_ALG_AVERAGE) {
     freq_nsymbols = 1;
   }
 
@@ -470,7 +470,7 @@ static void interpolate_pilots(srslte_chest_dl_t*     q,
                                     (fidx_offset) ? 1 : 2);
       }
     } else {
-      if (!cfg->interpolate_subframe && nsymbols > 1) {
+      if (cfg->estimator_alg == SRSLTE_ESTIMATOR_ALG_AVERAGE && nsymbols > 1) {
         fidx_offset = q->cell.id % 3;
         srslte_interp_linear_offset(
             &q->srslte_interp_lin_3, pilot_estimates, ce, fidx_offset, SRSLTE_NRE / 4 - fidx_offset);
@@ -487,7 +487,7 @@ static void interpolate_pilots(srslte_chest_dl_t*     q,
   }
 
   /* Now interpolate in the time domain between symbols */
-  if (sf->sf_type == SRSLTE_SF_NORM && (!cfg->interpolate_subframe || nsymbols < 3)) {
+  if (sf->sf_type == SRSLTE_SF_NORM && (cfg->estimator_alg == SRSLTE_ESTIMATOR_ALG_AVERAGE || nsymbols < 3)) {
     // If we average per subframe, just copy the estimates in the time domain
     for (uint32_t l = 1; l < 2 * SRSLTE_CP_NSYMB(q->cell.cp); l++) {
       memcpy(&ce[l * SRSLTE_NRE * q->cell.nof_prb], ce, sizeof(cf_t) * SRSLTE_NRE * q->cell.nof_prb);
@@ -547,7 +547,7 @@ static void average_pilots(srslte_chest_dl_t*     q,
   uint32_t nref = (sf->sf_type == SRSLTE_SF_MBSFN) ? 6 * q->cell.nof_prb : 2 * q->cell.nof_prb;
 
   // Average in the time domain if enabled
-  if (!cfg->interpolate_subframe && nsymbols > 1) {
+  if (cfg->estimator_alg == SRSLTE_ESTIMATOR_ALG_AVERAGE && nsymbols > 1) {
     cf_t* temp = output; // Use output as temporal buffer
 
     if (srslte_refsignal_cs_fidx(q->cell, 0, port_id, 0) < 3) {
@@ -645,24 +645,35 @@ static void chest_interpolate_noise_est(srslte_chest_dl_t*     q,
     q->noise_estimate[rxant_id][port_id] = estimate_noise_pilots(q, sf, port_id);
   }
 
-  if (q->wiener_dl && ch_mode == SRSLTE_SF_NORM) {
+  if (q->wiener_dl && ch_mode == SRSLTE_SF_NORM && cfg->estimator_alg == SRSLTE_ESTIMATOR_ALG_WIENER) {
+    bool     ready   = q->wiener_dl->ready;
     uint32_t nre   = q->cell.nof_prb * SRSLTE_NRE;
     uint32_t nref  = q->cell.nof_prb * 2;
     uint32_t shift = srslte_refsignal_cs_fidx(q->cell, 0, port_id, 0);
-    for (uint32_t m = 0, l = 0; m < 2 * SRSLTE_CP_NORM_NSYMB; m++) {
+    uint32_t nsymb   = srslte_refsignal_cs_nof_symbols(&q->csr_refs, sf, port_id);
-      uint32_t k = srslte_refsignal_cs_nsymbol(l, q->cell.cp, 0);
+    float    snr_lin = +INFINITY;
-      srslte_wiener_dl_run(q->wiener_dl,
-                           port_id,
-                           rxant_id,
-                           m,
-                           shift,
-                           &input[nref * l],
-                           &ce[nre * m],
-                           q->noise_estimate[rxant_id][port_id]);
 
-      if (l == k) {
+    if (isnormal(q->noise_estimate[rxant_id][port_id]) && isnormal(q->rsrp[rxant_id][port_id])) {
-        l++;
+      snr_lin = q->rsrp[rxant_id][port_id] / q->noise_estimate[rxant_id][port_id] / 2;
     }
+
+    for (uint32_t m = 0, l = 0; m < 2 * SRSLTE_CP_NORM_NSYMB + 4; m++) {
+      uint32_t ce_idx = 0;
+
+      if (m >= 4) {
+        ce_idx = (m - 4) * nre;
+      }
+
+      uint32_t k = srslte_refsignal_cs_nsymbol(l, q->cell.cp, port_id);
+      srslte_wiener_dl_run(
+          q->wiener_dl, port_id, rxant_id, m, shift, &q->pilot_estimates[nref * l], &ce[ce_idx], snr_lin);
+
+      if (m == k) {
+        l = (l + 1) % nsymb;
+      }
+    }
+    if (ready) {
+      return;
     }
   }
 
@@ -686,7 +697,7 @@ static void chest_interpolate_noise_est(srslte_chest_dl_t*     q,
         break;
     }
 
-    if (!cfg->interpolate_subframe && ch_mode == SRSLTE_SF_MBSFN) {
+    if (cfg->estimator_alg != SRSLTE_ESTIMATOR_ALG_INTERPOLATE && ch_mode == SRSLTE_SF_MBSFN) {
       ERROR("Warning: Subframe interpolation must be enabled in MBSFN subframes\n");
     }
 
@@ -970,3 +981,20 @@ int srslte_chest_dl_estimate_cfg(srslte_chest_dl_t*     q,
 
   return SRSLTE_SUCCESS;
 }
+
+srslte_chest_dl_estimator_alg_t srslte_chest_dl_str2estimator_alg(const char* str)
+{
+  srslte_chest_dl_estimator_alg_t ret = SRSLTE_ESTIMATOR_ALG_AVERAGE;
+
+  if (str) {
+    if (strcmp(str, "interpolate") == 0) {
+      ret = SRSLTE_ESTIMATOR_ALG_INTERPOLATE;
+    } else if (strcmp(str, "average") == 0) {
+      ret = SRSLTE_ESTIMATOR_ALG_AVERAGE;
+    } else if (strcmp(str, "wiener") == 0) {
+      ret = SRSLTE_ESTIMATOR_ALG_WIENER;
+    }
+  }
+
+  return ret;
+}

lib/src/phy/phch/test/pmch_file_test.c

@@ -195,7 +195,7 @@ int main(int argc, char** argv)
   // Special configuration for MBSFN channel estimation
   ue_dl_cfg.chest_cfg.filter_type          = SRSLTE_CHEST_FILTER_TRIANGLE;
   ue_dl_cfg.chest_cfg.filter_coef[0]       = 0.1;
-  ue_dl_cfg.chest_cfg.interpolate_subframe = true;
+  ue_dl_cfg.chest_cfg.estimator_alg        = SRSLTE_ESTIMATOR_ALG_INTERPOLATE;
   ue_dl_cfg.chest_cfg.noise_alg            = SRSLTE_NOISE_ALG_PSS;
 
   if ((ret = srslte_ue_dl_decode_fft_estimate(&ue_dl, &dl_sf, &ue_dl_cfg)) < 0) {

lib/src/phy/utils/mat.c

@@ -443,15 +443,23 @@ static inline void srslte_vec_sub_ccc_simd_inline(const cf_t* x, const cf_t* y,
 
 static inline cf_t reciprocal(cf_t x)
 {
-  cf_t y;
+  cf_t y = 0.0f;
 
   float mod  = __real__ x * __real__ x + __imag__ x * __imag__ x;
+
+  if (isnormal(mod)) {
     __real__ y = (__real__ x) / mod;
     __imag__ y = -(__imag__ x) / mod;
+  }
 
   return y;
 }
 
+static inline float _cabs2(cf_t x)
+{
+  return __real__ x * __real__ x + __imag__ x * __imag__ x;
+}
+
 void srslte_matrix_NxN_inv_run(srslte_matrix_NxN_inv_t* q, cf_t* in, cf_t* out)
 {
   if (q && in && out) {
@@ -476,8 +484,29 @@ void srslte_matrix_NxN_inv_run(srslte_matrix_NxN_inv_t* q, cf_t* in, cf_t* out)
 
     // 1) Forward elimination
     for (int i = 0; i < N - 1; i++) {
-      tmp_src_ptr = &q->matrix[N * 2 * (N - 1 - i)];
+      uint32_t row_i = N - i - 1;
-      cf_t b      = tmp_src_ptr[N - 1 - i];
+      uint32_t col_i = N - i - 1;
+      float    max_v = 0.0f;
+      uint32_t max_i = 0;
+
+      // Find Row with highest value in the column to reduce
+      for (int j = 0; j < N - i; j++) {
+        float v = _cabs2(q->matrix[(j + 1) * 2 * N - 1 - i]);
+        if (v > max_v) {
+          max_i = j;
+          max_v = v;
+        }
+      }
+
+      // Swap rows
+      if (max_i != row_i) {
+        memcpy(q->row_buffer, &q->matrix[row_i * N * 2], sizeof(cf_t) * N * 2);
+        memcpy(&q->matrix[row_i * N * 2], &q->matrix[max_i * N * 2], sizeof(cf_t) * N * 2);
+        memcpy(&q->matrix[max_i * N * 2], q->row_buffer, sizeof(cf_t) * N * 2);
+      }
+
+      tmp_src_ptr = &q->matrix[N * 2 * row_i]; // Select row
+      cf_t b      = tmp_src_ptr[col_i];
       srslte_vec_sc_prod_ccc_simd_inline(tmp_src_ptr, reciprocal(b), tmp_src_ptr, 2 * N);
 
       for (int j = 0; j < N - i - 1; j++) {

lib/test/phy/phy_dl_test.c

@@ -578,7 +578,7 @@ int main(int argc, char** argv)
     ue_dl_cfg.chest_cfg.filter_type          = SRSLTE_CHEST_FILTER_GAUSS;
     ue_dl_cfg.chest_cfg.noise_alg            = SRSLTE_NOISE_ALG_REFS;
     ue_dl_cfg.chest_cfg.rsrp_neighbour       = false;
-    ue_dl_cfg.chest_cfg.interpolate_subframe = false;
+    ue_dl_cfg.chest_cfg.estimator_alg        = SRSLTE_ESTIMATOR_ALG_AVERAGE;
     ue_dl_cfg.chest_cfg.cfo_estimate_enable  = false;
     ue_dl_cfg.chest_cfg.cfo_estimate_sf_mask = false;
     ue_dl_cfg.chest_cfg.sync_error_enable    = false;

srsue/src/phy/cc_worker.cc

@@ -87,7 +87,7 @@ cc_worker::cc_worker(uint32_t cc_idx_, uint32_t max_prb, srsue::phy_common* phy_
   // Define MBSFN subframes channel estimation and save default one
   chest_mbsfn_cfg.filter_type          = SRSLTE_CHEST_FILTER_TRIANGLE;
   chest_mbsfn_cfg.filter_coef[0]       = 0.1;
-  chest_mbsfn_cfg.interpolate_subframe = true;
+  chest_mbsfn_cfg.estimator_alg        = SRSLTE_ESTIMATOR_ALG_INTERPOLATE;
   chest_mbsfn_cfg.noise_alg            = SRSLTE_NOISE_ALG_PSS;
 
   chest_default_cfg = ue_dl_cfg.chest_cfg;
@@ -144,8 +144,8 @@ bool cc_worker::set_cell(srslte_cell_t cell_)
       return false;
     }
 
-    if (cell.frame_type == SRSLTE_TDD && !ue_dl_cfg.chest_cfg.interpolate_subframe) {
+    if (cell.frame_type == SRSLTE_TDD && ue_dl_cfg.chest_cfg.estimator_alg != SRSLTE_ESTIMATOR_ALG_INTERPOLATE) {
-      chest_default_cfg.interpolate_subframe = true;
+      chest_default_cfg.estimator_alg = SRSLTE_ESTIMATOR_ALG_INTERPOLATE;
       log_h->console("Enabling subframe interpolation for TDD cells (recommended setting)\n");
     }
 

srsue/src/phy/phy_common.cc

@@ -103,7 +103,8 @@ void phy_common::set_ue_dl_cfg(srslte_ue_dl_cfg_t* ue_dl_cfg)
 
   chest_cfg->rsrp_neighbour       = false;
   chest_cfg->sync_error_enable    = args->correct_sync_error;
-  chest_cfg->interpolate_subframe = args->interpolate_subframe_enabled;
+  chest_cfg->estimator_alg =
+      args->interpolate_subframe_enabled ? SRSLTE_ESTIMATOR_ALG_INTERPOLATE : SRSLTE_ESTIMATOR_ALG_AVERAGE;
   chest_cfg->cfo_estimate_enable  = args->cfo_ref_mask != 0;
   chest_cfg->cfo_estimate_sf_mask = args->cfo_ref_mask;
 }