nr_ul_channel_estimation.c 36.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
/*
 * 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.0  (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 <string.h>

#include "nr_ul_estimation.h"
26
#include "PHY/sse_intrin.h"
27
#include "PHY/NR_REFSIG/nr_refsig.h"
28 29
#include "PHY/NR_REFSIG/ptrs_nr.h"
#include "PHY/NR_TRANSPORT/nr_transport_proto.h"
30 31
#include "PHY/NR_UE_ESTIMATION/filt16a_32.h"

Francesco Mani's avatar
Francesco Mani committed
32 33
//#define DEBUG_CH
//#define DEBUG_PUSCH
34

35 36
#define dBc(x,y) (dB_fixed(((int32_t)(x))*(x) + ((int32_t)(y))*(y)))

37
int nr_pusch_channel_estimation(PHY_VARS_gNB *gNB,
38 39 40
                                unsigned char Ns,
                                unsigned short p,
                                unsigned char symbol,
41
                                int ul_id,
42
                                unsigned short bwp_start_subcarrier,
43 44
                                nfapi_nr_pusch_pdu_t *pusch_pdu) {

45
  int pilot[3280] __attribute__((aligned(16)));
46 47
  unsigned char aarx;
  unsigned short k;
48 49
  unsigned int pilot_cnt,re_cnt;
  int16_t ch[2],ch_r[2],ch_l[2],*pil,*rxF,*ul_ch;
50
  int16_t *fl,*fm,*fr,*fml,*fmr,*fmm,*fdcl,*fdcr,*fdclh,*fdcrh;
51 52
  int ch_offset,symbol_offset ;
  int32_t **ul_ch_estimates_time =  gNB->pusch_vars[ul_id]->ul_ch_estimates_time;
53
  __m128i *ul_ch_128;
54 55 56 57 58 59 60 61 62

#ifdef DEBUG_CH
  FILE *debug_ch_est;
  debug_ch_est = fopen("debug_ch_est.txt","w");
#endif

  //uint16_t Nid_cell = (eNB_offset == 0) ? gNB->frame_parms.Nid_cell : gNB->measurements.adj_cell_id[eNB_offset-1];

  uint8_t nushift;
63
  int **ul_ch_estimates  = gNB->pusch_vars[ul_id]->ul_ch_estimates;
64 65 66 67 68 69 70 71 72 73 74 75
  int **rxdataF = gNB->common_vars.rxdataF;

  nushift = (p>>1)&1;
  gNB->frame_parms.nushift = nushift;

  ch_offset     = gNB->frame_parms.ofdm_symbol_size*symbol;

  symbol_offset = gNB->frame_parms.ofdm_symbol_size*symbol;

  k = bwp_start_subcarrier;
  int re_offset = k;

76 77
  uint16_t nb_rb_pusch = pusch_pdu->rb_size;

78 79
/*
#ifdef DEBUG_CH
80
  printf("PUSCH Channel Estimation : ch_offset %d, symbol_offset %d OFDM size %d, Ncp=%d, l=%d, Ns=%d, k=%d symbol %d\n", ,ch_offset,symbol_offset,gNB->frame_parms.ofdm_symbol_size,
81 82 83 84 85 86 87 88 89 90 91
         gNB->frame_parms.Ncp,l,Ns,k, symbol);
#endif
*/
  switch (nushift) {
   case 0:
         fl = filt8_l0;
         fm = filt8_m0;
         fr = filt8_r0;
         fmm = filt8_mm0;
         fml = filt8_m0;
         fmr = filt8_mr0;
92 93 94 95
         fdcl = filt8_dcl0;
         fdcr = filt8_dcr0;
         fdclh = filt8_dcl0_h;
         fdcrh = filt8_dcr0_h;
96 97 98 99 100 101 102 103 104
         break;

   case 1:
         fl = filt8_l1;
         fm = filt8_m1;
         fr = filt8_r1;
         fmm = filt8_mm1;
         fml = filt8_ml1;
         fmr = filt8_m1;
105 106 107 108
         fdcl = filt8_dcl1;
         fdcr = filt8_dcr1;
         fdclh = filt8_dcl1_h;
         fdcrh = filt8_dcr1_h;
109 110 111
         break;

   default:
112 113 114 115 116
#ifdef DEBUG_CH
      if (debug_ch_est)
        fclose(debug_ch_est);

#endif
117 118 119 120 121 122 123
     return(-1);
     break;
   }


  //------------------generate DMRS------------------//

124
  if (pusch_pdu->transform_precoding==1) // if transform precoding is disabled
Raymond Knopp's avatar
Raymond Knopp committed
125
    nr_pusch_dmrs_rx(gNB, Ns, gNB->nr_gold_pusch_dmrs[pusch_pdu->scid][Ns][symbol], &pilot[0], 1000, 0, nb_rb_pusch, pusch_pdu->rb_start*NR_NB_SC_PER_RB, pusch_pdu->dmrs_config_type);
126 127
  else
    nr_pusch_dmrs_rx(gNB, Ns, gNB->nr_gold_pusch_dmrs[pusch_pdu->scid][Ns][symbol], &pilot[0], 1000, 0, nb_rb_pusch, 0, pusch_pdu->dmrs_config_type);
128 129

  //------------------------------------------------//
130 131 132 133
#ifdef DEBUG_PUSCH
  for (int i=0;i<(6*nb_rb_pusch);i++)
    printf("%d+j*(%d)\n",((int16_t*)pilot)[2*i],((int16_t*)pilot)[1+(2*i)]);
#endif
134 135
  for (aarx=0; aarx<gNB->frame_parms.nb_antennas_rx; aarx++) {

136 137 138
    pil   = (int16_t *)&pilot[0];
    rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+k+nushift)];
    ul_ch = (int16_t *)&ul_ch_estimates[aarx][ch_offset];
139

140
    memset(ul_ch,0,4*(gNB->frame_parms.ofdm_symbol_size));
141 142

#ifdef DEBUG_PUSCH
143
    printf("symbol_offset %d, nushift %d\n",symbol_offset,nushift);
144
    printf("ch est pilot addr %p RB_DL %d\n",&pilot[0], gNB->frame_parms.N_RB_UL);
145
    printf("bwp_start_subcarrier %d, k %d, first_carrier %d, nb_rb_pusch %d\n",bwp_start_subcarrier,k,gNB->frame_parms.first_carrier_offset,nb_rb_pusch);
146 147 148
    printf("rxF addr %p p %d\n", rxF,p);
    printf("ul_ch addr %p nushift %d\n",ul_ch,nushift);
#endif
149
    //if ((gNB->frame_parms.N_RB_UL&1)==0) {
150

151
    if (pusch_pdu->dmrs_config_type == pusch_dmrs_type1){
152

153 154 155 156 157 158
      // Treat first 2 pilots specially (left edge)
      ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
      ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);

#ifdef DEBUG_PUSCH
      printf("ch 0 %d\n",((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1]));
159 160
      printf("pilot 0 : rxF - > (%d,%d) (%d)  ch -> (%d,%d) (%d), pil -> (%d,%d) \n",rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
      printf("data 0 : rxF - > (%d,%d) (%d)\n",rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
161
#endif
162

163 164 165 166 167
      multadd_real_vector_complex_scalar(fl,
                                         ch,
                                         ul_ch,
                                         8);
      pil+=2;
168
      re_offset = (re_offset+2) % gNB->frame_parms.ofdm_symbol_size;
169 170 171 172 173 174 175 176
      rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
      //for (int i= 0; i<8; i++)
      //printf("ul_ch addr %p %d\n", ul_ch+i, *(ul_ch+i));

      ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
      ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);

#ifdef DEBUG_PUSCH
177 178
      printf("pilot 1 : rxF - > (%d,%d) (%d) ch -> (%d,%d) (%d), pil -> (%d,%d) \n",rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
      printf("data 1 : rxF - > (%d,%d) (%d)\n",rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
179 180 181 182 183 184
#endif
      multadd_real_vector_complex_scalar(fml,
                                         ch,
                                         ul_ch,
                                         8);
      pil+=2;
185
      re_offset = (re_offset+2) % gNB->frame_parms.ofdm_symbol_size;
186 187 188 189 190 191 192
      rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
      //printf("ul_ch addr %p\n",ul_ch);
      
      ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
      ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);

#ifdef DEBUG_PUSCH
193 194
      printf("pilot 2 : rxF - > (%d,%d) (%d) ch -> (%d,%d) (%d), pil -> (%d,%d) \n",rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
      printf("data 2 : rxF - > (%d,%d) (%d)\n",rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
195 196 197 198 199 200 201 202 203 204
#endif
      multadd_real_vector_complex_scalar(fmm,
                                         ch,
                                         ul_ch,
                                         8);
                                         
      //for (int i= 0; i<16; i++)
      //printf("ul_ch addr %p %d\n", ul_ch+i, *(ul_ch+i));
      
      pil+=2;
205
      re_offset = (re_offset+2) % gNB->frame_parms.ofdm_symbol_size;
206 207 208 209 210 211 212
      rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
      ul_ch+=8;

      for (pilot_cnt=3; pilot_cnt<(6*nb_rb_pusch-3); pilot_cnt+=2) {

        ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
        ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
213 214

  #ifdef DEBUG_PUSCH
215 216
        printf("pilot %u : rxF - > (%d,%d) (%d) ch -> (%d,%d) (%d), pil -> (%d,%d) \n",pilot_cnt,rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
	printf("data %u : rxF - > (%d,%d) (%d)\n",pilot_cnt,rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
217 218
  #endif
        multadd_real_vector_complex_scalar(fml,
219 220 221 222
                                           ch,
                                           ul_ch,
                                           8);
        pil+=2;
223
        re_offset = (re_offset+2) % gNB->frame_parms.ofdm_symbol_size;
224
        rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
225 226
        //printf("ul_ch addr %p\n",ul_ch);

227 228
        ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
        ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
229 230

  #ifdef DEBUG_PUSCH
231 232
        printf("pilot %u : rxF - > (%d,%d) (%d) ch -> (%d,%d) (%d), pil -> (%d,%d) \n",pilot_cnt+1,rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
	printf("data %u : rxF - > (%d,%d) (%d)\n",pilot_cnt+1,rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
233
  #endif
234 235 236 237
        multadd_real_vector_complex_scalar(fmm,
                                           ch,
                                           ul_ch,
                                           8);
238 239 240 241

        //for (int i= 0; i<16; i++)
        //printf("ul_ch addr %p %d\n", ul_ch+i, *(ul_ch+i));

242
        pil+=2;
243
        re_offset = (re_offset+2) % gNB->frame_parms.ofdm_symbol_size;
244 245 246 247 248 249
        rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
        ul_ch+=8;

      }
      
      // Treat first 2 pilots specially (right edge)
250
      ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
251 252
      ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
#ifdef DEBUG_PUSCH
253 254
      printf("pilot %u : rxF - > (%d,%d) (%d) ch -> (%d,%d) (%d), pil -> (%d,%d) \n",pilot_cnt,rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
      printf("data %u : rxF - > (%d,%d) (%d)\n",pilot_cnt,rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
255 256 257 258 259 260 261 262 263 264
#endif
      multadd_real_vector_complex_scalar(fm,
                                         ch,
                                         ul_ch,
                                         8);
                                         
      //for (int i= 0; i<8; i++)
      //printf("ul_ch addr %p %d\n", ul_ch+i, *(ul_ch+i));

      pil+=2;
265
      re_offset = (re_offset+2) % gNB->frame_parms.ofdm_symbol_size;
266 267 268 269 270 271
      rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
             
      ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
      ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
#ifdef DEBUG_PUSCH
      printf("ch 0 %d\n",((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1]));
272 273
      printf("pilot %u : rxF - > (%d,%d) (%d) ch -> (%d,%d) (%d), pil -> (%d,%d) \n",pilot_cnt+1,rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
      printf("data %u : rxF - > (%d,%d) (%d)\n",pilot_cnt+1,rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
274 275 276 277 278 279 280
#endif
      multadd_real_vector_complex_scalar(fmr,
                                         ch,
                                         ul_ch,
                                         8);
                                         
      pil+=2;
281
      re_offset = (re_offset+2) % gNB->frame_parms.ofdm_symbol_size;
282 283 284 285 286 287
      rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
      ul_ch+=8;
      
      ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
      ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
#ifdef DEBUG_PUSCH
288 289
      printf("pilot %u: rxF - > (%d,%d) (%d) ch -> (%d,%d) (%d), pil -> (%d,%d) \n",pilot_cnt+2,rxF[0],rxF[1],dBc(rxF[0],rxF[1]),ch[0],ch[1],dBc(ch[0],ch[1]),pil[0],pil[1]);
      printf("data %u : rxF - > (%d,%d) (%d)\n",pilot_cnt+2,rxF[2],rxF[3],dBc(rxF[2],rxF[3]));
290 291 292 293 294 295
#endif
      multadd_real_vector_complex_scalar(fr,
                                         ch,
                                         ul_ch,
                                         8);

296

297 298 299 300 301 302 303 304 305 306 307
      // check if PRB crosses DC and improve estimates around DC
      if ((bwp_start_subcarrier < gNB->frame_parms.ofdm_symbol_size) && (bwp_start_subcarrier+nb_rb_pusch*12 >= gNB->frame_parms.ofdm_symbol_size)) {
        ul_ch = (int16_t *)&ul_ch_estimates[aarx][ch_offset];
        uint16_t idxDC = 2*(gNB->frame_parms.ofdm_symbol_size - bwp_start_subcarrier);
        uint16_t idxPil = idxDC/2;
        re_offset = k;
        pil = (int16_t *)&pilot[0];
        pil += (idxPil-2);
        ul_ch += (idxDC-4);
        ul_ch = memset(ul_ch, 0, sizeof(int16_t)*10);
        re_offset = (re_offset+idxDC/2-2) % gNB->frame_parms.ofdm_symbol_size;
308 309 310
        rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
        ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
        ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
311 312 313 314 315 316 317 318

        // for proper allignment of SIMD vectors
        if((gNB->frame_parms.N_RB_UL&1)==0) {

          multadd_real_vector_complex_scalar(fdcl,
                                             ch,
                                             ul_ch-4,
                                             8);
319
        
320 321 322 323 324
          pil += 4;
          re_offset = (re_offset+4) % gNB->frame_parms.ofdm_symbol_size;
          rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
          ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
          ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
325
        
326 327 328 329 330 331 332 333 334 335
          multadd_real_vector_complex_scalar(fdcr,
                                             ch,
                                             ul_ch-4,
                                             8);
        }
        else {
          multadd_real_vector_complex_scalar(fdclh,
                                             ch,
                                             ul_ch,
                                             8);
336
        
337 338 339 340 341
          pil += 4;
          re_offset = (re_offset+4) % gNB->frame_parms.ofdm_symbol_size;
          rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];
          ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
          ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);
342
        
343 344 345 346 347
          multadd_real_vector_complex_scalar(fdcrh,
                                             ch,
                                             ul_ch,
                                             8);
        }
348
      }
349
#ifdef DEBUG_PUSCH
350 351 352 353 354 355
      ul_ch = (int16_t *)&ul_ch_estimates[aarx][ch_offset];
      for(uint16_t idxP=0; idxP<ceil((float)nb_rb_pusch*12/8); idxP++) {
        for(uint8_t idxI=0; idxI<16; idxI+=2) {
          printf("%d\t%d\t",ul_ch[idxP*16+idxI],ul_ch[idxP*16+idxI+1]);
        }
        printf("%d\n",idxP);
356 357
      }
#endif    
358 359
    }
    else { //pusch_dmrs_type2  |p_r,p_l,d,d,d,d,p_r,p_l,d,d,d,d|
360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442

      // Treat first DMRS specially (left edge)

        rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];

        ul_ch[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
        ul_ch[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);

        pil+=2;
        ul_ch+=2;
        re_offset = (re_offset + 1)%gNB->frame_parms.ofdm_symbol_size;
        ch_offset++;

        for (re_cnt = 1; re_cnt < (nb_rb_pusch*NR_NB_SC_PER_RB) - 5; re_cnt+=6){

          rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];

          ch_l[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
          ch_l[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);

          ul_ch[0] = ch_l[0];
          ul_ch[1] = ch_l[1];

          pil+=2;
          ul_ch+=2;
          ch_offset++;

          multadd_real_four_symbols_vector_complex_scalar(filt8_ml2,
                                                          ch_l,
                                                          ul_ch);

          re_offset = (re_offset+5)%gNB->frame_parms.ofdm_symbol_size;

          rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];

          ch_r[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
          ch_r[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);


          multadd_real_four_symbols_vector_complex_scalar(filt8_mr2,
                                                          ch_r,
                                                          ul_ch);

          //for (int re_idx = 0; re_idx < 8; re_idx+=2)
            //printf("ul_ch = %d + j*%d\n", ul_ch[re_idx], ul_ch[re_idx+1]);

          ul_ch+=8;
          ch_offset+=4;

          ul_ch[0] = ch_r[0];
          ul_ch[1] = ch_r[1];

          pil+=2;
          ul_ch+=2;
          ch_offset++;
          re_offset = (re_offset + 1)%gNB->frame_parms.ofdm_symbol_size;

        }

        // Treat last pilot specially (right edge)

        rxF   = (int16_t *)&rxdataF[aarx][(symbol_offset+nushift+re_offset)];

        ch_l[0] = (int16_t)(((int32_t)pil[0]*rxF[0] - (int32_t)pil[1]*rxF[1])>>15);
        ch_l[1] = (int16_t)(((int32_t)pil[0]*rxF[1] + (int32_t)pil[1]*rxF[0])>>15);

        ul_ch[0] = ch_l[0];
        ul_ch[1] = ch_l[1];

        ul_ch+=2;
        ch_offset++;

        multadd_real_four_symbols_vector_complex_scalar(filt8_rr1,
                                                        ch_l,
                                                        ul_ch);

        multadd_real_four_symbols_vector_complex_scalar(filt8_rr2,
                                                        ch_r,
                                                        ul_ch);

        ul_ch_128 = (__m128i *)&ul_ch_estimates[aarx][ch_offset];

        ul_ch_128[0] = _mm_slli_epi16 (ul_ch_128[0], 2);
443
    }
444 445


446
    // Convert to time domain
447

448
    switch (gNB->frame_parms.ofdm_symbol_size) {
449
        case 128:
frtabu's avatar
frtabu committed
450
          idft(IDFT_128,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
451 452 453 454 455
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        case 256:
frtabu's avatar
frtabu committed
456
          idft(IDFT_256,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
457 458 459 460 461
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        case 512:
frtabu's avatar
frtabu committed
462
          idft(IDFT_512,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
463 464 465 466 467
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        case 1024:
frtabu's avatar
frtabu committed
468
          idft(IDFT_1024,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
469 470 471 472 473
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        case 1536:
frtabu's avatar
frtabu committed
474
          idft(IDFT_1536,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
475 476 477 478 479
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        case 2048:
frtabu's avatar
frtabu committed
480
          idft(IDFT_2048,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
481 482 483 484 485
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        case 4096:
frtabu's avatar
frtabu committed
486
          idft(IDFT_4096,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
487 488 489 490 491
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        case 8192:
frtabu's avatar
frtabu committed
492
          idft(IDFT_8192,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
493 494 495 496 497
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;

        default:
frtabu's avatar
frtabu committed
498
          idft(IDFT_512,(int16_t*) &ul_ch_estimates[aarx][symbol_offset],
499 500 501 502 503
                 (int16_t*) ul_ch_estimates_time[aarx],
                 1);
          break;
      }

504 505 506 507 508 509 510
  }

#ifdef DEBUG_CH
  fclose(debug_ch_est);
#endif

  return(0);
511
}
512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567


/*******************************************************************
 *
 * 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;

  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];
568 569 570 571 572
    /* set the previous estimations to zero at first symbol */
    if(symbol == rel15_ul->start_symbol_index)
    {
      memset(phase_per_symbol,0,sizeof(int32_t)*frame_parms->symbols_per_slot);
    }
573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631
    /* 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);
    }
    /* 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 */

#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 );
632 633 634
      LOG_M("rxdataF_bf_ptrs_comp.m","bf_ptrs_cmp",
            &gNB->pusch_vars[0]->rxdataF_comp[aarx][rel15_ul->start_symbol_index * NR_NB_SC_PER_RB * rel15_ul->rb_size],
            rel15_ul->nr_of_symbols * NR_NB_SC_PER_RB * rel15_ul->rb_size,1,1);
635 636 637 638 639 640 641 642 643 644 645 646
#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],
647
                          (int16_t*)&gNB->pusch_vars[ulsch_id]->rxdataF_comp[aarx][(i * rel15_ul->rb_size * NR_NB_SC_PER_RB)],
648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685
                          (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;
686 687
  uint8_t               K_ptrs           = rel15_ul->pusch_ptrs.ptrs_freq_density;
  uint16_t              sc_per_symbol    = (rel15_ul->rb_size + K_ptrs - 1)/K_ptrs;
688 689 690 691 692 693 694 695 696
  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 */
697
  nr_gen_ref_conj_symbols(ptrs_gold_seq[Ns][symbol],sc_per_symbol*2,ptrs_p, NR_MOD_TABLE_QPSK_OFFSET,2);// 2 for QPSK
698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916

  /* 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;
}