Common Phase error estimation and compensation in uplink.

FEATURE STATEMENT:
- Introduce linear phase error noise model in Uplink at UE
- Perform common phase error (CPE) estimation and compensation at gNB

SOLUTION:
- A linear phase shift model is introduced in simulation.
- PTRS symbols are used to perform estimation of CPE from DMRS compensated signal
- The estimated values are interpolated in time domain and signal is compensated for the CPE.
- PTRS processing is done in Frequency Domain for each symbol in a slot and
  LLR's are calculated for each symbol accordingly.

IMPLEMENTATION:
* sim.h/channle_sim.c
- Linear Phase Noise Generation model definition.
* nr_ul_channel_estimation.c/nr_ul_estimation.h
- CPE estimation from PTRS and DMRS compensated signal.
* nr_dmrs_rx.c/nr_refsig.h
- Regenerate PTRS symbols at gNB.
* nr_ulsch_demodulation.c
- Removed old PTRS processing code and move to a common PTRS processing function
* defs_gNB.h/init.c
- New PTRS variables definition and initialization
* nr_ulsch_ue.c
- Corrected PTRS parameter to get new PTRS symbols for each OFDM symbol

TESTING
* ulsim.c
- Added Phase noise, Enable PTRS signal and verified the output.

VERIFICATION
- The LLR are rotated back with estimated CPE and no error is observed in scrambling/decoding

cmake_targets/CMakeLists.txt

@@ -2303,6 +2303,7 @@ add_library(LFDS7
 add_library(SIMU_COMMON
   ${OPENAIR1_DIR}/SIMULATION/TOOLS/random_channel.c
   ${OPENAIR1_DIR}/SIMULATION/TOOLS/rangen_double.c
+  ${OPENAIR1_DIR}/SIMULATION/TOOLS/phase_noise.c
   )
 
 # Simulation library

cmake_targets/build_oai

@@ -209,7 +209,7 @@ function main() {
                     GDB=0
                     CMAKE_BUILD_TYPE="RelWithDebInfo"
                     echo_info "Will Compile with gdb symbols"
-                    CMAKE_CMD="$CMAKE_CMD -DCMAKE_BUILD_TYPE=RelWithDebInfo"
+                    CMAKE_CMD="$CMAKE_CMD -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_EXPORT_COMPILE_COMMANDS=1"
                     shift
                     ;;
                 "MinSizeRel")

openair1/PHY/INIT/nr_init.c

@@ -247,10 +247,13 @@ int phy_init_nr_gNB(PHY_VARS_gNB *gNB,
     pusch_vars[ULSCH_id]->ul_ch_estimates_ext   = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
     pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates     = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
     pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates_ext = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
+    pusch_vars[ULSCH_id]->ptrs_phase_per_slot   = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
+    pusch_vars[ULSCH_id]->ptrs_valid_re_per_slot= (int16_t **)malloc16(Prx*sizeof(int16_t *) );
     pusch_vars[ULSCH_id]->ul_ch_estimates_time  = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
     pusch_vars[ULSCH_id]->rxdataF_comp          = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
-    pusch_vars[ULSCH_id]->ul_ch_mag0             = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
-    pusch_vars[ULSCH_id]->ul_ch_magb0            = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
+    pusch_vars[ULSCH_id]->rxdataF_ptrs_comp     = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
+    pusch_vars[ULSCH_id]->ul_ch_mag0            = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
+    pusch_vars[ULSCH_id]->ul_ch_magb0           = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
     pusch_vars[ULSCH_id]->ul_ch_mag             = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
     pusch_vars[ULSCH_id]->ul_ch_magb            = (int32_t **)malloc16(Prx*sizeof(int32_t *) );
     pusch_vars[ULSCH_id]->rho                   = (int32_t **)malloc16_clear(Prx*sizeof(int32_t*) );
@@ -264,9 +267,12 @@ int phy_init_nr_gNB(PHY_VARS_gNB *gNB,
       pusch_vars[ULSCH_id]->ul_ch_estimates_time[i]  = (int32_t *)malloc16_clear( 2*sizeof(int32_t)*fp->ofdm_symbol_size );
       pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates[i]       = (int32_t *)malloc16_clear( sizeof(int32_t)*fp->ofdm_symbol_size*2*fp->symbols_per_slot ); // max intensity in freq is 1 sc every 2 RBs
       pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates_ext[i]   = (int32_t *)malloc16_clear( sizeof(int32_t)*fp->ofdm_symbol_size*2*fp->symbols_per_slot );
+      pusch_vars[ULSCH_id]->ptrs_phase_per_slot[i]   = (int32_t *)malloc16_clear( sizeof(int32_t)*fp->symbols_per_slot); // symbols per slot
+      pusch_vars[ULSCH_id]->ptrs_valid_re_per_slot[i]= (int16_t *)malloc16_clear( sizeof(int16_t)*fp->symbols_per_slot);
       pusch_vars[ULSCH_id]->rxdataF_comp[i]          = (int32_t *)malloc16_clear( sizeof(int32_t)*N_RB_UL*12*fp->symbols_per_slot );
-      pusch_vars[ULSCH_id]->ul_ch_mag0[i]             = (int32_t *)malloc16_clear( fp->symbols_per_slot*sizeof(int32_t)*N_RB_UL*12 );
-      pusch_vars[ULSCH_id]->ul_ch_magb0[i]            = (int32_t *)malloc16_clear( fp->symbols_per_slot*sizeof(int32_t)*N_RB_UL*12 );
+      pusch_vars[ULSCH_id]->rxdataF_ptrs_comp[i]     = (int32_t *)malloc16_clear( sizeof(int32_t)*N_RB_UL*12*fp->symbols_per_slot );
+      pusch_vars[ULSCH_id]->ul_ch_mag0[i]            = (int32_t *)malloc16_clear( fp->symbols_per_slot*sizeof(int32_t)*N_RB_UL*12 );
+      pusch_vars[ULSCH_id]->ul_ch_magb0[i]           = (int32_t *)malloc16_clear( fp->symbols_per_slot*sizeof(int32_t)*N_RB_UL*12 );
       pusch_vars[ULSCH_id]->ul_ch_mag[i]             = (int32_t *)malloc16_clear( fp->symbols_per_slot*sizeof(int32_t)*N_RB_UL*12 );
       pusch_vars[ULSCH_id]->ul_ch_magb[i]            = (int32_t *)malloc16_clear( fp->symbols_per_slot*sizeof(int32_t)*N_RB_UL*12 );
       pusch_vars[ULSCH_id]->rho[i]                   = (int32_t *)malloc16_clear( sizeof(int32_t)*(fp->N_RB_UL*12*7*2) );
@@ -334,7 +340,10 @@ void phy_free_nr_gNB(PHY_VARS_gNB *gNB)
       free_and_zero(pusch_vars[ULSCH_id]->ul_ch_estimates_time[i]);
       free_and_zero(pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates[i]);
       free_and_zero(pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates_ext[i]);
+      free_and_zero(pusch_vars[ULSCH_id]->ptrs_phase_per_slot[i]);
+      free_and_zero(pusch_vars[ULSCH_id]->ptrs_valid_re_per_slot[i]);
       free_and_zero(pusch_vars[ULSCH_id]->rxdataF_comp[i]);
+      free_and_zero(pusch_vars[ULSCH_id]->rxdataF_ptrs_comp[i]);
       free_and_zero(pusch_vars[ULSCH_id]->ul_ch_mag0[i]);
       free_and_zero(pusch_vars[ULSCH_id]->ul_ch_magb0[i]);
       free_and_zero(pusch_vars[ULSCH_id]->ul_ch_mag[i]);
@@ -349,7 +358,10 @@ void phy_free_nr_gNB(PHY_VARS_gNB *gNB)
     free_and_zero(pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates);
     free_and_zero(pusch_vars[ULSCH_id]->ul_ch_ptrs_estimates_ext);
     free_and_zero(pusch_vars[ULSCH_id]->ul_ch_estimates_time);
+    free_and_zero(pusch_vars[ULSCH_id]->ptrs_phase_per_slot);
+    free_and_zero(pusch_vars[ULSCH_id]->ptrs_valid_re_per_slot);
     free_and_zero(pusch_vars[ULSCH_id]->rxdataF_comp);
+    free_and_zero(pusch_vars[ULSCH_id]->rxdataF_ptrs_comp);
     free_and_zero(pusch_vars[ULSCH_id]->ul_ch_mag0);
     free_and_zero(pusch_vars[ULSCH_id]->ul_ch_magb0);
     free_and_zero(pusch_vars[ULSCH_id]->ul_ch_mag);

openair1/PHY/NR_ESTIMATION/nr_ul_channel_estimation.c

@@ -25,6 +25,8 @@
 #include "nr_ul_estimation.h"
 #include "PHY/sse_intrin.h"
 #include "PHY/NR_REFSIG/nr_refsig.h"
+#include "PHY/NR_REFSIG/ptrs_nr.h"
+#include "PHY/NR_TRANSPORT/nr_transport_proto.h"
 #include "PHY/NR_UE_ESTIMATION/filt16a_32.h"
 
 //#define DEBUG_CH
@@ -501,3 +503,406 @@ int nr_pusch_channel_estimation(PHY_VARS_gNB *gNB,
 
   return(0);
 }
+
+
+/*******************************************************************
+ *
+ * NAME :         nr_pusch_ptrs_processing
+ *
+ * PARAMETERS :   gNB         : gNB data structure
+ *                rel15_ul    : UL parameters
+ *                UE_id       : UE ID
+ *                nr_tti_rx   : slot rx TTI
+ *            dmrs_symbol_flag: DMRS Symbol Flag
+ *                symbol      : OFDM Symbol
+ *                nb_re_pusch : PUSCH RE's
+ *                nb_re_pusch : PUSCH RE's
+ *
+ * RETURN :       nothing
+ *
+ * DESCRIPTION :
+ *  If ptrs is enabled process the symbol accordingly
+ *  1) Estimate phase noise per PTRS symbol
+ *  2) Interpolate PTRS estimated value in TD after all PTRS symbols
+ *  3) Compensated DMRS based estimated signal with PTRS estimation for slot
+ *********************************************************************/
+void nr_pusch_ptrs_processing(PHY_VARS_gNB *gNB,
+                              nfapi_nr_pusch_pdu_t *rel15_ul,
+                              uint8_t ulsch_id,
+                              uint8_t nr_tti_rx,
+                              uint8_t dmrs_symbol_flag,
+                              unsigned char symbol,
+                              uint32_t nb_re_pusch)
+{
+  NR_DL_FRAME_PARMS *frame_parms = &gNB->frame_parms;
+  int16_t *phase_per_symbol;
+  int16_t *no_re_per_symbol;
+
+  uint8_t         L_ptrs          = 0;
+  uint8_t         right_side_ref  = 0;
+  uint8_t         left_side_ref   = 0;
+  uint8_t         nb_dmrs_in_slot = 0;
+
+  //#define DEBUG_UL_PTRS 1
+  /* First symbol calculate PTRS symbol index for slot & set the variables */
+  if(symbol == rel15_ul->start_symbol_index)
+  {
+    gNB->pusch_vars[ulsch_id]->ptrs_symbols = 0;
+    L_ptrs = 1<<(rel15_ul->pusch_ptrs.ptrs_time_density);
+    set_ptrs_symb_idx(&gNB->pusch_vars[ulsch_id]->ptrs_symbols,
+                      rel15_ul->nr_of_symbols,
+                      rel15_ul->start_symbol_index,
+                      L_ptrs,
+                      rel15_ul->ul_dmrs_symb_pos);
+  }/* First symbol check */
+
+  /* loop over antennas */
+  for (int aarx=0; aarx< frame_parms->nb_antennas_rx; aarx++)
+  {
+    phase_per_symbol = (int16_t*)gNB->pusch_vars[ulsch_id]->ptrs_phase_per_slot[aarx];
+    no_re_per_symbol = gNB->pusch_vars[ulsch_id]->ptrs_valid_re_per_slot[aarx];
+    no_re_per_symbol[symbol] = nb_re_pusch;
+    /* if not PTRS symbol set current ptrs symbol index to zero*/
+    gNB->pusch_vars[ulsch_id]->ptrs_symbol_index = 0;
+    gNB->pusch_vars[ulsch_id]->ptrs_sc_per_ofdm_symbol = 0;
+    /* Check if current symbol contains PTRS */
+    if(is_ptrs_symbol(symbol, gNB->pusch_vars[ulsch_id]->ptrs_symbols))
+    {
+      gNB->pusch_vars[ulsch_id]->ptrs_symbol_index = symbol;
+      /*------------------------------------------------------------------------------------------------------- */
+      /* 1) Estimate phase noise per PTRS symbol                                                                */
+      /*------------------------------------------------------------------------------------------------------- */
+      nr_pusch_phase_estimation(frame_parms,
+                                rel15_ul,
+                                (int16_t *)&gNB->pusch_vars[ulsch_id]->ul_ch_ptrs_estimates_ext[aarx][symbol*nb_re_pusch],
+                                nr_tti_rx,
+                                symbol,
+                                (int16_t*)&gNB->pusch_vars[ulsch_id]->rxdataF_comp[aarx][(symbol * nb_re_pusch)],
+                                gNB->nr_gold_pusch_dmrs[rel15_ul->scid],
+                                &phase_per_symbol[2* symbol],
+                                &gNB->pusch_vars[ulsch_id]->ptrs_sc_per_ofdm_symbol);
+      /*  Subtract total PTRS RE's in the symbol from PUSCH RE's */
+      no_re_per_symbol[symbol] = nb_re_pusch - gNB->pusch_vars[ulsch_id]->ptrs_sc_per_ofdm_symbol;
+    }
+    /* DMRS Symbol channel estimates extraction */
+    else if(dmrs_symbol_flag)
+    {
+      phase_per_symbol[2* symbol]= (int16_t)((1<<15)-1); // 32767
+      phase_per_symbol[2* symbol +1]= 0;// no angle
+    }
+    /* For last OFDM symbol at each antenna perform interpolation and compensation for the slot*/
+    if(symbol == (rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols -1))
+    {
+      nb_dmrs_in_slot = get_dmrs_symbols_in_slot(rel15_ul->ul_dmrs_symb_pos,(rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols));
+      for(uint8_t dmrs_sym = 0; dmrs_sym < nb_dmrs_in_slot;  dmrs_sym ++)
+      {
+        if(dmrs_sym == 0)
+        {
+          /* get first DMRS position */
+          left_side_ref = get_next_dmrs_symbol_in_slot(rel15_ul->ul_dmrs_symb_pos, rel15_ul->start_symbol_index, (rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols));
+          /* get first DMRS position is not at start symbol position then we need to extrapolate left side  */
+          if(left_side_ref > rel15_ul->start_symbol_index)
+          {
+            left_side_ref = rel15_ul->start_symbol_index;
+          }
+        }
+        /* get the next symbol from left_side_ref value */
+        right_side_ref = get_next_dmrs_symbol_in_slot(rel15_ul->ul_dmrs_symb_pos, left_side_ref+1, (rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols));
+        /* if no symbol found then interpolate till end of slot*/
+        if(right_side_ref == 0)
+        {
+          right_side_ref = (rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols);
+        }
+        /*------------------------------------------------------------------------------------------------------- */
+        /* 2) Interpolate PTRS estimated value in TD */
+        /*------------------------------------------------------------------------------------------------------- */
+        nr_pusch_phase_interpolation(phase_per_symbol,left_side_ref,right_side_ref);
+        /* set left to last dmrs */
+        left_side_ref = right_side_ref;
+      } /*loop over dmrs positions */
+
+      //nr_pusch_phase_interpolation(phase_per_symbol,rel15_ul->start_symbol_index,frame_parms->symbols_per_slot );
+#ifdef DEBUG_UL_PTRS
+      LOG_M("ptrsEst.m","est",gNB->pusch_vars[ulsch_id]->ptrs_phase_per_slot[aarx],frame_parms->symbols_per_slot,1,1 );
+#endif
+
+      /*------------------------------------------------------------------------------------------------------- */
+      /* 3) Compensated DMRS based estimated signal with PTRS estimation                                        */
+      /*--------------------------------------------------------------------------------------------------------*/
+      for(uint8_t i =rel15_ul->start_symbol_index; i< (rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols);i++)
+      {
+#ifdef DEBUG_UL_PTRS
+        printf("PTRS: Rotate Symbol %2d with  %d + j* %d\n", i, phase_per_symbol[2* i],phase_per_symbol[(2* i) +1]);
+#endif
+        rotate_cpx_vector((int16_t*)&gNB->pusch_vars[ulsch_id]->rxdataF_comp[aarx][(i * rel15_ul->rb_size * NR_NB_SC_PER_RB)],
+                          &phase_per_symbol[2* i],
+                          (int16_t*)&gNB->pusch_vars[ulsch_id]->rxdataF_ptrs_comp[aarx][(i * rel15_ul->rb_size * NR_NB_SC_PER_RB)],
+                          (rel15_ul->rb_size * NR_NB_SC_PER_RB),
+                          15);
+      }// symbol loop
+    }//interpolation and compensation
+  }// Antenna loop
+}
+
+/*******************************************************************
+ *
+ * NAME :         nr_pusch_phase_estimation
+ *
+ * PARAMETERS :   frame_parms  : UL frame parameters
+ *                rel15_ul     : UL PDU Structure
+ *                Ns           :
+ *                Symbol       : OFDM symbol index
+ *                rxF          : Channel compensated signal
+ *                ptrs_gold_seq: Gold sequence for PTRS regeneration
+ *                error_est    : Estimated error output vector [Re Im]
+ * RETURN :       nothing
+ *
+ * DESCRIPTION :
+ *  perform phase estimation from regenerated PTRS SC and channel compensated
+ *  signal
+ *********************************************************************/
+void nr_pusch_phase_estimation(NR_DL_FRAME_PARMS *frame_parms,
+                               nfapi_nr_pusch_pdu_t *rel15_ul,
+                               int16_t *ptrs_ch_p,
+                               unsigned char Ns,
+                               unsigned char symbol,
+                               int16_t *rxF_comp,
+                               uint32_t ***ptrs_gold_seq,
+                               int16_t *error_est,
+                               uint16_t *ptrs_sc)
+{
+  uint8_t               is_ptrs_re       = 0;
+  uint16_t              re_cnt           = 0;
+  uint16_t              cnt              = 0;
+  unsigned short        nb_re_pusch      = NR_NB_SC_PER_RB * rel15_ul->rb_size;
+  uint8_t               K_ptrs           = (rel15_ul->pusch_ptrs.ptrs_freq_density)?4:2;
+  uint16_t              sc_per_symbol    = rel15_ul->rb_size/K_ptrs;
+  int16_t              *ptrs_p           = (int16_t *)malloc(sizeof(int32_t)*(sc_per_symbol));
+  int16_t              *dmrs_comp_p      = (int16_t *)malloc(sizeof(int32_t)*(sc_per_symbol));
+  double                abs              = 0.0;
+  double                real             = 0.0;
+  double                imag             = 0.0;
+#ifdef DEBUG_UL_PTRS
+  double                alpha            = 0;
+#endif
+  /* generate PTRS RE for the symbol */
+  nr_ptrs_rx_gen(ptrs_gold_seq[Ns][symbol],rel15_ul->rb_size,ptrs_p);
+
+  /* loop over all sub carriers to get compensated RE on ptrs symbols*/
+  for (int re = 0; re < nb_re_pusch; re++)
+  {
+    is_ptrs_re = is_ptrs_subcarrier(re,
+                                    rel15_ul->rnti,
+                                    0,
+                                    rel15_ul->dmrs_config_type,
+                                    K_ptrs,
+                                    rel15_ul->rb_size,
+                                    rel15_ul->pusch_ptrs.ptrs_ports_list[0].ptrs_re_offset,
+                                    0,// start_re is 0 here
+                                    frame_parms->ofdm_symbol_size);
+    if(is_ptrs_re)
+    {
+      dmrs_comp_p[re_cnt*2]     = rxF_comp[re *2];
+      dmrs_comp_p[(re_cnt*2)+1] = rxF_comp[(re *2)+1];
+      re_cnt++;
+    }
+    else
+    {
+      /* Skip PTRS symbols and keep data in a continuous vector */
+      rxF_comp[cnt *2]= rxF_comp[re *2];
+      rxF_comp[(cnt *2)+1]= rxF_comp[(re *2)+1];
+      cnt++;
+    }
+  }/* RE loop */
+  /* update the total ptrs RE in a symbol */
+  *ptrs_sc = re_cnt;
+
+  /*Multiple compensated data with conj of PTRS */
+  mult_cpx_vector(dmrs_comp_p, ptrs_p, ptrs_ch_p,(1 + sc_per_symbol/4)*4,15); // 2^15 shifted
+
+  /* loop over all ptrs sub carriers in a symbol */
+  /* sum the error vector */
+  for(int i = 0;i < sc_per_symbol; i++)
+  {
+    real+= ptrs_ch_p[(2*i)];
+    imag+= ptrs_ch_p[(2*i)+1];
+  }
+#ifdef DEBUG_UL_PTRS
+    alpha = atan(imag/real);
+    printf("PTRS: Symbol  %d atan(Im,real):= %f \n",symbol, alpha );
+#endif
+  /* mean */
+  real /= sc_per_symbol;
+  imag /= sc_per_symbol;
+  /* absolute calculation */
+  abs = sqrt(((real * real) + (imag *  imag)));
+  /* normalized error estimation */
+  error_est[0]= (real / abs)*(1<<15);
+  /* compensation in given by conjugate of estimated phase (e^-j*2*pi*fd*t)*/
+  error_est[1]= (-1)*(imag / abs)*(1<<15);
+#ifdef DEBUG_UL_PTRS
+    printf("PTRS: Estimated Symbol  %d -> %d + j* %d \n",symbol, error_est[0], error_est[1] );
+#endif
+  /* free vectors */
+  free(ptrs_p);
+  free(dmrs_comp_p);
+}
+
+
+/*******************************************************************
+ *
+ * NAME :         nr_pusch_phase_interpolation
+ *
+ * PARAMETERS :   *error_est    : Data Pointer [Re Im Re Im ...]
+ *                 start_symbol : Start Symbol
+ *                 end_symbol   : End Symbol
+ * RETURN :       nothing
+ *
+ * DESCRIPTION :
+ * Perform Interpolation, extrapolation based upon the estimation
+ * location between the data Pointer Array.
+ *
+ *********************************************************************/
+void nr_pusch_phase_interpolation(int16_t *error_est,
+                                  uint8_t start_symbol,
+                                  uint8_t end_symbol
+                                  )
+{
+
+  int next = 0, prev = 0, candidates= 0, distance=0, leftEdge= 0, rightEdge = 0, getDiff =0 ;
+  double weight = 0.0;
+  double scale  = 0.125 ; // to avoid saturation due to fixed point multiplication
+#ifdef DEBUG_UL_PTRS
+  printf("PTRS: INT: Left limit %d, Right limit %d, Loop over %d Symbols \n",
+         start_symbol,end_symbol-1, (end_symbol -start_symbol)-1);
+#endif
+  for(int i =start_symbol; i< end_symbol;i++)
+  {
+    /* Only update when an estimation is found */
+    if( error_est[i*2] != 0 )
+    {
+      /* if found a symbol then set next symbol also */
+      next = nr_ptrs_find_next_estimate(error_est, i, end_symbol);
+      /* left extrapolation, if first estimate value is zero */
+      if( error_est[i*2] == 0 )
+      {
+        leftEdge = 1;
+      }
+      /* right extrapolation, if next is 0 before end symbol */
+      if((next == 0) && (end_symbol > i))
+      {
+        rightEdge = 1;
+        /* special case as no right extrapolation possible with DMRS on left */
+        /* In this case take mean of most recent 2 estimated points */
+        if(prev ==0)
+        {
+          prev = start_symbol -1;
+          next = start_symbol -2;
+          getDiff =1;
+        }else
+        {
+          /* for right edge  previous is second last from right side */
+          next = prev;
+          /* Set the current as recent estimation reference */
+          prev = i;
+        }
+      }
+      /* update  current symbol as prev  for next symbol */
+      if (rightEdge==0)
+        /* Set the current as recent estimation reference */
+        prev = i;
+    }
+    /*extrapolation left side*/
+    if(leftEdge)
+    {
+      distance = next - prev;
+      weight = 1.0/distance;
+      candidates = i;
+      for(int j = 1; j <= candidates; j++)
+      {
+        error_est[(i-j)*2]    = 8 *(((double)(error_est[prev*2]) * scale * (distance + j) * weight) -
+                                    ((double)(error_est[next*2]) * scale * j * weight));
+        error_est[((i-j)*2)+1]= 8 *(((double)(error_est[(prev*2)+1]) * scale* (distance + j) * weight) -
+                                    ((double)(error_est[((next*2)+1)]) * scale * j * weight));
+#ifdef DEBUG_UL_PTRS
+        printf("PTRS: INT: Left Edge i= %d weight= %f %d + j*%d, Prev %d Next %d \n",
+               (i-j),weight, error_est[(i-j)*2],error_est[((i-j)*2)+1], prev,next);
+#endif
+      }
+      leftEdge = 0;
+    }
+    /* extrapolation at right side */
+    else if (rightEdge )
+    {
+      if(getDiff)
+      {
+        error_est[(i+1)*2]    = ((1<<15) +(error_est[prev*2]) - error_est[next*2]);
+        error_est[((i+1)*2)+1]= error_est[(prev*2)+1] - error_est[(next*2)+1];
+#ifdef DEBUG_UL_PTRS
+        printf("PTRS: INT: Right Edge Special Case i= %d weight= %f %d + j*%d, Prev %d Next %d \n",
+               (i+1),weight, error_est[(i+1)*2],error_est[((i+1)*2)+1], prev,next);
+#endif
+        i++;
+      }
+      else
+      {
+        distance = prev - next;
+        candidates = (end_symbol -1) - i;
+        weight = 1.0/distance;
+        for(int j = 1; j <= candidates; j++)
+        {
+          error_est[(i+j)*2]    =  8 *(((double)(error_est[prev*2]) * scale * (distance + j) * weight) -
+                                       ((double)(error_est[next*2]) * scale * j * weight));
+          error_est[((i+j)*2)+1]=  8 *(((double)(error_est[(prev*2)+1]) * scale * (distance + j) * weight) -
+                                       ((double)(error_est[((next*2)+1)]) * scale *j * weight));
+#ifdef DEBUG_UL_PTRS
+          printf("PTRS: INT: Right Edge i= %d weight= %f %d + j*%d, Prev %d Next %d \n",
+                 (i+j),weight, error_est[(i+j)*2],error_est[((i+j)*2)+1], prev,next);
+#endif
+        }
+        if(candidates > 1)
+        {
+          i+=candidates;
+        }
+      }
+    }
+    /* Interpolation between 2 estimated points */
+    else if(next != 0 && ( error_est[2*i] == 0 ))
+    {
+      distance = next - prev;
+      weight = 1.0/distance;
+      candidates = next - i ;
+      for(int j = 0; j < candidates; j++)
+      {
+
+        error_est[(i+j)*2]    = 8 *(((double)(error_est[prev*2]) * scale * (distance - (j+1)) * weight) +
+                                    ((double)(error_est[next*2]) * scale * (j+1) * weight));
+        error_est[((i+j)*2)+1]= 8 *(((double)(error_est[(prev*2)+1]) * scale *(distance - (j+1)) * weight) +
+                                    ((double)(error_est[((next*2)+1)]) * scale *(j+1) * weight));
+#ifdef DEBUG_UL_PTRS
+        printf("PTRS: INT: Interpolation i= %d weight= %f %d + j*%d, Prev %d Next %d\n",
+               (i+j),weight, error_est[(i+j)*2],error_est[((i+j)*2)+1],prev,next);
+#endif
+      }
+      if(candidates > 1)
+      {
+        i+=candidates-1;
+      }
+    }// interpolation
+  }// symbol loop
+}
+
+/* Find the next non zero Real value in a complex vector */
+int nr_ptrs_find_next_estimate(int16_t *error_est,
+                               uint8_t counter,
+                               uint8_t end_symbol)
+{
+  for (int i = counter +1 ; i< end_symbol; i++)
+  {
+    if( error_est[2*i] != 0)
+    {
+      return i;
+    }
+  }
+  return 0;
+}

openair1/PHY/NR_ESTIMATION/nr_ul_estimation.h

@@ -48,4 +48,30 @@ int nr_pusch_channel_estimation(PHY_VARS_gNB *gNB,
 
 int nr_est_timing_advance_pusch(PHY_VARS_gNB* phy_vars_gNB, int UE_id);
 
+
+void nr_pusch_ptrs_processing(PHY_VARS_gNB *gNB,
+                              nfapi_nr_pusch_pdu_t *rel15_ul,
+                              uint8_t ulsch_id,
+                              uint8_t nr_tti_rx,
+                              uint8_t dmrs_symbol_flag,
+                              unsigned char symbol,
+                              uint32_t nb_re_pusch);
+
+void nr_pusch_phase_estimation(NR_DL_FRAME_PARMS *frame_parms,
+                               nfapi_nr_pusch_pdu_t *rel15_ul,
+                               int16_t *ptrs_ch_p,
+                               unsigned char Ns,
+                               unsigned char symbol,
+                               int16_t *rxF_comp,
+                               uint32_t ***ptrs_gold_seq,
+                               int16_t *error_est,
+                               uint16_t  *ptrs_sc);
+
+void nr_pusch_phase_interpolation(int16_t *error_est,
+                                  uint8_t start_symbol,
+                                  uint8_t end_symbol);
+
+int nr_ptrs_find_next_estimate(int16_t *error_est,
+                               uint8_t counter,
+                               uint8_t end_symbol);
 #endif

openair1/PHY/NR_REFSIG/dmrs_nr.c

@@ -362,4 +362,25 @@ void generate_dmrs_pbch(uint32_t dmrs_pbch_bitmap[DMRS_PBCH_I_SSB][DMRS_PBCH_N_H
   free(dmrs_sequence);
 #endif
 }
+/* return the position of next dmrs symbol in a slot */
+int get_next_dmrs_symbol_in_slot(uint16_t  ul_dmrs_symb_pos, uint8_t counter, uint8_t end_symbol)
+{
+  for(uint8_t symbol = counter; symbol < end_symbol; symbol++)
+    if((ul_dmrs_symb_pos >>symbol)&0x01 )
+    {
+      return symbol;
+    }
+  return 0;
+}
+
 
+/* return the total number of dmrs symbol in a slot */
+uint8_t get_dmrs_symbols_in_slot(uint16_t l_prime_mask,  uint16_t nb_symb)
+{
+  uint8_t tmp = 0;
+  for (int i = 0; i < nb_symb; i++)
+  {
+    tmp += (l_prime_mask >> i) & 0x01;
+  }
+  return tmp;
+}

openair1/PHY/NR_REFSIG/nr_dmrs_rx.c

@@ -223,3 +223,30 @@ int nr_pbch_dmrs_rx(int symbol,
   
   return(0);
 }
+
+/*!
+  \brief This function generate gold ptrs sequence for each OFDM symbol
+  \param *in gold sequence for ptrs per OFDM symbol
+  \param length is number of RE in a OFDM symbol
+  \param *output pointer to all ptrs RE in a OFDM symbol
+*/
+void nr_ptrs_rx_gen(uint32_t *in, uint32_t length, int16_t *output)
+{
+  uint16_t offset;
+  uint8_t idx, b_idx;
+  int mod_order = 2;
+  offset = NR_MOD_TABLE_QPSK_OFFSET;
+  for (int i=0; i<length/mod_order; i++)
+    {
+      idx = 0;
+      for (int j=0; j<mod_order; j++)
+        {
+          b_idx = (i*mod_order+j)&0x1f;
+          if (i && (!b_idx))
+            in++;
+          idx ^= (((*in)>>b_idx)&1)<<(mod_order-j-1);
+        }
+      output[i<<1] = nr_rx_mod_table[(offset+idx)<<1];
+      output[(i<<1)+1] =  nr_rx_mod_table[((offset+idx)<<1)+1];
+    }
+}

openair1/PHY/NR_REFSIG/nr_refsig.h

@@ -52,6 +52,9 @@ int nr_pusch_dmrs_rx(PHY_VARS_gNB *gNB,
                      uint8_t dmrs_type);
 
 void init_scrambling_luts(void);
+void nr_ptrs_rx_gen(uint32_t *in, uint32_t length, int16_t *output);
+uint8_t get_next_dmrs_symbol_in_slot(uint16_t  ul_dmrs_symb_pos, uint8_t counter, uint8_t end_symbol);
+uint8_t get_dmrs_symbols_in_slot(uint16_t l_prime_mask,  uint16_t nb_symb);
 
 extern __m64 byte2m64_re[256];
 extern __m64 byte2m64_im[256];

openair1/PHY/NR_TRANSPORT/nr_ulsch_demodulation.c

@@ -232,18 +232,13 @@ void nr_ulsch_extract_rbs_single(int32_t **rxdataF,
 
   unsigned short start_re, re, nb_re_pusch;
   unsigned char aarx;
-  uint8_t K_ptrs = 0;
   uint32_t rxF_ext_index = 0;
   uint32_t ul_ch0_ext_index = 0;
   uint32_t ul_ch0_index = 0;
-  uint32_t ul_ch0_ptrs_ext_index = 0;
-  uint32_t ul_ch0_ptrs_index = 0;
-  uint8_t is_ptrs_symbol_flag,k_prime;
-  uint16_t n=0, num_ptrs_symbols;
+  uint8_t k_prime;
+  uint16_t n=0;
   int16_t *rxF,*rxF_ext;
   int *ul_ch0,*ul_ch0_ext;
-  int *ul_ch0_ptrs,*ul_ch0_ptrs_ext;
-  uint16_t n_rnti = pusch_pdu->rnti;
   uint8_t delta = 0;
 
 #ifdef DEBUG_RB_EXT
@@ -252,15 +247,13 @@ void nr_ulsch_extract_rbs_single(int32_t **rxdataF,
   printf("--------------------ch_ext_index = %d-----------------------\n", symbol*NR_NB_SC_PER_RB * pusch_pdu->rb_size);
 
 #endif
-  
+
   uint8_t is_dmrs_re;
   start_re = (frame_parms->first_carrier_offset + (pusch_pdu->rb_start * NR_NB_SC_PER_RB))%frame_parms->ofdm_symbol_size;
   nb_re_pusch = NR_NB_SC_PER_RB * pusch_pdu->rb_size;
-  is_ptrs_symbol_flag = 0;
-  num_ptrs_symbols = 0;
 
   for (aarx = 0; aarx < frame_parms->nb_antennas_rx; aarx++) {
-    
+
     rxF       = (int16_t *)&rxdataF[aarx][symbol * frame_parms->ofdm_symbol_size];
     rxF_ext   = (int16_t *)&pusch_vars->rxdataF_ext[aarx][symbol * nb_re_pusch]; // [hna] rxdataF_ext isn't contiguous in order to solve an alignment problem ib llr computation in case of mod_order = 4, 6
 
@@ -268,10 +261,6 @@ void nr_ulsch_extract_rbs_single(int32_t **rxdataF,
 
     ul_ch0_ext = &pusch_vars->ul_ch_estimates_ext[aarx][symbol*nb_re_pusch];
 
-    ul_ch0_ptrs     = &pusch_vars->ul_ch_ptrs_estimates[aarx][pusch_vars->ptrs_symbol_index*frame_parms->ofdm_symbol_size]; // update channel estimates if new dmrs symbol are available
-
-    ul_ch0_ptrs_ext = &pusch_vars->ul_ch_ptrs_estimates_ext[aarx][symbol*nb_re_pusch];
-
     n = 0;
     k_prime = 0;
 
@@ -282,54 +271,28 @@ void nr_ulsch_extract_rbs_single(int32_t **rxdataF,
       else
         is_dmrs_re = 0;
 
-      if (pusch_pdu->pdu_bit_map & PUSCH_PDU_BITMAP_PUSCH_PTRS) {
-        K_ptrs = (pusch_pdu->pusch_ptrs.ptrs_freq_density)?4:2;
-        if(is_ptrs_symbol(symbol, pusch_vars->ptrs_symbols))
-            is_ptrs_symbol_flag = is_ptrs_subcarrier((start_re + re) % frame_parms->ofdm_symbol_size,
-                                                     n_rnti,
-                                                     aarx,
-                                                     pusch_pdu->dmrs_config_type,
-                                                     K_ptrs,
-                                                     pusch_pdu->rb_size,
-                                                     pusch_pdu->pusch_ptrs.ptrs_ports_list[0].ptrs_re_offset,
-                                                     start_re,
-                                                     frame_parms->ofdm_symbol_size);
-
-        if (is_ptrs_symbol_flag == 1)
-          num_ptrs_symbols++;
-
-      }
-
-  #ifdef DEBUG_RB_EXT
-      printf("re = %d, is_ptrs_symbol_flag = %d, symbol = %d\n", re, is_ptrs_symbol_flag, symbol);
-  #endif
-
-      if ( is_dmrs_re == 0 && is_ptrs_symbol_flag == 0) {
+      /* save only data and respective channel estimates */
+      if ( is_dmrs_re == 0) {
 
-        rxF_ext[rxF_ext_index]     = (rxF[ ((start_re + re)*2)      % (frame_parms->ofdm_symbol_size*2)]);
-        rxF_ext[rxF_ext_index + 1] = (rxF[(((start_re + re)*2) + 1) % (frame_parms->ofdm_symbol_size*2)]);
+        rxF_ext[rxF_ext_index]       = (rxF[ ((start_re + re)*2)      % (frame_parms->ofdm_symbol_size*2)]);
+        rxF_ext[rxF_ext_index + 1]   = (rxF[(((start_re + re)*2) + 1) % (frame_parms->ofdm_symbol_size*2)]);
         ul_ch0_ext[ul_ch0_ext_index] = ul_ch0[ul_ch0_index];
-        ul_ch0_ptrs_ext[ul_ch0_ptrs_ext_index] = ul_ch0_ptrs[ul_ch0_ptrs_index];
 
-  #ifdef DEBUG_RB_EXT
+#ifdef DEBUG_RB_EXT
         printf("rxF_ext[%d] = (%d,%d)\n", rxF_ext_index>>1, rxF_ext[rxF_ext_index],rxF_ext[rxF_ext_index+1]);
-  #endif
+#endif
 
         ul_ch0_ext_index++;
-        ul_ch0_ptrs_ext_index++;
         rxF_ext_index +=2;
-      } else {
+      } else
+      {
         k_prime++;
         k_prime&=1;
         n+=(k_prime)?0:1;
       }
       ul_ch0_index++;
-      ul_ch0_ptrs_index++;
     }
   }
-
-  pusch_vars->ptrs_sc_per_ofdm_symbol = num_ptrs_symbols;
-
 }
 
 void nr_ulsch_scale_channel(int **ul_ch_estimates_ext,
@@ -1018,32 +981,19 @@ int nr_rx_pusch(PHY_VARS_gNB *gNB,
                 unsigned char harq_pid)
 {
 
-  uint8_t aarx, aatx, dmrs_symbol_flag, ptrs_symbol_flag; // dmrs_symbol_flag, a flag to indicate DMRS REs in current symbol
+  uint8_t aarx, aatx, dmrs_symbol_flag; // dmrs_symbol_flag, a flag to indicate DMRS REs in current symbol
   uint32_t nb_re_pusch, bwp_start_subcarrier;
-  uint8_t L_ptrs = 0; // PTRS parameter
   int avgs;
   int avg[4];
   NR_DL_FRAME_PARMS *frame_parms = &gNB->frame_parms;
   nfapi_nr_pusch_pdu_t *rel15_ul = &gNB->ulsch[ulsch_id][0]->harq_processes[harq_pid]->ulsch_pdu;
 
   dmrs_symbol_flag = 0;
-  ptrs_symbol_flag = 0;
 
   if(symbol == rel15_ul->start_symbol_index){
     gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset = 0;
     gNB->pusch_vars[ulsch_id]->dmrs_symbol = 0;
     gNB->pusch_vars[ulsch_id]->cl_done = 0;
-    gNB->pusch_vars[ulsch_id]->ptrs_symbols = 0;
-
-    if (rel15_ul->pdu_bit_map & PUSCH_PDU_BITMAP_PUSCH_PTRS) {  // if there is ptrs pdu
-      L_ptrs = 1<<(rel15_ul->pusch_ptrs.ptrs_time_density);
-
-      set_ptrs_symb_idx(&gNB->pusch_vars[ulsch_id]->ptrs_symbols,
-                        rel15_ul->nr_of_symbols,
-                        rel15_ul->start_symbol_index,
-                        L_ptrs,
-                        rel15_ul->ul_dmrs_symb_pos);
-    }
   }
 
   bwp_start_subcarrier = (rel15_ul->rb_start*NR_NB_SC_PER_RB + frame_parms->first_carrier_offset) % frame_parms->ofdm_symbol_size;
@@ -1067,15 +1017,6 @@ int nr_rx_pusch(PHY_VARS_gNB *gNB,
     nb_re_pusch = rel15_ul->rb_size * NR_NB_SC_PER_RB;
   }
 
-  if (rel15_ul->pdu_bit_map & PUSCH_PDU_BITMAP_PUSCH_PTRS) {  // if there is ptrs pdu
-    if(is_ptrs_symbol(symbol, gNB->pusch_vars[ulsch_id]->ptrs_symbols))
-    {
-      ptrs_symbol_flag = 1;
-      gNB->pusch_vars[ulsch_id]->ptrs_symbol_index = symbol;
-    }
-  }
-
-
   //----------------------------------------------------------
   //--------------------- Channel estimation ---------------------
   //----------------------------------------------------------
@@ -1103,11 +1044,9 @@ int nr_rx_pusch(PHY_VARS_gNB *gNB,
                                 frame_parms);
     stop_meas(&gNB->ulsch_rbs_extraction_stats);
 
-    if(ptrs_symbol_flag == 1)
-    {
-      /*  Subtract total PTRS RE's in the smybol from PUSCH RE's  */
-      nb_re_pusch -= gNB->pusch_vars[ulsch_id]->ptrs_sc_per_ofdm_symbol;
-    }
+    //----------------------------------------------------------
+    //--------------------- Channel Scaling --------------------
+    //----------------------------------------------------------
 
     nr_ulsch_scale_channel(gNB->pusch_vars[ulsch_id]->ul_ch_estimates_ext,
                            frame_parms,
@@ -1137,6 +1076,10 @@ int nr_rx_pusch(PHY_VARS_gNB *gNB,
       gNB->pusch_vars[ulsch_id]->cl_done = 1;
     }
 
+    //----------------------------------------------------------
+    //--------------------- Channel Compensation ---------------
+    //----------------------------------------------------------
+
     start_meas(&gNB->ulsch_channel_compensation_stats);
     nr_ulsch_channel_compensation(gNB->pusch_vars[ulsch_id]->rxdataF_ext,
                                   gNB->pusch_vars[ulsch_id]->ul_ch_estimates_ext,
@@ -1162,9 +1105,57 @@ int nr_rx_pusch(PHY_VARS_gNB *gNB,
     nr_idft(&((uint32_t*)gNB->pusch_vars[ulsch_id]->rxdataF_ext[0])[symbol * rel15_ul->rb_size * NR_NB_SC_PER_RB], nb_re_pusch);
 #endif
 
-  //----------------------------------------------------------
-  //-------------------- LLRs computation --------------------
-  //----------------------------------------------------------
+
+
+    //----------------------------------------------------------
+    //--------------------- PTRS Processing --------------------
+    //----------------------------------------------------------
+
+    /* In case PTRS is enabled then LLR will be calculated after PTRS symbols are processed *
+     * otherwise LLR are calculated for each symbol based upon DMRS channel estimates only. */
+    if (rel15_ul->pdu_bit_map & PUSCH_PDU_BITMAP_PUSCH_PTRS)
+    {
+      start_meas(&gNB->ulsch_ptrs_processing_stats);
+      nr_pusch_ptrs_processing(gNB,
+                               rel15_ul,
+                               ulsch_id,
+                               nr_tti_rx,
+                               dmrs_symbol_flag,
+                               symbol,
+                               nb_re_pusch);
+      stop_meas(&gNB->ulsch_ptrs_processing_stats);
+      /*---------------------------------------------------------------------------------------------------- */
+      /*  Calculate LLR based upon new compensated data after all PTRS symbols are processed                 */
+      /*-----------------------------------------------------------------------------------------------------*/
+      if(symbol == (rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols -1))
+      {
+        for (int aarx=0; aarx< frame_parms->nb_antennas_rx; aarx++)
+        {
+          for(uint8_t i =rel15_ul->start_symbol_index; i< (rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols);i++)
+          {
+            start_meas(&gNB->ulsch_llr_stats);
+            nr_ulsch_compute_llr(&gNB->pusch_vars[ulsch_id]->rxdataF_ptrs_comp[aarx][i * (rel15_ul->rb_size * NR_NB_SC_PER_RB)],
+                                 gNB->pusch_vars[ulsch_id]->ul_ch_mag0,
+                                 gNB->pusch_vars[ulsch_id]->ul_ch_magb0,
+                                 &gNB->pusch_vars[ulsch_id]->llr[gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset * rel15_ul->qam_mod_order],
+                                 rel15_ul->rb_size,
+                                 gNB->pusch_vars[ulsch_id]->ptrs_valid_re_per_slot[aarx][i],
+                                 i,
+                                 rel15_ul->qam_mod_order);
+            stop_meas(&gNB->ulsch_llr_stats);
+            gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset = gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset +  gNB->pusch_vars[ulsch_id]->ptrs_valid_re_per_slot[aarx][i];
+          }// symbol loop
+          /* For scope library rxdataF_comp is used so in PTRS case we point it to correct data pointer */
+          gNB->pusch_vars[ulsch_id]->rxdataF_comp[aarx] = gNB->pusch_vars[ulsch_id]->rxdataF_ptrs_comp[aarx];
+        }// antenna loop
+      }// last symbol check
+    }// PTRS processing
+    else // if no ptrs is enabled
+    {
+    //----------------------------------------------------------
+    //-------------------- LLRs computation --------------------
+    //----------------------------------------------------------
+
     start_meas(&gNB->ulsch_llr_stats);
     AssertFatal(gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset * rel15_ul->qam_mod_order+nb_re_pusch*rel15_ul->qam_mod_order < (8*((3*8*6144)+12)) , "Mysterious llr buffer size check");
     nr_ulsch_compute_llr(&gNB->pusch_vars[ulsch_id]->rxdataF_comp[0][symbol * rel15_ul->rb_size * NR_NB_SC_PER_RB],
@@ -1177,9 +1168,9 @@ int nr_rx_pusch(PHY_VARS_gNB *gNB,
                          rel15_ul->qam_mod_order);
     stop_meas(&gNB->ulsch_llr_stats);
 
+    gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset = gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset +  nb_re_pusch;
+    }/* End for PTRS Condition */
   }
-
-  gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset = gNB->pusch_vars[ulsch_id]->rxdataF_ext_offset +  nb_re_pusch;
   return (0);
 
 }

openair1/PHY/NR_UE_TRANSPORT/nr_ulsch_ue.c

@@ -372,7 +372,7 @@ void nr_ue_ulsch_procedures(PHY_VARS_NR_UE *UE,
     uint16_t m=0, n=0, dmrs_idx=0, ptrs_idx = 0;
 
     for (l=start_symbol; l<start_symbol+number_of_symbols; l++) {
-
+      ptrs_idx = 0; // every PTRS symbol has a new gold sequence
       k = start_sc;
       n = 0;
       dmrs_idx = 0;

openair1/PHY/defs_gNB.h

@@ -403,6 +403,10 @@ typedef struct {
   /// - first index: rx antenna id [0..nb_antennas_rx[
   /// - second index: ? [0..12*N_RB_UL*frame_parms->symbols_per_tti[
   int32_t **rxdataF_comp;
+  /// \brief Holds the PTRS compensated signal.
+  /// - first index: rx antenna id [0..nb_antennas_rx[
+  /// - second index: ? [0..12*N_RB_UL*frame_parms->symbols_per_tti[
+  int32_t **rxdataF_ptrs_comp;
   /// \brief Magnitude of the UL channel estimates. Used for 2nd-bit level thresholds in LLR computation
   /// - first index: rx antenna id [0..nb_antennas_rx[
   /// - second index: ? [0..12*N_RB_UL*frame_parms->symbols_per_tti[
@@ -446,6 +450,14 @@ typedef struct {
   uint16_t ptrs_symbols;
   // PTRS subcarriers per OFDM symbol
   uint16_t ptrs_sc_per_ofdm_symbol;
+  /// \brief Estimated phase error based upon PTRS on each symbol .
+  /// - first index: ? [0..7] Number of Antenna
+  /// - second index: ? [0...14] smybol per slot
+  int32_t **ptrs_phase_per_slot;
+  /// \brief Total RE count after PTRS RE's are extracted from respective symbol.
+  /// - first index: ? [0..7] Number of Antenna
+  /// - second index: ? [0...14] smybol per slot
+  int16_t **ptrs_valid_re_per_slot;
   /// flag to verify if channel level computation is done
   uint8_t cl_done;
 } NR_gNB_PUSCH;
@@ -775,9 +787,11 @@ typedef struct PHY_VARS_gNB_s {
   time_stats_t ulsch_deinterleaving_stats;
   time_stats_t ulsch_unscrambling_stats;
   time_stats_t ulsch_channel_estimation_stats;
+  time_stats_t ulsch_ptrs_processing_stats;
   time_stats_t ulsch_channel_compensation_stats;
   time_stats_t ulsch_rbs_extraction_stats;
   time_stats_t ulsch_llr_stats;
+
   /*
   time_stats_t rx_dft_stats;
   time_stats_t ulsch_freq_offset_estimation_stats;

openair1/SIMULATION/NR_PHY/ulsim.c

@@ -447,7 +447,7 @@ int main(int argc, char **argv)
 				  0);
   fix_scc(scc,ssb_bitmap);
 
-  xer_fprint(stdout, &asn_DEF_NR_CellGroupConfig, (const void*)secondaryCellGroup);
+  // xer_fprint(stdout, &asn_DEF_NR_CellGroupConfig, (const void*)secondaryCellGroup);
 
   AssertFatal((gNB->if_inst         = NR_IF_Module_init(0))!=NULL,"Cannot register interface");
 
@@ -562,14 +562,13 @@ int main(int argc, char **argv)
   uint8_t ptrs_time_density = get_L_ptrs(ptrs_mcs1, ptrs_mcs2, ptrs_mcs3, Imcs, mcs_table);
   uint8_t ptrs_freq_density = get_K_ptrs(n_rb0, n_rb1, nb_rb);
   int ptrs_symbols = 0; // to calculate total PTRS RE's in a slot
+  double ts = 1.0/(frame_parms->subcarrier_spacing * frame_parms->ofdm_symbol_size);
+
+  number_dmrs_symbols = get_dmrs_symbols_in_slot(l_prime_mask, nb_symb_sch);
 
   if(1<<ptrs_time_density >= nb_symb_sch)
     pdu_bit_map &= ~PUSCH_PDU_BITMAP_PUSCH_PTRS; // disable PUSCH PTRS
 
-  for (i = 0; i < nb_symb_sch; i++) {
-    number_dmrs_symbols += (l_prime_mask >> i) & 0x01;
-  }
-
   mod_order      = nr_get_Qm_ul(Imcs, 0);
   code_rate      = nr_get_code_rate_ul(Imcs, 0);
 
@@ -586,6 +585,7 @@ int main(int argc, char **argv)
       reset_meas(&gNB->ulsch_ldpc_decoding_stats);
       reset_meas(&gNB->ulsch_unscrambling_stats);
       reset_meas(&gNB->ulsch_channel_estimation_stats);
+      reset_meas(&gNB->ulsch_ptrs_processing_stats);
       reset_meas(&gNB->ulsch_llr_stats);
       reset_meas(&gNB->ulsch_channel_compensation_stats);
       reset_meas(&gNB->ulsch_rbs_extraction_stats);
@@ -765,19 +765,24 @@ int main(int argc, char **argv)
 
       int error_flag;
       for (trial = 0; trial < n_trials; trial++) {
-	
-	error_flag = 0;
-	  //----------------------------------------------------------
-	  //------------------------ add noise -----------------------
-	  //----------------------------------------------------------
-	if (input_fd == NULL ) {
-	  for (i=0; i<slot_length; i++) {
-	    for (ap=0; ap<frame_parms->nb_antennas_rx; ap++) {
-	      ((int16_t*) &gNB->common_vars.rxdata[ap][slot_offset])[(2*i) + (delay*2)]   = (int16_t)((double)(((int16_t *)&UE->common_vars.txdata[ap][slot_offset])[(i<<1)])/sqrt(scale)   + (sqrt(sigma/2)*gaussdouble(0.0,1.0))); // convert to fixed point
-	      ((int16_t*) &gNB->common_vars.rxdata[ap][slot_offset])[(2*i)+1 + (delay*2)]   = (int16_t)((double)(((int16_t *)&UE->common_vars.txdata[ap][slot_offset])[(i<<1)+1])/sqrt(scale) + (sqrt(sigma/2)*gaussdouble(0.0,1.0)));
-	    }
-	  }
-	}
+
+        error_flag = 0;
+        //----------------------------------------------------------
+        //------------------------ add noise -----------------------
+        //----------------------------------------------------------
+        if (input_fd == NULL ) {
+          for (i=0; i<slot_length; i++) {
+            for (ap=0; ap<frame_parms->nb_antennas_rx; ap++) {
+              /* Add phase noise if enabled */
+              if (pdu_bit_map & PUSCH_PDU_BITMAP_PUSCH_PTRS) {
+                phase_noise(ts, &((int16_t*)&UE->common_vars.txdata[ap][slot_offset])[(i<<1)],
+                            &((int16_t*)&UE->common_vars.txdata[ap][slot_offset])[(i<<1)+1]);
+              }
+              ((int16_t*) &gNB->common_vars.rxdata[ap][slot_offset])[(2*i) + (delay*2)]   = (int16_t)((double)(((int16_t *)&UE->common_vars.txdata[ap][slot_offset])[(i<<1)])/sqrt(scale)   + (sqrt(sigma/2)*gaussdouble(0.0,1.0))); // convert to fixed point
+              ((int16_t*) &gNB->common_vars.rxdata[ap][slot_offset])[(2*i)+1 + (delay*2)]   = (int16_t)((double)(((int16_t *)&UE->common_vars.txdata[ap][slot_offset])[(i<<1)+1])/sqrt(scale) + (sqrt(sigma/2)*gaussdouble(0.0,1.0)));
+            }
+          }
+        }
 	else {
 	  AssertFatal(frame_parms->nb_antennas_rx == 1, "nb_ant != 1\n");
 	  // 800 samples is N_TA_OFFSET for FR1 @ 30.72 Ms/s,
@@ -825,6 +830,10 @@ int main(int argc, char **argv)
 		&gNB->pusch_vars[0]->ul_ch_estimates_ext[0][(start_symbol+1)*NR_NB_SC_PER_RB * pusch_pdu->rb_size],(nb_symb_sch-1)*NR_NB_SC_PER_RB * pusch_pdu->rb_size,1,1);
 	  LOG_M("rxsigF0_comp.m","rxsF0_comp",
 		&gNB->pusch_vars[0]->rxdataF_comp[0][(start_symbol+1)*NR_NB_SC_PER_RB * pusch_pdu->rb_size],(nb_symb_sch-1)*NR_NB_SC_PER_RB * pusch_pdu->rb_size,1,1);
+#ifdef DEBUG_UL_PTRS
+          LOG_M("rxdataF_ptrs_comp.m","ptrs_comp",
+                &gNB->pusch_vars[0]->rxdataF_ptrs_comp[0][(start_symbol+1)*NR_NB_SC_PER_RB * pusch_pdu->rb_size],(nb_symb_sch-1)*NR_NB_SC_PER_RB * pusch_pdu->rb_size,1,1);
+#endif
 	}
         start_meas(&gNB->phy_proc_rx);
         ////////////////////////////////////////////////////////////
@@ -902,6 +911,7 @@ int main(int argc, char **argv)
       if (print_perf==1) {
         printDistribution(&gNB->phy_proc_rx,table_rx,"Total PHY proc rx");
         printStatIndent(&gNB->ulsch_channel_estimation_stats,"ULSCH channel estimation time");
+        printStatIndent(&gNB->ulsch_ptrs_processing_stats,"ULSCH PTRS Processing time");
         printStatIndent(&gNB->ulsch_rbs_extraction_stats,"ULSCH rbs extraction time");
         printStatIndent(&gNB->ulsch_channel_compensation_stats,"ULSCH channel compensation time");
         printStatIndent(&gNB->ulsch_llr_stats,"ULSCH llr computation");

openair1/SIMULATION/TOOLS/phase_noise.c

+/*
+ * Licensed to the OpenAirInterface (OAI) Software Alliance under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The OpenAirInterface Software Alliance licenses this file to You under
+ * the OAI Public License, Version 1.1  (the "License"); you may not use this file
+ * except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.openairinterface.org/?page_id=698
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *-------------------------------------------------------------------------------
+ * For more information about the OpenAirInterface (OAI) Software Alliance:
+ *      contact@openairinterface.org
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include  "sim.h"
+
+
+/* linear phase noise model */
+void phase_noise(double ts, int16_t * InRe, int16_t * InIm)
+{
+  double fd = 300;//0.01*30000
+  static double i=0;
+  double real, imag,x ,y;
+  i++;
+  real = cos(fd * 2 * M_PI * i * ts);
+  imag = sin (fd * 2 * M_PI * i * ts);
+  x = ((real * (double)InRe[0]) - (imag * (double)InIm[0])) ;
+  y= ((real * (double)InIm[0]) + (imag * (double)InRe[0])) ;
+  InRe[0]= (int16_t)(x);
+  InIm[0]= (int16_t)(y);
+}

openair1/SIMULATION/TOOLS/sim.h

@@ -440,6 +440,14 @@ void init_channelmod(void) ;
 double N_RB2sampling_rate(uint16_t N_RB);
 double N_RB2channel_bandwidth(uint16_t N_RB);
 
+/* Linear phase noise model */
+/*!
+  \brief This function produce phase noise and add to input signal
+  \param ts Sampling time 
+  \param *Re *Im Real and Imag part of the signal
+*/
+void phase_noise(double ts, int16_t * InRe, int16_t * InIm);
+
 #include "targets/RT/USER/rfsim.h"
 
 void do_DL_sig(sim_t *sim,