nr_ul_channel_estimation.c 38.7 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
#include "PHY/NR_REFSIG/dmrs_nr.h"
29 30
#include "PHY/NR_REFSIG/ptrs_nr.h"
#include "PHY/NR_TRANSPORT/nr_transport_proto.h"
31
#include "PHY/NR_UE_TRANSPORT/srs_modulation_nr.h"
32 33
#include "PHY/NR_UE_ESTIMATION/filt16a_32.h"

dir's avatar
dir committed
34
#include "PHY/NR_REFSIG/ul_ref_seq_nr.h"
35
#include "executables/softmodem-common.h"
dir's avatar
dir committed
36 37


Francesco Mani's avatar
Francesco Mani committed
38 39
//#define DEBUG_CH
//#define DEBUG_PUSCH
40
//#define SRS_DEBUG
41

42
#define NO_INTERP 1
43 44
#define dBc(x,y) (dB_fixed(((int32_t)(x))*(x) + ((int32_t)(y))*(y)))

45 46 47 48 49 50 51 52 53 54 55
__attribute__((always_inline)) inline c16_t c32x16cumulVectVectWithSteps(c16_t *in1,
                                                                         int *offset1,
                                                                         const int step1,
                                                                         c16_t *in2,
                                                                         int *offset2,
                                                                         const int step2,
                                                                         const int modulo2,
                                                                         const int N) {

  int localOffset1=*offset1;
  int localOffset2=*offset2;
56
  c32_t cumul={0};
57 58 59 60 61 62 63 64 65 66
  for (int i=0; i<N; i++) {
    cumul=c32x16maddShift(in1[localOffset1], in2[localOffset2], cumul, 15);
    localOffset1+=step1;
    localOffset2= (localOffset2 + step2) % modulo2;
  }
  *offset1=localOffset1;
  *offset2=localOffset2;
  return c16x32div(cumul, N);
}

67
int nr_pusch_channel_estimation(PHY_VARS_gNB *gNB,
68 69 70
                                unsigned char Ns,
                                unsigned short p,
                                unsigned char symbol,
71
                                int ul_id,
72
                                unsigned short bwp_start_subcarrier,
73
                                nfapi_nr_pusch_pdu_t *pusch_pdu,
74 75
                                int *max_ch,
                                uint32_t *nvar) {
76

77
  c16_t pilot[3280] __attribute__((aligned(32)));
78
  const int chest_freq = gNB->chest_freq;
79 80 81 82 83 84

#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];
francescomani's avatar
francescomani committed
85 86
  NR_gNB_PUSCH *pusch_vars = &gNB->pusch_vars[ul_id];
  c16_t **ul_ch_estimates = (c16_t **)pusch_vars->ul_ch_estimates;
87 88
  const int symbolSize = gNB->frame_parms.ofdm_symbol_size;
  const int soffset = (Ns&3)*gNB->frame_parms.symbols_per_slot*symbolSize;
89
  const int nushift = (p >> 1) & 1;
90 91
  int ch_offset     = symbolSize*symbol;
  const int symbol_offset = symbolSize*symbol;
92

93 94
  const int k0 = bwp_start_subcarrier;
  const int nb_rb_pusch = pusch_pdu->rb_size;
95

96
  LOG_D(PHY, "In %s: ch_offset %d, soffset %d, symbol_offset %d, OFDM size %d, Ns = %d, k0 = %d, symbol %d\n",
cig's avatar
cig committed
97
        __FUNCTION__,
Eurecom's avatar
Eurecom committed
98
        ch_offset, soffset,
cig's avatar
cig committed
99
        symbol_offset,
100
        symbolSize,
cig's avatar
cig committed
101
        Ns,
102
        k0,
cig's avatar
cig committed
103 104
        symbol);

105 106
  //------------------generate DMRS------------------//

107 108 109 110
  if(pusch_pdu->ul_dmrs_scrambling_id != gNB->pusch_gold_init[pusch_pdu->scid])  {
    gNB->pusch_gold_init[pusch_pdu->scid] = pusch_pdu->ul_dmrs_scrambling_id;
    nr_gold_pusch(gNB, pusch_pdu->scid, pusch_pdu->ul_dmrs_scrambling_id);
  }
111

112
  if (pusch_pdu->transform_precoding == transformPrecoder_disabled) {
113 114 115 116
    // Note: pilot returned by the following function is already the complex conjugate of the transmitted DMRS
    nr_pusch_dmrs_rx(gNB,
                     Ns,
                     gNB->nr_gold_pusch_dmrs[pusch_pdu->scid][Ns][symbol],
117
                     pilot,
118 119 120 121 122
                     (1000 + p),
                     0,
                     nb_rb_pusch,
                     (pusch_pdu->bwp_start + pusch_pdu->rb_start) * NR_NB_SC_PER_RB,
                     pusch_pdu->dmrs_config_type);
123
  } else { // if transform precoding or SC-FDMA is enabled in Uplink
dir's avatar
dir committed
124
    // NR_SC_FDMA supports type1 DMRS so only 6 DMRS REs per RB possible
125
    const int index = get_index_for_dmrs_lowpapr_seq(nb_rb_pusch * (NR_NB_SC_PER_RB / 2));
126 127
    const uint8_t u = pusch_pdu->dfts_ofdm.low_papr_group_number;
    const uint8_t v = pusch_pdu->dfts_ofdm.low_papr_sequence_number;
128
    c16_t *dmrs_seq = gNB_dmrs_lowpaprtype1_sequence[u][v][index];
129
    LOG_D(PHY,"Transform Precoding params. u: %d, v: %d, index for dmrsseq: %d\n", u, v, index);
dir's avatar
dir committed
130 131
    AssertFatal(index >= 0, "Num RBs not configured according to 3GPP 38.211 section 6.3.1.4. For PUSCH with transform precoding, num RBs cannot be multiple of any other primenumber other than 2,3,5\n");
    AssertFatal(dmrs_seq != NULL, "DMRS low PAPR seq not found, check if DMRS sequences are generated");
132
    nr_pusch_lowpaprtype1_dmrs_rx(gNB, Ns, dmrs_seq, pilot, 1000, 0, nb_rb_pusch, 0, pusch_pdu->dmrs_config_type);
133 134 135 136
#ifdef DEBUG_PUSCH
    printf ("NR_UL_CHANNEL_EST: index %d, u %d,v %d\n", index, u, v);
    LOG_M("gNb_DMRS_SEQ.m","gNb_DMRS_SEQ", dmrs_seq,6*nb_rb_pusch,1,1);
#endif
dir's avatar
dir committed
137
  }
138
  //------------------------------------------------//
dir's avatar
dir committed
139

140
#ifdef DEBUG_PUSCH
141

cig's avatar
cig committed
142
  for (int i = 0; i < (6 * nb_rb_pusch); i++) {
143
    LOG_I(PHY, "In %s: %d + j*(%d)\n", __FUNCTION__, pilot[i].r,pilot[i].i);
cig's avatar
cig committed
144 145
  }

146
#endif
147

148 149
  int nest_count = 0;
  uint64_t noise_amp2 = 0;
150
  c16_t ul_ls_est[symbolSize] __attribute__((aligned(32)));
151
  memset(ul_ls_est, 0, sizeof(c16_t) * symbolSize);
152
  delay_t *delay = &gNB->ulsch[ul_id].delay;
francescomani's avatar
francescomani committed
153
  memset(delay, 0, sizeof(*delay));
154

155 156 157
  for (int aarx=0; aarx<gNB->frame_parms.nb_antennas_rx; aarx++) {
    c16_t *rxdataF = (c16_t *)&gNB->common_vars.rxdataF[aarx][symbol_offset];
    c16_t *ul_ch = &ul_ch_estimates[p*gNB->frame_parms.nb_antennas_rx+aarx][ch_offset];
158

159
    memset(ul_ch,0,sizeof(*ul_ch)*symbolSize);
160
#ifdef DEBUG_PUSCH
cig's avatar
cig committed
161
    LOG_I(PHY, "In %s symbol_offset %d, nushift %d\n", __FUNCTION__, symbol_offset, nushift);
162
    LOG_I(PHY, "In %s ch est pilot, N_RB_UL %d\n", __FUNCTION__, gNB->frame_parms.N_RB_UL);
163
    LOG_I(PHY, "In %s bwp_start_subcarrier %d, k0 %d, first_carrier %d, nb_rb_pusch %d\n", __FUNCTION__, bwp_start_subcarrier, k0, gNB->frame_parms.first_carrier_offset, nb_rb_pusch);
cig's avatar
cig committed
164
    LOG_I(PHY, "In %s ul_ch addr %p nushift %d\n", __FUNCTION__, ul_ch, nushift);
165 166
#endif

167
    if (pusch_pdu->dmrs_config_type == pusch_dmrs_type1 && chest_freq == 0) {
168
      c16_t *pil   = pilot;
169
      int re_offset = k0;
170
      LOG_D(PHY,"PUSCH estimation DMRS type 1, Freq-domain interpolation");
171 172
      // For configuration type 1: k = 4*n + 2*k' + delta,
      // where k' is 0 or 1, and delta is in Table 6.4.1.1.3-1 from TS 38.211
173
      int pilot_cnt = 0;
174
      int delta = nr_pusch_dmrs_delta(pusch_dmrs_type1, p);
175

176
      for (int n = 0; n < 3 * nb_rb_pusch; n++) {
177
        // LS estimation
178
        c32_t ch = {0};
179

180
        for (int k_line = 0; k_line <= 1; k_line++) {
181
          re_offset = (k0 + (n << 2) + (k_line << 1) + delta) % symbolSize;
182
          ch = c32x16maddShift(*pil, rxdataF[soffset + re_offset], ch, 16);
183
          pil++;
184
        }
185

186
        c16_t ch16 = {.r = (int16_t)ch.r, .i = (int16_t)ch.i};
187
        *max_ch = max(*max_ch, max(abs(ch.r), abs(ch.i)));
188 189 190
        for (int k = pilot_cnt << 1; k < (pilot_cnt << 1) + 4; k++) {
          ul_ls_est[k] = ch16;
        }
191 192 193
        pilot_cnt += 2;
      }

194
      nr_est_delay(gNB->frame_parms.ofdm_symbol_size, ul_ls_est, (c16_t *)pusch_vars->ul_ch_estimates_time[aarx], delay);
195 196
      int delay_idx = get_delay_idx(delay->est_delay, MAX_DELAY_COMP);
      c16_t *ul_delay_table = gNB->frame_parms.delay_table[delay_idx];
197 198

#ifdef DEBUG_PUSCH
199
      printf("Estimated delay = %i\n", delay->est_delay >> 1);
200 201 202 203 204
#endif

      pilot_cnt = 0;
      for (int n = 0; n < 3*nb_rb_pusch; n++) {

205 206
        // Channel interpolation
        for (int k_line = 0; k_line <= 1; k_line++) {
207 208 209 210 211

          // Apply delay
          int k = pilot_cnt << 1;
          c16_t ch16 = c16mulShift(ul_ls_est[k], ul_delay_table[k], 8);

212
#ifdef DEBUG_PUSCH
213 214
          re_offset = (k0 + (n << 2) + (k_line << 1)) % symbolSize;
          c16_t *rxF = &rxdataF[soffset + re_offset];
215 216
          printf("pilot %4d: pil -> (%6d,%6d), rxF -> (%4d,%4d), ch -> (%4d,%4d)\n",
                 pilot_cnt, pil->r, pil->i, rxF->r, rxF->i, ch.r, ch.i);
217
#endif
218

219
          if (pilot_cnt == 0) {
220
            c16multaddVectRealComplex(filt16_ul_p0, &ch16, ul_ch, 16);
221
          } else if (pilot_cnt == 1 || pilot_cnt == 2) {
222
            c16multaddVectRealComplex(filt16_ul_p1p2, &ch16, ul_ch, 16);
223
          } else if (pilot_cnt == (6 * nb_rb_pusch - 1)) {
224
            c16multaddVectRealComplex(filt16_ul_last, &ch16, ul_ch, 16);
225
          } else {
226
            c16multaddVectRealComplex(filt16_ul_middle, &ch16, ul_ch, 16);
227 228 229
            if (pilot_cnt % 2 == 0) {
              ul_ch += 4;
            }
230
          }
231

232
          pilot_cnt++;
233
        }
234
      }
235

236 237 238
      // Revert delay
      pilot_cnt = 0;
      ul_ch = &ul_ch_estimates[p * gNB->frame_parms.nb_antennas_rx + aarx][ch_offset];
239 240
      int inv_delay_idx = get_delay_idx(-delay->est_delay, MAX_DELAY_COMP);
      c16_t *ul_inv_delay_table = gNB->frame_parms.delay_table[inv_delay_idx];
241 242 243 244 245
      for (int n = 0; n < 3 * nb_rb_pusch; n++) {
        for (int k_line = 0; k_line <= 1; k_line++) {
          int k = pilot_cnt << 1;
          ul_ch[k] = c16mulShift(ul_ch[k], ul_inv_delay_table[k], 8);
          ul_ch[k + 1] = c16mulShift(ul_ch[k + 1], ul_inv_delay_table[k + 1], 8);
246 247 248
          noise_amp2 += c16amp2(c16sub(ul_ls_est[k], ul_ch[k]));
          noise_amp2 += c16amp2(c16sub(ul_ls_est[k + 1], ul_ch[k + 1]));

249
#ifdef DEBUG_PUSCH
250 251 252 253 254
          re_offset = (k0 + (n << 2) + (k_line << 1)) % symbolSize;
          c16_t *rxF = &rxdataF[soffset + re_offset];
          printf("ch -> (%4d,%4d), ch_inter -> (%4d,%4d)\n", ul_ls_est[k].r, ul_ls_est[k].i, ul_ch[k].r, ul_ch[k].i);
#endif
          pilot_cnt++;
255
          nest_count += 2;
256
        }
257
      }
258

259 260 261 262 263 264
    } else if (pusch_pdu->dmrs_config_type == pusch_dmrs_type2 && chest_freq == 0) { // pusch_dmrs_type2  |p_r,p_l,d,d,d,d,p_r,p_l,d,d,d,d|
      LOG_D(PHY, "PUSCH estimation DMRS type 2, Freq-domain interpolation\n");
      c16_t *pil = pilot;
      c16_t *rx = &rxdataF[soffset + nushift];
      for (int n = 0; n < nb_rb_pusch * NR_NB_SC_PER_RB; n += 6) {
        c16_t ch0 = c16mulShift(*pil, rx[(k0 + n) % symbolSize], 15);
265
        pil++;
266
        c16_t ch1 = c16mulShift(*pil, rx[(k0 + n + 1) % symbolSize], 15);
267
        pil++;
268 269
        c16_t ch = c16addShift(ch0, ch1, 1);
        *max_ch = max(*max_ch, max(abs(ch.r), abs(ch.i)));
270 271 272
        multadd_real_four_symbols_vector_complex_scalar(filt8_rep4, &ch, &ul_ls_est[n]);
        ul_ls_est[n + 4] = ch;
        ul_ls_est[n + 5] = ch;
273 274
        noise_amp2 += c16amp2(c16sub(ch0, ch));
        nest_count++;
275
      }
276 277

      // Delay compensation
278
      nr_est_delay(gNB->frame_parms.ofdm_symbol_size, ul_ls_est, (c16_t *)pusch_vars->ul_ch_estimates_time[aarx], delay);
279 280
      int delay_idx = get_delay_idx(-delay->est_delay, MAX_DELAY_COMP);
      c16_t *ul_delay_table = gNB->frame_parms.delay_table[delay_idx];
281 282 283 284
      for (int n = 0; n < nb_rb_pusch * NR_NB_SC_PER_RB; n++) {
        ul_ch[n] = c16mulShift(ul_ls_est[n], ul_delay_table[n % 6], 8);
      }

285
    }
286

287 288 289 290 291 292
    else if (pusch_pdu->dmrs_config_type == pusch_dmrs_type1) { // this is case without frequency-domain linear interpolation, just take average of LS channel estimates of 6 DMRS REs and use a common value for the whole PRB
      LOG_D(PHY,"PUSCH estimation DMRS type 1, no Freq-domain interpolation\n");
      c16_t *rxF   =  &rxdataF[soffset + nushift];
      int pil_offset = 0;
      int re_offset = k0;
      c16_t ch;
293

294
      // First PRB
295
      ch=c32x16cumulVectVectWithSteps(pilot, &pil_offset, 1, rxF, &re_offset, 2, symbolSize, 6);
296

297
#if NO_INTERP
298 299
      for (c16_t *end=ul_ch+12; ul_ch<end; ul_ch++)
        *ul_ch=ch;
300
#else
301
      c16multaddVectRealComplex(filt8_avlip0, &ch, ul_ch, 8);
302
      ul_ch += 8;
303
      c16multaddVectRealComplex(filt8_avlip1, &ch, ul_ch, 8);
304
      ul_ch += 8;
305
      c16multaddVectRealComplex(filt8_avlip2, &ch, ul_ch, 8);
306
      ul_ch -= 12;
307
#endif
308

309 310
      for (int pilot_cnt=6; pilot_cnt<6*(nb_rb_pusch-1); pilot_cnt += 6) {
        ch=c32x16cumulVectVectWithSteps(pilot, &pil_offset, 1, rxF, &re_offset, 2, symbolSize, 6);
311
        *max_ch = max(*max_ch, max(abs(ch.r), abs(ch.i)));
312

313
#if NO_INTERP
314 315
      for (c16_t *end=ul_ch+12; ul_ch<end; ul_ch++)
          *ul_ch=ch;
316
#else
317 318 319 320
        ul_ch[3].r += (ch.r * 1365)>>15; // 1/12*16384
        ul_ch[3].i += (ch.i * 1365)>>15; // 1/12*16384

        ul_ch += 4;
321
        c16multaddVectRealComplex(filt8_avlip3, &ch, ul_ch, 8);
322
        ul_ch += 8;
323
        c16multaddVectRealComplex(filt8_avlip4, &ch, ul_ch, 8);
324
        ul_ch += 8;
325
        c16multaddVectRealComplex(filt8_avlip5, &ch, ul_ch, 8);
326
        ul_ch -= 8;
327
#endif
328
      }
329
      // Last PRB
330
      ch=c32x16cumulVectVectWithSteps(pilot, &pil_offset, 1, rxF, &re_offset, 2, symbolSize, 6);
331

332
#if NO_INTERP
333 334
      for (c16_t *end=ul_ch+12; ul_ch<end; ul_ch++)
        *ul_ch=ch;
335
#else
336 337
      ul_ch[3].r += (ch.r * 1365)>>15; // 1/12*16384
      ul_ch[3].i += (ch.i * 1365)>>15; // 1/12*16384
338

339
      ul_ch += 4;
340
      c16multaddVectRealComplex(filt8_avlip3, &ch, ul_ch, 8);
341
      ul_ch += 8;
342
      c16multaddVectRealComplex(filt8_avlip6, &ch, ul_ch, 8);
343
#endif
344
    } else  { // this is case without frequency-domain linear interpolation, just take average of LS channel estimates of 4 DMRS REs and use a common value for the whole PRB
345
      LOG_D(PHY,"PUSCH estimation DMRS type 2, no Freq-domain interpolation");
346 347 348
      c16_t *pil   = pilot;
      int re_offset = k0;
      c32_t ch0={0};
349
      //First PRB
350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
      ch0=c32x16mulShift(*pil,
                         rxdataF[soffset + nushift + re_offset],
                         15);
      pil++;
      re_offset = (re_offset+1) % symbolSize;
      ch0=c32x16maddShift(*pil,
                          rxdataF[nushift+re_offset],
                          ch0,
                          15);
      pil++;
      re_offset = (re_offset+5) % symbolSize;
      ch0=c32x16maddShift(*pil,
                          rxdataF[nushift+re_offset],
                          ch0,
                          15);
      re_offset = (re_offset+1) % symbolSize;
      ch0=c32x16maddShift(*pil,
                          rxdataF[nushift+re_offset],
                          ch0,
                          15);
      pil++;
      re_offset = (re_offset+5) % symbolSize;

      c16_t ch=c16x32div(ch0, 4);
374
#if NO_INTERP
375 376
      for (c16_t *end=ul_ch+12; ul_ch<end; ul_ch++)
        *ul_ch=ch;
377
#else
378 379 380 381 382 383
      c16multaddVectRealComplex(filt8_avlip0, &ch, ul_ch, 8);
      ul_ch += 8;
      c16multaddVectRealComplex(filt8_avlip1, &ch, ul_ch, 8);
      ul_ch += 8;
      c16multaddVectRealComplex(filt8_avlip2, &ch, ul_ch, 8);
      ul_ch -= 12;
384
#endif
385

386 387 388 389 390
      for (int pilot_cnt=4; pilot_cnt<4*(nb_rb_pusch-1); pilot_cnt += 4) {
        c32_t ch0;
        ch0=c32x16mulShift(*pil, rxdataF[nushift+re_offset], 15);
        pil++;
        re_offset = (re_offset+1) % symbolSize;
391

392 393 394
        ch0=c32x16maddShift(*pil, rxdataF[nushift+re_offset], ch0, 15);
        pil++;
        re_offset = (re_offset+5) % symbolSize;
395

396 397 398
        ch0=c32x16maddShift(*pil, rxdataF[nushift+re_offset], ch0, 15);
        pil++;
        re_offset = (re_offset+1) % symbolSize;
399

400 401 402
        ch0=c32x16maddShift(*pil, rxdataF[nushift+re_offset], ch0, 15);
        pil++;
        re_offset = (re_offset+5) % symbolSize;
403

404
        ch=c16x32div(ch0, 4);
405
        *max_ch = max(*max_ch, max(abs(ch.r), abs(ch.i)));
406

407
#if NO_INTERP
408 409
        for (c16_t *end=ul_ch+12; ul_ch<end; ul_ch++)
          *ul_ch=ch;
410
#else
411
        ul_ch[3] = c16maddShift(ch, (c16_t){1365, 1365}, (c16_t){0, 0}, 15); // 1365 = 1/12*16384 (full range is +/- 32768)
412 413
        ul_ch += 4;
        c16multaddVectRealComplex(filt8_avlip3, &ch, ul_ch, 8);
414
        ul_ch += 8;
415
        c16multaddVectRealComplex(filt8_avlip4, &ch, ul_ch, 8);
416
        ul_ch += 8;
417 418
        c16multaddVectRealComplex(filt8_avlip5, &ch, ul_ch, 8);
        ul_ch -= 8;
419
#endif
420
      }
421

422 423 424 425
      // Last PRB
      ch0=c32x16mulShift(*pil, rxdataF[nushift+re_offset], 15);
      pil++;
      re_offset = (re_offset+1) % symbolSize;
426

427 428 429
      ch0=c32x16maddShift(*pil, rxdataF[nushift+re_offset], ch0, 15);
      pil++;
      re_offset = (re_offset+5) % symbolSize;
430

431 432 433
      ch0=c32x16maddShift(*pil, rxdataF[nushift+re_offset], ch0, 15);
      pil++;
      re_offset = (re_offset+1) % symbolSize;
434

435 436 437
      ch0=c32x16maddShift(*pil, rxdataF[nushift+re_offset], ch0, 15);
      pil++;
      re_offset = (re_offset+5) % symbolSize;
438

439
      ch=c16x32div(ch0, 4);
440
#if NO_INTERP
441 442
      for (c16_t *end=ul_ch+12; ul_ch<end; ul_ch++)
          *ul_ch=ch;
443
#else
444
      ul_ch[3] = c16maddShift(ch, (c16_t){1365, 1365}, (c16_t){0, 0}, 15); // 1365 = 1/12*16384 (full range is +/- 32768)
445 446
      ul_ch += 4;
      c16multaddVectRealComplex(filt8_avlip3, &ch, ul_ch, 8);
447
      ul_ch += 8;
448
      c16multaddVectRealComplex(filt8_avlip6, &ch, ul_ch, 8);
449
#endif
450
    }
451

452
#ifdef DEBUG_PUSCH
453 454 455 456
    ul_ch = &ul_ch_estimates[p * gNB->frame_parms.nb_antennas_rx + aarx][ch_offset];
    for (int idxP = 0; idxP < ceil((float)nb_rb_pusch * 12 / 8); idxP++) {
      for (int idxI = 0; idxI < 8; idxI++) {
          printf("%d\t%d\t", ul_ch[idxP * 8 + idxI].r, ul_ch[idxP * 8 + idxI].i);
457
      }
458
      printf("%d\n", idxP);
459
    }
460
#endif
461

462 463 464 465 466
  }

#ifdef DEBUG_CH
  fclose(debug_ch_est);
#endif
467 468 469 470 471 472

  if (nvar && nest_count > 0) {
    *nvar = (uint32_t)(noise_amp2 / nest_count);
  }

  return 0;
473
}
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496


/*******************************************************************
 *
 * 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
 *********************************************************************/
Tsung-Yu Chan's avatar
Tsung-Yu Chan committed
497
// #define DEBUG_UL_PTRS
498
void nr_pusch_ptrs_processing(PHY_VARS_gNB *gNB,
499
                              NR_DL_FRAME_PARMS *frame_parms,
500 501 502 503
                              nfapi_nr_pusch_pdu_t *rel15_ul,
                              uint8_t ulsch_id,
                              uint8_t nr_tti_rx,
                              unsigned char symbol,
francescomani's avatar
francescomani committed
504 505 506
                              uint32_t nb_re_pusch)
{
  NR_gNB_PUSCH *pusch_vars = &gNB->pusch_vars[ulsch_id];
507 508 509 510 511 512 513 514
  int32_t *ptrs_re_symbol   = NULL;
  int8_t   ret = 0;
  uint8_t  symbInSlot       = rel15_ul->start_symbol_index + rel15_ul->nr_of_symbols;
  uint8_t *startSymbIndex   = &rel15_ul->start_symbol_index;
  uint8_t *nbSymb           = &rel15_ul->nr_of_symbols;
  uint8_t  *L_ptrs          = &rel15_ul->pusch_ptrs.ptrs_time_density;
  uint8_t  *K_ptrs          = &rel15_ul->pusch_ptrs.ptrs_freq_density;
  uint16_t *dmrsSymbPos     = &rel15_ul->ul_dmrs_symb_pos;
francescomani's avatar
francescomani committed
515 516
  uint16_t *ptrsSymbPos = &pusch_vars->ptrs_symbols;
  uint8_t *ptrsSymbIdx = &pusch_vars->ptrs_symbol_index;
517 518
  uint16_t *nb_rb           = &rel15_ul->rb_size;
  uint8_t  *ptrsReOffset    = &rel15_ul->pusch_ptrs.ptrs_ports_list[0].ptrs_re_offset;
519

520
  /* loop over antennas */
521
  for (int aarx=0; aarx<frame_parms->nb_antennas_rx; aarx++) {
francescomani's avatar
francescomani committed
522 523
    c16_t *phase_per_symbol = (c16_t *)pusch_vars->ptrs_phase_per_slot[aarx];
    ptrs_re_symbol = &pusch_vars->ptrs_re_per_slot;
524
    *ptrs_re_symbol = 0;
525
    phase_per_symbol[symbol].i = 0; 
526
    /* set DMRS estimates to 0 angle with magnitude 1 */
527
    if(is_dmrs_symbol(symbol,*dmrsSymbPos)) {
528
      /* set DMRS real estimation to 32767 */
529
      phase_per_symbol[symbol].r=INT16_MAX; // 32767
530
#ifdef DEBUG_UL_PTRS
531
      printf("[PHY][PTRS]: DMRS Symbol %d -> %4d + j*%4d\n", symbol, phase_per_symbol[symbol].r,phase_per_symbol[symbol].i);
532 533
#endif
    }
534
    else {// real ptrs value is set to 0
535
      phase_per_symbol[symbol].r = 0; 
536 537
    }

538
    if(symbol == *startSymbIndex) {
539 540 541 542 543 544
      *ptrsSymbPos = 0;
      set_ptrs_symb_idx(ptrsSymbPos,
                        *nbSymb,
                        *startSymbIndex,
                        1<< *L_ptrs,
                        *dmrsSymbPos);
545
    }
546

547
    /* if not PTRS symbol set current ptrs symbol index to zero*/
548
    *ptrsSymbIdx = 0;
549

550
    /* Check if current symbol contains PTRS */
551
    if(is_ptrs_symbol(symbol, *ptrsSymbPos)) {
552
      *ptrsSymbIdx = symbol;
553
      /*------------------------------------------------------------------------------------------------------- */
554
      /* 1) Estimate common phase error per PTRS symbol                                                                */
555
      /*------------------------------------------------------------------------------------------------------- */
francescomani's avatar
francescomani committed
556 557 558
      nr_ptrs_cpe_estimation(*K_ptrs,
                             *ptrsReOffset,
                             *nb_rb,
559 560
                             rel15_ul->rnti,
                             nr_tti_rx,
francescomani's avatar
francescomani committed
561 562 563
                             symbol,
                             frame_parms->ofdm_symbol_size,
                             (int16_t *)&pusch_vars->rxdataF_comp[aarx][(symbol * nb_re_pusch)],
564
                             gNB->nr_gold_pusch_dmrs[rel15_ul->scid][nr_tti_rx][symbol],
francescomani's avatar
francescomani committed
565
                             (int16_t *)&phase_per_symbol[symbol],
566
                             ptrs_re_symbol);
567
    }
568

569
    /* For last OFDM symbol at each antenna perform interpolation and compensation for the slot*/
570
    if(symbol == (symbInSlot -1)) {
571 572 573 574
      /*------------------------------------------------------------------------------------------------------- */
      /* 2) Interpolate PTRS estimated value in TD */
      /*------------------------------------------------------------------------------------------------------- */
      /* If L-PTRS is > 0 then we need interpolation */
575
      if(*L_ptrs > 0) {
576
        ret = nr_ptrs_process_slot(*dmrsSymbPos, *ptrsSymbPos, (int16_t*)phase_per_symbol, *startSymbIndex, *nbSymb);
577
        if(ret != 0) {
578
          LOG_W(PHY,"[PTRS] Compensation is skipped due to error in PTRS slot processing !!\n");
579
        }
580
      }
581

582 583 584
      /*------------------------------------------------------------------------------------------------------- */
      /* 3) Compensated DMRS based estimated signal with PTRS estimation                                        */
      /*--------------------------------------------------------------------------------------------------------*/
585
      for(uint8_t i = *startSymbIndex; i < symbInSlot; i++) {
586 587
        /* DMRS Symbol has 0 phase so no need to rotate the respective symbol */
        /* Skip rotation if the slot processing is wrong */
588
        if((!is_dmrs_symbol(i,*dmrsSymbPos)) && (ret == 0)) {
589
#ifdef DEBUG_UL_PTRS
590
          printf("[PHY][UL][PTRS]: Rotate Symbol %2d with  %d + j* %d\n", i, phase_per_symbol[i].r,phase_per_symbol[i].i);
591
#endif
Tsung-Yu Chan's avatar
Tsung-Yu Chan committed
592
          rotate_cpx_vector((c16_t *)&pusch_vars->rxdataF_comp[aarx][i * nb_re_pusch],
593
                            &phase_per_symbol[i],
Tsung-Yu Chan's avatar
Tsung-Yu Chan committed
594
                            (c16_t *)&pusch_vars->rxdataF_comp[aarx][i * nb_re_pusch],
francescomani's avatar
francescomani committed
595 596
                            ((*nb_rb) * NR_NB_SC_PER_RB),
                            15);
Tsung-Yu Chan's avatar
Tsung-Yu Chan committed
597 598 599 600
        } // if not DMRS Symbol
      } // symbol loop
    } // last symbol check
  } // Antenna loop
601
}
602

603
uint32_t calc_power(const int16_t *x, const uint32_t size) {
604 605
  int64_t sum_x = 0;
  int64_t sum_x2 = 0;
606 607 608 609
  for(int k = 0; k<size; k++) {
    sum_x = sum_x + x[k];
    sum_x2 = sum_x2 + x[k]*x[k];
  }
610

611 612 613
  return sum_x2/size - (sum_x/size)*(sum_x/size);
}

614 615 616 617 618 619 620 621 622 623 624 625 626 627 628
int nr_srs_channel_estimation(
    const PHY_VARS_gNB *gNB,
    const int frame,
    const int slot,
    const nfapi_nr_srs_pdu_t *srs_pdu,
    const nr_srs_info_t *nr_srs_info,
    const c16_t **srs_generated_signal,
    int32_t srs_received_signal[][gNB->frame_parms.ofdm_symbol_size * (1 << srs_pdu->num_symbols)],
    int32_t srs_estimated_channel_freq[][1 << srs_pdu->num_ant_ports]
                                      [gNB->frame_parms.ofdm_symbol_size * (1 << srs_pdu->num_symbols)],
    int32_t srs_estimated_channel_time[][1 << srs_pdu->num_ant_ports][gNB->frame_parms.ofdm_symbol_size],
    int32_t srs_estimated_channel_time_shifted[][1 << srs_pdu->num_ant_ports][gNB->frame_parms.ofdm_symbol_size],
    int8_t *snr_per_rb,
    int8_t *snr)
{
629 630 631 632
#ifdef SRS_DEBUG
  LOG_I(NR_PHY,"Calling %s function\n", __FUNCTION__);
#endif

633
  const NR_DL_FRAME_PARMS *frame_parms = &gNB->frame_parms;
634
  const uint64_t subcarrier_offset = frame_parms->first_carrier_offset + srs_pdu->bwp_start*NR_NB_SC_PER_RB;
635

636 637 638 639
  const uint8_t N_ap = 1<<srs_pdu->num_ant_ports;
  const uint8_t K_TC = 2<<srs_pdu->comb_size;
  const uint16_t m_SRS_b = srs_bandwidth_config[srs_pdu->config_index][srs_pdu->bandwidth_index][0];
  const uint16_t M_sc_b_SRS = m_SRS_b * NR_NB_SC_PER_RB/K_TC;
640 641 642 643 644
  uint8_t fd_cdm = N_ap;
  if (N_ap == 4 && ((K_TC == 2 && srs_pdu->cyclic_shift >= 4) || (K_TC == 4 && srs_pdu->cyclic_shift >= 6))) {
    fd_cdm = 2;
  }

645
  c16_t srs_ls_estimated_channel[frame_parms->ofdm_symbol_size*(1<<srs_pdu->num_symbols)];
646
  uint32_t noise_power_per_rb[srs_pdu->bwp_size];
647 648 649 650 651 652 653 654 655 656 657 658 659 660 661

  const uint32_t arr_len = frame_parms->nb_antennas_rx * N_ap * M_sc_b_SRS;

  int16_t ch_real[arr_len];
  memset(ch_real, 0, arr_len * sizeof(int16_t));
  
  int16_t ch_imag[arr_len];
  memset(ch_imag, 0, arr_len * sizeof(int16_t));
 
  int16_t noise_real[arr_len];
  memset(noise_real, 0, arr_len * sizeof(int16_t));
 
  int16_t noise_imag[arr_len];
  memset(noise_imag, 0, arr_len * sizeof(int16_t));

662 663
  int16_t ls_estimated[2];

664
  uint8_t mem_offset = ((16 - ((long)&srs_estimated_channel_freq[0][0][subcarrier_offset + nr_srs_info->k_0_p[0][0]])) & 0xF) >> 2; // >> 2 <=> /sizeof(int32_t)
665

666 667 668 669 670
  // filt16_end is {4096,8192,8192,8192,12288,16384,16384,16384,0,0,0,0,0,0,0,0}
  // The End of OFDM symbol corresponds to the position of last 16384 in the filter
  // The multadd_real_vector_complex_scalar applies the remaining 8 zeros of filter, therefore, to avoid a buffer overflow,
  // we added 8 in the array size
  int32_t srs_est[frame_parms->ofdm_symbol_size*(1<<srs_pdu->num_symbols) + mem_offset + 8] __attribute__ ((aligned(32)));
671

672 673
  for (int ant = 0; ant < frame_parms->nb_antennas_rx; ant++) {

674 675
    for (int p_index = 0; p_index < N_ap; p_index++) {

676
      memset(srs_ls_estimated_channel, 0, frame_parms->ofdm_symbol_size*(1<<srs_pdu->num_symbols)*sizeof(c16_t));
677
      memset(srs_est, 0, (frame_parms->ofdm_symbol_size*(1<<srs_pdu->num_symbols) + mem_offset)*sizeof(int32_t));
678 679 680 681 682 683 684 685 686 687

#ifdef SRS_DEBUG
      LOG_I(NR_PHY,"====================== UE port %d --> gNB Rx antenna %i ======================\n", p_index, ant);
#endif

      uint16_t subcarrier = subcarrier_offset + nr_srs_info->k_0_p[p_index][0];
      if (subcarrier>frame_parms->ofdm_symbol_size) {
        subcarrier -= frame_parms->ofdm_symbol_size;
      }

688
      int16_t *srs_estimated_channel16 = (int16_t *)&srs_est[subcarrier + mem_offset];
689 690 691 692 693 694 695 696 697 698

      for (int k = 0; k < M_sc_b_SRS; k++) {

        if (k%fd_cdm==0) {

          ls_estimated[0] = 0;
          ls_estimated[1] = 0;
          uint16_t subcarrier_cdm = subcarrier;

          for (int cdm_idx = 0; cdm_idx < fd_cdm; cdm_idx++) {
699 700
            int16_t generated_real = srs_generated_signal[p_index][subcarrier_cdm].r;
            int16_t generated_imag = srs_generated_signal[p_index][subcarrier_cdm].i;
701

702 703
            int16_t received_real = ((c16_t*)srs_received_signal[ant])[subcarrier_cdm].r;
            int16_t received_imag = ((c16_t*)srs_received_signal[ant])[subcarrier_cdm].i;
704 705 706 707 708 709 710 711 712 713 714 715 716 717

            // We know that nr_srs_info->srs_generated_signal_bits bits are enough to represent the generated_real and generated_imag.
            // So we only need a nr_srs_info->srs_generated_signal_bits shift to ensure that the result fits into 16 bits.
            ls_estimated[0] += (int16_t)(((int32_t)generated_real*received_real + (int32_t)generated_imag*received_imag)>>nr_srs_info->srs_generated_signal_bits);
            ls_estimated[1] += (int16_t)(((int32_t)generated_real*received_imag - (int32_t)generated_imag*received_real)>>nr_srs_info->srs_generated_signal_bits);

            // Subcarrier increment
            subcarrier_cdm += K_TC;
            if (subcarrier_cdm >= frame_parms->ofdm_symbol_size) {
              subcarrier_cdm=subcarrier_cdm-frame_parms->ofdm_symbol_size;
            }
          }
        }

718 719
        srs_ls_estimated_channel[subcarrier].r = ls_estimated[0];
        srs_ls_estimated_channel[subcarrier].i = ls_estimated[1];
720 721 722 723 724 725 726 727 728 729 730 731

#ifdef SRS_DEBUG
        int subcarrier_log = subcarrier-subcarrier_offset;
        if(subcarrier_log < 0) {
          subcarrier_log = subcarrier_log + frame_parms->ofdm_symbol_size;
        }
        if(subcarrier_log%12 == 0) {
          LOG_I(NR_PHY,"------------------------------------ %d ------------------------------------\n", subcarrier_log/12);
          LOG_I(NR_PHY,"\t  __genRe________genIm__|____rxRe_________rxIm__|____lsRe________lsIm_\n");
        }
        LOG_I(NR_PHY,"(%4i) %6i\t%6i  |  %6i\t%6i  |  %6i\t%6i\n",
              subcarrier_log,
732 733
              ((c16_t*)srs_generated_signal[p_index])[subcarrier].r, ((c16_t*)srs_generated_signal[p_index])[subcarrier].i,
              ((c16_t*)srs_received_signal[ant])[subcarrier].r, ((c16_t*)srs_received_signal[ant])[subcarrier].i,
734 735 736
              ls_estimated[0], ls_estimated[1]);
#endif

737 738
        const uint16_t sc_offset = subcarrier + mem_offset;

739 740 741 742 743 744 745
        // Channel interpolation
        if(srs_pdu->comb_size == 0) {
          if(k == 0) { // First subcarrier case
            // filt8_start is {12288,8192,4096,0,0,0,0,0}
            multadd_real_vector_complex_scalar(filt8_start, ls_estimated, srs_estimated_channel16, 8);
          } else if(subcarrier < K_TC) { // Start of OFDM symbol case
            // filt8_start is {12288,8192,4096,0,0,0,0,0}
746
            srs_estimated_channel16 = (int16_t *)&srs_est[subcarrier];
laurent's avatar
laurent committed
747
            const short *filter = mem_offset == 0 ? filt8_start : filt8_start_shift2;
748
            multadd_real_vector_complex_scalar(filter, ls_estimated, srs_estimated_channel16, 8);
749 750
          } else if((subcarrier+K_TC)>=frame_parms->ofdm_symbol_size || k == (M_sc_b_SRS-1)) { // End of OFDM symbol or last subcarrier cases
            // filt8_end is {4096,8192,12288,16384,0,0,0,0}
laurent's avatar
laurent committed
751
            const short *filter = mem_offset == 0 || k == (M_sc_b_SRS - 1) ? filt8_end : filt8_end_shift2;
752
            multadd_real_vector_complex_scalar(filter, ls_estimated, srs_estimated_channel16, 8);
753 754 755 756 757 758
          } else if(k%2 == 1) { // 1st middle case
            // filt8_middle2 is {4096,8192,8192,8192,4096,0,0,0}
            multadd_real_vector_complex_scalar(filt8_middle2, ls_estimated, srs_estimated_channel16, 8);
          } else if(k%2 == 0) { // 2nd middle case
            // filt8_middle4 is {0,0,4096,8192,8192,8192,4096,0}
            multadd_real_vector_complex_scalar(filt8_middle4, ls_estimated, srs_estimated_channel16, 8);
759
            srs_estimated_channel16 = (int16_t *)&srs_est[sc_offset];
760 761 762 763 764 765
          }
        } else {
          if(k == 0) { // First subcarrier case
            // filt16_start is {12288,8192,8192,8192,4096,0,0,0,0,0,0,0,0,0,0,0}
            multadd_real_vector_complex_scalar(filt16_start, ls_estimated, srs_estimated_channel16, 16);
          } else if(subcarrier < K_TC) { // Start of OFDM symbol case
766
            srs_estimated_channel16 = (int16_t *)&srs_est[sc_offset];
767 768 769 770 771 772 773 774
            // filt16_start is {12288,8192,8192,8192,4096,0,0,0,0,0,0,0,0,0,0,0}
            multadd_real_vector_complex_scalar(filt16_start, ls_estimated, srs_estimated_channel16, 16);
          } else if((subcarrier+K_TC)>=frame_parms->ofdm_symbol_size || k == (M_sc_b_SRS-1)) { // End of OFDM symbol or last subcarrier cases
            // filt16_end is {4096,8192,8192,8192,12288,16384,16384,16384,0,0,0,0,0,0,0,0}
            multadd_real_vector_complex_scalar(filt16_end, ls_estimated, srs_estimated_channel16, 16);
          } else { // Middle case
            // filt16_middle4 is {4096,8192,8192,8192,8192,8192,8192,8192,4096,0,0,0,0,0,0,0}
            multadd_real_vector_complex_scalar(filt16_middle4, ls_estimated, srs_estimated_channel16, 16);
775
            srs_estimated_channel16 = (int16_t *)&srs_est[sc_offset];
776 777 778 779 780 781 782
          }
        }

        // Subcarrier increment
        subcarrier += K_TC;
        if (subcarrier >= frame_parms->ofdm_symbol_size) {
          subcarrier=subcarrier-frame_parms->ofdm_symbol_size;
783
        }
784 785 786

      } // for (int k = 0; k < M_sc_b_SRS; k++)

787 788 789 790
      memcpy(&srs_estimated_channel_freq[ant][p_index][0],
             &srs_est[mem_offset],
             (frame_parms->ofdm_symbol_size*(1<<srs_pdu->num_symbols))*sizeof(int32_t));

791 792 793 794 795 796 797
      // Compute noise
      subcarrier = subcarrier_offset + nr_srs_info->k_0_p[p_index][0];
      if (subcarrier>frame_parms->ofdm_symbol_size) {
        subcarrier -= frame_parms->ofdm_symbol_size;
      }
      uint16_t base_idx = ant*N_ap*M_sc_b_SRS + p_index*M_sc_b_SRS;
      for (int k = 0; k < M_sc_b_SRS; k++) {
798 799
        ch_real[base_idx+k] = ((c16_t*)srs_estimated_channel_freq[ant][p_index])[subcarrier].r;
        ch_imag[base_idx+k] = ((c16_t*)srs_estimated_channel_freq[ant][p_index])[subcarrier].i;
800 801
        noise_real[base_idx+k] = abs(srs_ls_estimated_channel[subcarrier].r - ch_real[base_idx+k]);
        noise_imag[base_idx+k] = abs(srs_ls_estimated_channel[subcarrier].i - ch_imag[base_idx+k]);
802 803 804
        subcarrier += K_TC;
        if (subcarrier >= frame_parms->ofdm_symbol_size) {
          subcarrier=subcarrier-frame_parms->ofdm_symbol_size;
805
        }
806
      }
807 808

#ifdef SRS_DEBUG
809 810 811
      subcarrier = subcarrier_offset + nr_srs_info->k_0_p[p_index][0];
      if (subcarrier>frame_parms->ofdm_symbol_size) {
        subcarrier -= frame_parms->ofdm_symbol_size;
812
      }
813 814 815 816 817 818 819 820 821 822 823 824 825 826 827

      for (int k = 0; k < K_TC*M_sc_b_SRS; k++) {

        int subcarrier_log = subcarrier-subcarrier_offset;
        if(subcarrier_log < 0) {
          subcarrier_log = subcarrier_log + frame_parms->ofdm_symbol_size;
        }

        if(subcarrier_log%12 == 0) {
          LOG_I(NR_PHY,"------------------------------------- %d -------------------------------------\n", subcarrier_log/12);
          LOG_I(NR_PHY,"\t  __lsRe__________lsIm__|____intRe_______intIm__|____noiRe_______noiIm__\n");
        }

        LOG_I(NR_PHY,"(%4i) %6i\t%6i  |  %6i\t%6i  |  %6i\t%6i\n",
              subcarrier_log,
828 829 830 831
              srs_ls_estimated_channel[subcarrier].r,
              srs_ls_estimated_channel[subcarrier].i,
              ((c16_t*)srs_estimated_channel_freq[ant][p_index])[subcarrier].r,
              ((c16_t*)srs_estimated_channel_freq[ant][p_index])[subcarrier].i,
832 833 834 835 836 837 838
              noise_real[base_idx+(k/K_TC)], noise_imag[base_idx+(k/K_TC)]);

        // Subcarrier increment
        subcarrier++;
        if (subcarrier >= frame_parms->ofdm_symbol_size) {
          subcarrier=subcarrier-frame_parms->ofdm_symbol_size;
        }
839 840
      }
#endif
841

842 843 844 845
      // Convert to time domain
      freq2time(gNB->frame_parms.ofdm_symbol_size,
                (int16_t*) srs_estimated_channel_freq[ant][p_index],
                (int16_t*) srs_estimated_channel_time[ant][p_index]);
846

847 848 849
      memcpy(&srs_estimated_channel_time_shifted[ant][p_index][0],
             &srs_estimated_channel_time[ant][p_index][gNB->frame_parms.ofdm_symbol_size>>1],
             (gNB->frame_parms.ofdm_symbol_size>>1)*sizeof(int32_t));
rmagueta's avatar
rmagueta committed
850

851 852 853
      memcpy(&srs_estimated_channel_time_shifted[ant][p_index][gNB->frame_parms.ofdm_symbol_size>>1],
             &srs_estimated_channel_time[ant][p_index][0],
             (gNB->frame_parms.ofdm_symbol_size>>1)*sizeof(int32_t));
rmagueta's avatar
rmagueta committed
854

855 856
    } // for (int p_index = 0; p_index < N_ap; p_index++)
  } // for (int ant = 0; ant < frame_parms->nb_antennas_rx; ant++)
857

858
  // Compute signal power
859
  uint32_t signal_power = calc_power(ch_real, arr_len) + calc_power(ch_imag, arr_len);
860 861

#ifdef SRS_DEBUG
862
  LOG_I(NR_PHY,"signal_power = %u\n", signal_power);
863 864
#endif

865
  if (signal_power == 0) {
866 867 868 869
    LOG_W(NR_PHY, "Received SRS signal power is 0\n");
    return -1;
  }

870 871
  // Compute noise power

872
  const uint8_t signal_power_bits = log2_approx(signal_power);
873 874
  const uint8_t factor_bits = signal_power_bits < 32 ? 32 - signal_power_bits : 0; // 32 due to input of dB_fixed(uint32_t x)
  const int32_t factor_dB = dB_fixed(1<<factor_bits);
875

876
  const uint8_t srs_symbols_per_rb = srs_pdu->comb_size == 0 ? 6 : 3;
877
  const uint8_t n_noise_est = frame_parms->nb_antennas_rx*N_ap*srs_symbols_per_rb;
878 879 880 881 882
  uint64_t sum_re = 0;
  uint64_t sum_re2 = 0;
  uint64_t sum_im = 0;
  uint64_t sum_im2 = 0;

883
  for (int rb = 0; rb < m_SRS_b; rb++) {
884

885 886 887 888
    sum_re = 0;
    sum_re2 = 0;
    sum_im = 0;
    sum_im2 = 0;
889 890

    for (int ant = 0; ant < frame_parms->nb_antennas_rx; ant++) {
891 892 893 894 895 896 897 898 899 900 901
      for (int p_index = 0; p_index < N_ap; p_index++) {
        uint16_t base_idx = ant*N_ap*M_sc_b_SRS + p_index*M_sc_b_SRS + rb*srs_symbols_per_rb;
        for (int srs_symb = 0; srs_symb < srs_symbols_per_rb; srs_symb++) {
          sum_re = sum_re + noise_real[base_idx+srs_symb];
          sum_re2 = sum_re2 + noise_real[base_idx+srs_symb]*noise_real[base_idx+srs_symb];
          sum_im = sum_im + noise_imag[base_idx+srs_symb];
          sum_im2 = sum_im2 + noise_imag[base_idx+srs_symb]*noise_imag[base_idx+srs_symb];
        } // for (int srs_symb = 0; srs_symb < srs_symbols_per_rb; srs_symb++)
      } // for (int p_index = 0; p_index < N_ap; p_index++)
    } // for (int ant = 0; ant < frame_parms->nb_antennas_rx; ant++)

902 903
    noise_power_per_rb[rb] = max(sum_re2 / n_noise_est - (sum_re / n_noise_est) * (sum_re / n_noise_est) +
                                 sum_im2 / n_noise_est - (sum_im / n_noise_est) * (sum_im / n_noise_est), 1);
904
    snr_per_rb[rb] = dB_fixed((int32_t)((signal_power<<factor_bits)/noise_power_per_rb[rb])) - factor_dB;
905 906

#ifdef SRS_DEBUG
907
    LOG_I(NR_PHY,"noise_power_per_rb[%i] = %i, snr_per_rb[%i] = %i dB\n", rb, noise_power_per_rb[rb], rb, snr_per_rb[rb]);
908 909
#endif

910
  } // for (int rb = 0; rb < m_SRS_b; rb++)
911

912
  const uint32_t noise_power = max(calc_power(noise_real, arr_len) + calc_power(noise_imag, arr_len), 1);
913

914
  *snr = dB_fixed((int32_t)((signal_power<<factor_bits)/(noise_power))) - factor_dB;
915

916
#ifdef SRS_DEBUG
917
  LOG_I(NR_PHY,"noise_power = %u, SNR = %i dB\n", noise_power, *snr);
918 919
#endif

920
  return 0;
921
}