Noise power estimation and MMSE computation for DL at UE side

openair1/PHY/NR_UE_ESTIMATION/nr_dl_channel_estimation.c

@@ -1303,7 +1303,8 @@ void NFAPI_NR_DMRS_TYPE1_linear_interp(NR_DL_FRAME_PARMS *frame_parms,
                                        unsigned short bwp_start_subcarrier,
                                        unsigned short nb_rb_pdsch,
                                        int8_t delta,
-                                       delay_t *delay)
+                                       delay_t *delay,
+                                       uint32_t *nvar)
 {
   c16_t *dl_ch0 = dl_ch;
   int re_offset = bwp_start_subcarrier % frame_parms->ofdm_symbol_size;
@@ -1311,6 +1312,9 @@ void NFAPI_NR_DMRS_TYPE1_linear_interp(NR_DL_FRAME_PARMS *frame_parms,
   c16_t dl_ls_est[frame_parms->ofdm_symbol_size] __attribute__((aligned(32)));
   memset(dl_ls_est, 0, sizeof(dl_ls_est));
 
+  int nest_count = 0;
+  uint64_t noise_amp2 = 0;
+
   for (int pilot_cnt = 0; pilot_cnt < 6 * nb_rb_pdsch; pilot_cnt++) {
     if (pilot_cnt % 2 == 0) {
       c16_t ch = c16mulShift(*pil, rxF[re_offset], 15);
@@ -1366,6 +1370,12 @@ void NFAPI_NR_DMRS_TYPE1_linear_interp(NR_DL_FRAME_PARMS *frame_parms,
   c16_t *dl_inv_delay_table = frame_parms->delay_table[inv_delay_idx];
   for (int k = 0; k < 12 * nb_rb_pdsch; k++) {
     dl_ch[k] = c16mulShift(dl_ch[k], dl_inv_delay_table[k], 8);
+    noise_amp2 += c16amp2(c16sub(dl_ls_est[k], dl_ch[k]));
+    nest_count++;
+  }
+
+  if (nvar && nest_count > 0) {
+    *nvar = (uint32_t)(noise_amp2 / (nest_count * frame_parms->nb_antennas_rx));
   }
 }
 
@@ -1462,7 +1472,8 @@ void NFAPI_NR_DMRS_TYPE2_linear_interp(NR_DL_FRAME_PARMS *frame_parms,
                                        unsigned short nb_rb_pdsch,
                                        int8_t delta,
                                        unsigned short p,
-                                       delay_t *delay)
+                                       delay_t *delay,
+                                       uint32_t *nvar)
 {
   int re_offset = bwp_start_subcarrier % frame_parms->ofdm_symbol_size;
   c16_t *dl_ch0 = dl_ch;
@@ -1470,6 +1481,9 @@ void NFAPI_NR_DMRS_TYPE2_linear_interp(NR_DL_FRAME_PARMS *frame_parms,
   c16_t dl_ls_est[frame_parms->ofdm_symbol_size] __attribute__((aligned(32)));
   memset(dl_ls_est, 0, sizeof(dl_ls_est));
 
+  int nest_count = 0;
+  uint64_t noise_amp2 = 0;
+
   for (int pilot_cnt = 0; pilot_cnt < 4 * nb_rb_pdsch; pilot_cnt += 2) {
     c16_t ch_l = c16mulShift(*pil, rxF[re_offset], 15);
 #ifdef DEBUG_PDSCH
@@ -1518,6 +1532,12 @@ void NFAPI_NR_DMRS_TYPE2_linear_interp(NR_DL_FRAME_PARMS *frame_parms,
   c16_t *dl_inv_delay_table = frame_parms->delay_table[inv_delay_idx];
   for (int k = 0; k < 12 * nb_rb_pdsch; k++) {
     dl_ch[k] = c16mulShift(dl_ch[k], dl_inv_delay_table[k], 8);
+    noise_amp2 += c16amp2(c16sub(dl_ls_est[k], dl_ch[k]));
+    nest_count++;
+  }
+
+  if (nvar && nest_count > 0) {
+    *nvar = (uint32_t)(noise_amp2 / (nest_count * frame_parms->nb_antennas_rx));
   }
 }
 
@@ -1619,7 +1639,8 @@ int nr_pdsch_channel_estimation(PHY_VARS_NR_UE *ue,
                                 uint32_t pdsch_est_size,
                                 int32_t dl_ch_estimates[][pdsch_est_size],
                                 int rxdataFsize,
-                                c16_t rxdataF[][rxdataFsize])
+                                c16_t rxdataF[][rxdataFsize],
+                                uint32_t *nvar)
 {
   int gNB_id = proc->gNB_id;
   int Ns = proc->nr_slot_rx;
@@ -1679,7 +1700,8 @@ int nr_pdsch_channel_estimation(PHY_VARS_NR_UE *ue,
                                         bwp_start_subcarrier,
                                         nb_rb_pdsch,
                                         delta,
-                                        &delay);
+                                        &delay,
+                                        nvar);
 
     } else if (config_type == NFAPI_NR_DMRS_TYPE2 && ue->chest_freq == 0) {
       NFAPI_NR_DMRS_TYPE2_linear_interp(&ue->frame_parms,
@@ -1690,7 +1712,8 @@ int nr_pdsch_channel_estimation(PHY_VARS_NR_UE *ue,
                                         nb_rb_pdsch,
                                         delta,
                                         p,
-                                        &delay);
+                                        &delay,
+                                        nvar);
 
     } else if (config_type == NFAPI_NR_DMRS_TYPE1) {
       NFAPI_NR_DMRS_TYPE1_average_prb(&ue->frame_parms,

openair1/PHY/NR_UE_ESTIMATION/nr_estimation.h

@@ -89,7 +89,8 @@ int nr_pdsch_channel_estimation(PHY_VARS_NR_UE *ue,
                                 uint32_t pdsch_est_size,
                                 int32_t dl_ch_estimates[][pdsch_est_size],
                                 int rxdataFsize,
-                                c16_t rxdataF[][rxdataFsize]);
+                                c16_t rxdataF[][rxdataFsize],
+                                uint32_t *nvar);
 
 void nr_adjust_synch_ue(NR_DL_FRAME_PARMS *frame_parms,
                         PHY_VARS_NR_UE *ue,

openair1/PHY/NR_UE_TRANSPORT/nr_dlsch_demodulation.c

@@ -81,20 +81,21 @@ unsigned char offset_mumimo_llr_drange[29][3]={{8,8,8},{7,7,7},{7,7,7},{7,7,7},{
 #define print_ints(s,x) printf("%s = %d %d %d %d\n",s,(x)[0],(x)[1],(x)[2],(x)[3])
 #define print_shorts(s,x) printf("%s = [%d+j*%d, %d+j*%d, %d+j*%d, %d+j*%d]\n",s,(x)[0],(x)[1],(x)[2],(x)[3],(x)[4],(x)[5],(x)[6],(x)[7])
 
-/* compute H_h_H matrix inversion up to 4x4 matrices */
-uint8_t nr_zero_forcing_rx(uint32_t rx_size_symbol,
-                           unsigned char n_rx,
-                           unsigned char n_tx,//number of layer
-                           int32_t rxdataF_comp[][n_rx][rx_size_symbol * NR_SYMBOLS_PER_SLOT],
-                           int32_t dl_ch_mag[][n_rx][rx_size_symbol],
-                           int32_t dl_ch_magb[][n_rx][rx_size_symbol],
-                           int32_t dl_ch_magr[][n_rx][rx_size_symbol],
-                           int32_t dl_ch_estimates_ext[][rx_size_symbol],
-                           unsigned short nb_rb,
-                           unsigned char mod_order,
-                           int shift,
-                           unsigned char symbol,
-                           int length);
+/* compute the MMSE up to 4x4 matrices */
+uint8_t nr_dlsch_mmse(uint32_t rx_size_symbol,
+                      unsigned char n_rx,
+                      unsigned char n_tx, // number of layer
+                      int32_t rxdataF_comp[][n_rx][rx_size_symbol * NR_SYMBOLS_PER_SLOT],
+                      int32_t dl_ch_mag[][n_rx][rx_size_symbol],
+                      int32_t dl_ch_magb[][n_rx][rx_size_symbol],
+                      int32_t dl_ch_magr[][n_rx][rx_size_symbol],
+                      int32_t dl_ch_estimates_ext[][rx_size_symbol],
+                      unsigned short nb_rb,
+                      unsigned char mod_order,
+                      int shift,
+                      unsigned char symbol,
+                      int length,
+                      uint32_t noise_var);
 
 /* Apply layer demapping */
 static void nr_dlsch_layer_demapping(int16_t *llr_cw[2],
@@ -249,7 +250,8 @@ int nr_rx_pdsch(PHY_VARS_NR_UE *ue,
                 int nbRx,
                 int32_t rxdataF_comp[][nbRx][rx_size_symbol * NR_SYMBOLS_PER_SLOT],
                 c16_t ptrs_phase_per_slot[][NR_SYMBOLS_PER_SLOT],
-                int32_t ptrs_re_per_slot[][NR_SYMBOLS_PER_SLOT])
+                int32_t ptrs_re_per_slot[][NR_SYMBOLS_PER_SLOT],
+                uint32_t nvar)
 {
   const int nl = dlsch[0].Nl;
   const int matrixSz = ue->frame_parms.nb_antennas_rx * nl;
@@ -528,20 +530,23 @@ int nr_rx_pdsch(PHY_VARS_NR_UE *ue,
                            symbol,
                            nb_rb_pdsch,
                            nb_re_pdsch);
-    if (nl >= 2) // Apply zero forcing for 2, 3, and 4 Tx layers
-      nr_zero_forcing_rx(rx_size_symbol,
-                         n_rx,
-                         nl,
-                         rxdataF_comp,
-                         dl_ch_mag,
-                         dl_ch_magb,
-                         dl_ch_magr,
-                         dl_ch_estimates_ext,
-                         nb_rb_pdsch,
-                         dlsch[0].dlsch_config.qamModOrder,
-                         *log2_maxh,
-                         symbol,
-                         nb_re_pdsch);
+    // Apply zero forcing for 2, 3, and 4 Tx layers
+    if (nl >= 2) {
+      nr_dlsch_mmse(rx_size_symbol,
+                    n_rx,
+                    nl,
+                    rxdataF_comp,
+                    dl_ch_mag,
+                    dl_ch_magb,
+                    dl_ch_magr,
+                    dl_ch_estimates_ext,
+                    nb_rb_pdsch,
+                    dlsch[0].dlsch_config.qamModOrder,
+                    *log2_maxh,
+                    symbol,
+                    nb_re_pdsch,
+                    nvar);
+    }
   }
 
   if (meas_enabled) {
@@ -1668,23 +1673,23 @@ void nr_conjch0_mult_ch1(int *ch0,
   simde_m_empty();
 }
 
-/* Zero Forcing Rx function: up to 4 layers
- *
- *
- * */
-uint8_t nr_zero_forcing_rx(uint32_t rx_size_symbol,
-                           unsigned char n_rx,
-                           unsigned char n_tx,//number of layer
-                           int32_t rxdataF_comp[][n_rx][rx_size_symbol * NR_SYMBOLS_PER_SLOT],
-                           int32_t dl_ch_mag[][n_rx][rx_size_symbol],
-                           int32_t dl_ch_magb[][n_rx][rx_size_symbol],
-                           int32_t dl_ch_magr[][n_rx][rx_size_symbol],
-                           int32_t dl_ch_estimates_ext[][rx_size_symbol],
-                           unsigned short nb_rb,
-                           unsigned char mod_order,
-                           int shift,
-                           unsigned char symbol,
-                           int length)
+/*
+ * MMSE Rx function: up to 4 layers
+ */
+uint8_t nr_dlsch_mmse(uint32_t rx_size_symbol,
+                      unsigned char n_rx,
+                      unsigned char n_tx, // number of layer
+                      int32_t rxdataF_comp[][n_rx][rx_size_symbol * NR_SYMBOLS_PER_SLOT],
+                      int32_t dl_ch_mag[][n_rx][rx_size_symbol],
+                      int32_t dl_ch_magb[][n_rx][rx_size_symbol],
+                      int32_t dl_ch_magr[][n_rx][rx_size_symbol],
+                      int32_t dl_ch_estimates_ext[][rx_size_symbol],
+                      unsigned short nb_rb,
+                      unsigned char mod_order,
+                      int shift,
+                      unsigned char symbol,
+                      int length,
+                      uint32_t noise_var)
 {
   int *ch0r, *ch0c;
   uint32_t nb_rb_0 = length/12 + ((length%12)?1:0);
@@ -1716,6 +1721,18 @@ uint8_t nr_zero_forcing_rx(uint32_t rx_size_symbol,
     }
   }
 
+  // Add noise_var such that: H^h * H + noise_var * I
+  if (noise_var != 0) {
+    __m128i nvar_128i = simde_mm_set1_epi32(noise_var >> 3);
+    for (int p = 0; p < n_tx; p++) {
+      __m128i *conjH_H_128i = (__m128i *)conjH_H_elements[0][p][p];
+      for (int k = 0; k < 3 * nb_rb_0; k++) {
+        conjH_H_128i[0] = simde_mm_add_epi32(conjH_H_128i[0], nvar_128i);
+        conjH_H_128i++;
+      }
+    }
+  }
+
   //Compute the inverse and determinant of the H^*H matrix
   //Allocate the inverse matrix
   c16_t *inv_H_h_H[n_tx][n_tx];

openair1/PHY/NR_UE_TRANSPORT/nr_transport_proto_ue.h

@@ -412,12 +412,13 @@ int nr_rx_pdsch(PHY_VARS_NR_UE *ue,
                 uint32_t dl_valid_re[NR_SYMBOLS_PER_SLOT],
                 c16_t rxdataF[][ue->frame_parms.samples_per_slot_wCP],
                 uint32_t llr_offset[NR_SYMBOLS_PER_SLOT],
-                int32_t *log2_maxhrx_size_symbol,
+                int32_t *log2_maxh,
                 int rx_size_symbol,
                 int nbRx,
                 int32_t rxdataF_comp[][nbRx][rx_size_symbol * NR_SYMBOLS_PER_SLOT],
                 c16_t ptrs_phase_per_slot[][NR_SYMBOLS_PER_SLOT],
-                int32_t ptrs_re_per_slot[][NR_SYMBOLS_PER_SLOT]);
+                int32_t ptrs_re_per_slot[][NR_SYMBOLS_PER_SLOT],
+                uint32_t nvar);
 
 int32_t generate_nr_prach(PHY_VARS_NR_UE *ue, uint8_t gNB_id, int frame, uint8_t slot);
 

openair1/SCHED_NR_UE/phy_procedures_nr_ue.c

@@ -530,24 +530,30 @@ static int nr_ue_pdsch_procedures(PHY_VARS_NR_UE *ue,
     __attribute__((aligned(32))) int32_t rxdataF_comp[dlsch[0].Nl][ue->frame_parms.nb_antennas_rx][rx_size_symbol * NR_SYMBOLS_PER_SLOT];
     memset(rxdataF_comp, 0, sizeof(rxdataF_comp));
 
+    uint32_t nvar = 0;
+
     for (m = s0; m < (s0 +s1); m++) {
       if (dlsch0->dlsch_config.dlDmrsSymbPos & (1 << m)) {
         for (uint8_t aatx=0; aatx<dlsch0->Nl; aatx++) {//for MIMO Config: it shall loop over no_layers
           LOG_D(PHY,"PDSCH Channel estimation gNB id %d, PDSCH antenna port %d, slot %d, symbol %d\n",0,aatx,nr_slot_rx,m);
+          uint32_t nvar_tmp = 0;
           nr_pdsch_channel_estimation(ue,
                                       proc,
-                                      get_dmrs_port(aatx,dlsch0->dlsch_config.dmrs_ports),
+                                      get_dmrs_port(aatx, dlsch0->dlsch_config.dmrs_ports),
                                       m,
                                       dlsch0->dlsch_config.nscid,
                                       dlsch0->dlsch_config.dlDmrsScramblingId,
                                       BWPStart,
                                       dlsch0->dlsch_config.dmrsConfigType,
                                       dlsch0->dlsch_config.rb_offset,
-                                      ue->frame_parms.first_carrier_offset+(BWPStart + pdsch_start_rb)*12,
+                                      ue->frame_parms.first_carrier_offset + (BWPStart + pdsch_start_rb) * 12,
                                       pdsch_nb_rb,
                                       pdsch_est_size,
                                       pdsch_dl_ch_estimates,
-                                      ue->frame_parms.samples_per_slot_wCP, rxdataF);
+                                      ue->frame_parms.samples_per_slot_wCP,
+                                      rxdataF,
+                                      &nvar_tmp);
+          nvar += nvar_tmp;
 #if 0
           ///LOG_M: the channel estimation
           int nr_frame_rx = proc->frame_rx;
@@ -561,6 +567,9 @@ static int nr_ue_pdsch_procedures(PHY_VARS_NR_UE *ue,
         }
       }
     }
+
+    nvar /= (s1 * dlsch0->Nl * ue->frame_parms.nb_antennas_rx);
+
     nr_ue_measurement_procedures(2, ue, proc, &dlsch[0], pdsch_est_size, pdsch_dl_ch_estimates);
 
     if (ue->chest_time == 1) { // averaging time domain channel estimates
@@ -619,7 +628,8 @@ static int nr_ue_pdsch_procedures(PHY_VARS_NR_UE *ue,
                       ue->frame_parms.nb_antennas_rx,
                       rxdataF_comp,
                       ptrs_phase_per_slot,
-                      ptrs_re_per_slot)
+                      ptrs_re_per_slot,
+                      nvar)
           < 0)
         return -1;