lte_ul_channel_estimation.c 29.7 KB
Newer Older
1
2
3
4
5
/*
 * 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
Cedric Roux's avatar
Cedric Roux committed
6
 * the OAI Public License, Version 1.1  (the "License"); you may not use this file
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 * 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
 */
ghaddab's avatar
ghaddab committed
21

22
23
#include "PHY/defs_eNB.h"
#include "PHY/phy_extern.h"
24
#include "PHY/sse_intrin.h"
25
//#define DEBUG_CH
26
#include "common/utils/LOG/log.h"
27
#include "PHY/LTE_TRANSPORT/transport_common_proto.h"
knopp's avatar
knopp committed
28
#include "lte_estimation.h"
29
30

// round(exp(sqrt(-1)*(pi/2)*[0:1:N-1]/N)*pow2(15))
gauthier's avatar
gauthier committed
31
static int16_t ru_90[2*128] = {32767, 0,32766, 402,32758, 804,32746, 1206,32729, 1608,32706, 2009,32679, 2411,32647, 2811,32610, 3212,32568, 3612,32522, 4011,32470, 4410,32413, 4808,32352, 5205,32286, 5602,32214, 5998,32138, 6393,32058, 6787,31972, 7180,31881, 7571,31786, 7962,31686, 8351,31581, 8740,31471, 9127,31357, 9512,31238, 9896,31114, 10279,30986, 10660,30853, 11039,30715, 11417,30572, 11793,30425, 12167,30274, 12540,30118, 12910,29957, 13279,29792, 13646,29622, 14010,29448, 14373,29269, 14733,29086, 15091,28899, 15447,28707, 15800,28511, 16151,28311, 16500,28106, 16846,27897, 17190,27684, 17531,27467, 17869,27246, 18205,27020, 18538,26791, 18868,26557, 19195,26320, 19520,26078, 19841,25833, 20160,25583, 20475,25330, 20788,25073, 21097,24812, 21403,24548, 21706,24279, 22006,24008, 22302,23732, 22595,23453, 22884,23170, 23170,22884, 23453,22595, 23732,22302, 24008,22006, 24279,21706, 24548,21403, 24812,21097, 25073,20788, 25330,20475, 25583,20160, 25833,19841, 26078,19520, 26320,19195, 26557,18868, 26791,18538, 27020,18205, 27246,17869, 27467,17531, 27684,17190, 27897,16846, 28106,16500, 28311,16151, 28511,15800, 28707,15447, 28899,15091, 29086,14733, 29269,14373, 29448,14010, 29622,13646, 29792,13279, 29957,12910, 30118,12540, 30274,12167, 30425,11793, 30572,11417, 30715,11039, 30853,10660, 30986,10279, 31114,9896, 31238,9512, 31357,9127, 31471,8740, 31581,8351, 31686,7962, 31786,7571, 31881,7180, 31972,6787, 32058,6393, 32138,5998, 32214,5602, 32286,5205, 32352,4808, 32413,4410, 32470,4011, 32522,3612, 32568,3212, 32610,2811, 32647,2411, 32679,2009, 32706,1608, 32729,1206, 32746,804, 32758,402, 32766};
32

gauthier's avatar
gauthier committed
33
static int16_t ru_90c[2*128] = {32767, 0,32766, -402,32758, -804,32746, -1206,32729, -1608,32706, -2009,32679, -2411,32647, -2811,32610, -3212,32568, -3612,32522, -4011,32470, -4410,32413, -4808,32352, -5205,32286, -5602,32214, -5998,32138, -6393,32058, -6787,31972, -7180,31881, -7571,31786, -7962,31686, -8351,31581, -8740,31471, -9127,31357, -9512,31238, -9896,31114, -10279,30986, -10660,30853, -11039,30715, -11417,30572, -11793,30425, -12167,30274, -12540,30118, -12910,29957, -13279,29792, -13646,29622, -14010,29448, -14373,29269, -14733,29086, -15091,28899, -15447,28707, -15800,28511, -16151,28311, -16500,28106, -16846,27897, -17190,27684, -17531,27467, -17869,27246, -18205,27020, -18538,26791, -18868,26557, -19195,26320, -19520,26078, -19841,25833, -20160,25583, -20475,25330, -20788,25073, -21097,24812, -21403,24548, -21706,24279, -22006,24008, -22302,23732, -22595,23453, -22884,23170, -23170,22884, -23453,22595, -23732,22302, -24008,22006, -24279,21706, -24548,21403, -24812,21097, -25073,20788, -25330,20475, -25583,20160, -25833,19841, -26078,19520, -26320,19195, -26557,18868, -26791,18538, -27020,18205, -27246,17869, -27467,17531, -27684,17190, -27897,16846, -28106,16500, -28311,16151, -28511,15800, -28707,15447, -28899,15091, -29086,14733, -29269,14373, -29448,14010, -29622,13646, -29792,13279, -29957,12910, -30118,12540, -30274,12167, -30425,11793, -30572,11417, -30715,11039, -30853,10660, -30986,10279, -31114,9896, -31238,9512, -31357,9127, -31471,8740, -31581,8351, -31686,7962, -31786,7571, -31881,7180, -31972,6787, -32058,6393, -32138,5998, -32214,5602, -32286,5205, -32352,4808, -32413,4410, -32470,4011, -32522,3612, -32568,3212, -32610,2811, -32647,2411, -32679,2009, -32706,1608, -32729,1206, -32746,804, -32758,402, -32766};
34
35
36

#define SCALE 0x3FFF

37
int32_t lte_ul_channel_estimation(PHY_VARS_eNB *eNB,
38
				  L1_rxtx_proc_t *proc,
39
                                  module_id_t UE_id,
40
                                  unsigned char l,
41
                                  unsigned char Ns) {
42

43
44
  LTE_DL_FRAME_PARMS *frame_parms = &eNB->frame_parms;
  LTE_eNB_PUSCH *pusch_vars = eNB->pusch_vars[UE_id];
45
46
47
  int32_t **ul_ch_estimates=pusch_vars->drs_ch_estimates;
  int32_t **ul_ch_estimates_time=  pusch_vars->drs_ch_estimates_time;
  int32_t **rxdataF_ext=  pusch_vars->rxdataF_ext;
knopp's avatar
knopp committed
48
  int subframe = proc->subframe_rx;
49
  uint8_t harq_pid; 
gauthier's avatar
gauthier committed
50
51
52
53
  int16_t delta_phase = 0;
  int16_t *ru1 = ru_90;
  int16_t *ru2 = ru_90;
  int16_t current_phase1,current_phase2;
54

gauthier's avatar
gauthier committed
55
56
  uint16_t aa,Msc_RS,Msc_RS_idx;
  uint16_t * Msc_idx_ptr;
57
  int k,pilot_pos1 = 3 - frame_parms->Ncp, pilot_pos2 = 10 - 2*frame_parms->Ncp;
gauthier's avatar
gauthier committed
58
59
  int32_t *ul_ch1=NULL, *ul_ch2=NULL;
  int16_t ul_ch_estimates_re,ul_ch_estimates_im;
60

Xiwen JIANG's avatar
Xiwen JIANG committed
61
  //uint8_t nb_antennas_rx = frame_parms->nb_antenna_ports_eNB;
gauthier's avatar
gauthier committed
62
  uint8_t nb_antennas_rx = frame_parms->nb_antennas_rx;
63
  uint8_t cyclic_shift;
64

gauthier's avatar
gauthier committed
65
66
67
  uint32_t alpha_ind;
  uint32_t u=frame_parms->pusch_config_common.ul_ReferenceSignalsPUSCH.grouphop[Ns+(subframe<<1)];
  uint32_t v=frame_parms->pusch_config_common.ul_ReferenceSignalsPUSCH.seqhop[Ns+(subframe<<1)];
68

69

70
  int symbol_offset,i;
71
72
73
74

  //debug_msg("lte_ul_channel_estimation: cyclic shift %d\n",cyclicShift);


gauthier's avatar
gauthier committed
75
76
  int16_t alpha_re[12] = {32767, 28377, 16383,     0,-16384,  -28378,-32768,-28378,-16384,    -1, 16383, 28377};
  int16_t alpha_im[12] = {0,     16383, 28377, 32767, 28377,   16383,     0,-16384,-28378,-32768,-28378,-16384};
77

78
#if defined(__x86_64__) || defined(__i386__)
79
80
  __m128i *rxdataF128,*ul_ref128,*ul_ch128;
  __m128i mmtmpU0,mmtmpU1,mmtmpU2,mmtmpU3;
81
82
83
84
#elif defined(__arm__)
  int16x8_t *rxdataF128,*ul_ref128,*ul_ch128;
  int32x4_t mmtmp0,mmtmp1,mmtmp_re,mmtmp_im;
#endif
85

86
87
88
89
90
91
92
93
94
95
96
97
98
  int32_t temp_in_ifft_0[2048*2] __attribute__((aligned(32)));


#ifdef Rel14
  if (eNB->ulsch[UE_id]->ue_type > 0) harq_pid = 0;
  else
#endif
    {
      harq_pid = subframe2harq_pid(frame_parms,proc->frame_rx,subframe);
    }

  uint16_t N_rb_alloc = eNB->ulsch[UE_id]->harq_processes[harq_pid]->nb_rb;
  int32_t tmp_estimates[N_rb_alloc*12] __attribute__((aligned(16)));
99

100
101
102
  Msc_RS = N_rb_alloc*12;

  cyclic_shift = (frame_parms->pusch_config_common.ul_ReferenceSignalsPUSCH.cyclicShift +
103
                  eNB->ulsch[UE_id]->harq_processes[harq_pid]->n_DMRS2 +
104
                  frame_parms->pusch_config_common.ul_ReferenceSignalsPUSCH.nPRS[(subframe<<1)+Ns]) % 12;
105

106
  Msc_idx_ptr = (uint16_t*) bsearch(&Msc_RS, dftsizes, 34, sizeof(uint16_t), compareints);
107

108
109
110
  if (Msc_idx_ptr)
    Msc_RS_idx = Msc_idx_ptr - dftsizes;
  else {
111
    LOG_E(PHY,"lte_ul_channel_estimation: index for Msc_RS=%d not found\n",Msc_RS);
112
113
114
    return(-1);
  }

115
  LOG_D(PHY,"subframe %d, Ns %d, l %d, Msc_RS = %d, Msc_RS_idx = %d, u %d, v %d, cyclic_shift %d\n",subframe,Ns,l,Msc_RS, Msc_RS_idx,u,v,cyclic_shift);
knopp's avatar
   
knopp committed
116
#ifdef DEBUG_CH
117
118

  if (Ns==0)
Mongazon's avatar
Mongazon committed
119
    LOG_M("drs_seq0.m","drsseq0",ul_ref_sigs_rx[u][v][Msc_RS_idx],2*Msc_RS,2,1);
120
  else
Mongazon's avatar
Mongazon committed
121
    LOG_M("drs_seq1.m","drsseq1",ul_ref_sigs_rx[u][v][Msc_RS_idx],2*Msc_RS,2,1);
122

123
124
125
126
127
128
129
#endif


  if (l == (3 - frame_parms->Ncp)) {

    symbol_offset = frame_parms->N_RB_UL*12*(l+((7-frame_parms->Ncp)*(Ns&1)));

130
    for (aa=0; aa<nb_antennas_rx; aa++) {
131
132
      //           msg("Componentwise prod aa %d, symbol_offset %d,ul_ch_estimates %p,ul_ch_estimates[aa] %p,ul_ref_sigs_rx[0][0][Msc_RS_idx] %p\n",aa,symbol_offset,ul_ch_estimates,ul_ch_estimates[aa],ul_ref_sigs_rx[0][0][Msc_RS_idx]);

133
#if defined(__x86_64__) || defined(__i386__)
134
135
136
      rxdataF128 = (__m128i *)&rxdataF_ext[aa][symbol_offset];
      ul_ch128   = (__m128i *)&ul_ch_estimates[aa][symbol_offset];
      ul_ref128  = (__m128i *)ul_ref_sigs_rx[u][v][Msc_RS_idx];
137
138
139
140
141
#elif defined(__arm__)
      rxdataF128 = (int16x8_t *)&rxdataF_ext[aa][symbol_offset];
      ul_ch128   = (int16x8_t *)&ul_ch_estimates[aa][symbol_offset];
      ul_ref128  = (int16x8_t *)ul_ref_sigs_rx[u][v][Msc_RS_idx];
#endif
142

143
      for (i=0; i<Msc_RS/12; i++) {
144
#if defined(__x86_64__) || defined(__i386__)
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
        // multiply by conjugated channel
        mmtmpU0 = _mm_madd_epi16(ul_ref128[0],rxdataF128[0]);
        // mmtmpU0 contains real part of 4 consecutive outputs (32-bit)
        mmtmpU1 = _mm_shufflelo_epi16(ul_ref128[0],_MM_SHUFFLE(2,3,0,1));
        mmtmpU1 = _mm_shufflehi_epi16(mmtmpU1,_MM_SHUFFLE(2,3,0,1));
        mmtmpU1 = _mm_sign_epi16(mmtmpU1,*(__m128i*)&conjugate[0]);
        mmtmpU1 = _mm_madd_epi16(mmtmpU1,rxdataF128[0]);
        // mmtmpU1 contains imag part of 4 consecutive outputs (32-bit)
        mmtmpU0 = _mm_srai_epi32(mmtmpU0,15);
        mmtmpU1 = _mm_srai_epi32(mmtmpU1,15);
        mmtmpU2 = _mm_unpacklo_epi32(mmtmpU0,mmtmpU1);
        mmtmpU3 = _mm_unpackhi_epi32(mmtmpU0,mmtmpU1);

        ul_ch128[0] = _mm_packs_epi32(mmtmpU2,mmtmpU3);
        //  printf("rb %d ch: %d %d\n",i,((int16_t*)ul_ch128)[0],((int16_t*)ul_ch128)[1]);
        // multiply by conjugated channel
        mmtmpU0 = _mm_madd_epi16(ul_ref128[1],rxdataF128[1]);
        // mmtmpU0 contains real part of 4 consecutive outputs (32-bit)
        mmtmpU1 = _mm_shufflelo_epi16(ul_ref128[1],_MM_SHUFFLE(2,3,0,1));
        mmtmpU1 = _mm_shufflehi_epi16(mmtmpU1,_MM_SHUFFLE(2,3,0,1));
        mmtmpU1 = _mm_sign_epi16(mmtmpU1,*(__m128i*)conjugate);
        mmtmpU1 = _mm_madd_epi16(mmtmpU1,rxdataF128[1]);
        // mmtmpU1 contains imag part of 4 consecutive outputs (32-bit)
        mmtmpU0 = _mm_srai_epi32(mmtmpU0,15);
        mmtmpU1 = _mm_srai_epi32(mmtmpU1,15);
        mmtmpU2 = _mm_unpacklo_epi32(mmtmpU0,mmtmpU1);
        mmtmpU3 = _mm_unpackhi_epi32(mmtmpU0,mmtmpU1);

        ul_ch128[1] = _mm_packs_epi32(mmtmpU2,mmtmpU3);

        mmtmpU0 = _mm_madd_epi16(ul_ref128[2],rxdataF128[2]);
        // mmtmpU0 contains real part of 4 consecutive outputs (32-bit)
        mmtmpU1 = _mm_shufflelo_epi16(ul_ref128[2],_MM_SHUFFLE(2,3,0,1));
        mmtmpU1 = _mm_shufflehi_epi16(mmtmpU1,_MM_SHUFFLE(2,3,0,1));
        mmtmpU1 = _mm_sign_epi16(mmtmpU1,*(__m128i*)conjugate);
        mmtmpU1 = _mm_madd_epi16(mmtmpU1,rxdataF128[2]);
        // mmtmpU1 contains imag part of 4 consecutive outputs (32-bit)
        mmtmpU0 = _mm_srai_epi32(mmtmpU0,15);
        mmtmpU1 = _mm_srai_epi32(mmtmpU1,15);
        mmtmpU2 = _mm_unpacklo_epi32(mmtmpU0,mmtmpU1);
        mmtmpU3 = _mm_unpackhi_epi32(mmtmpU0,mmtmpU1);

        ul_ch128[2] = _mm_packs_epi32(mmtmpU2,mmtmpU3);
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
#elif defined(__arm__)
      mmtmp0 = vmull_s16(((int16x4_t*)ul_ref128)[0],((int16x4_t*)rxdataF128)[0]);
      mmtmp1 = vmull_s16(((int16x4_t*)ul_ref128)[1],((int16x4_t*)rxdataF128)[1]);
      mmtmp_re = vcombine_s32(vpadd_s32(vget_low_s32(mmtmp0),vget_high_s32(mmtmp0)),
                              vpadd_s32(vget_low_s32(mmtmp1),vget_high_s32(mmtmp1)));
      mmtmp0 = vmull_s16(vrev32_s16(vmul_s16(((int16x4_t*)ul_ref128)[0],*(int16x4_t*)conjugate)), ((int16x4_t*)rxdataF128)[0]);
      mmtmp1 = vmull_s16(vrev32_s16(vmul_s16(((int16x4_t*)ul_ref128)[1],*(int16x4_t*)conjugate)), ((int16x4_t*)rxdataF128)[1]);
      mmtmp_im = vcombine_s32(vpadd_s32(vget_low_s32(mmtmp0),vget_high_s32(mmtmp0)),
                              vpadd_s32(vget_low_s32(mmtmp1),vget_high_s32(mmtmp1)));

      ul_ch128[0] = vcombine_s16(vmovn_s32(mmtmp_re),vmovn_s32(mmtmp_im));
      ul_ch128++;
      ul_ref128++;
      rxdataF128++;
      mmtmp0 = vmull_s16(((int16x4_t*)ul_ref128)[0],((int16x4_t*)rxdataF128)[0]);
      mmtmp1 = vmull_s16(((int16x4_t*)ul_ref128)[1],((int16x4_t*)rxdataF128)[1]);
      mmtmp_re = vcombine_s32(vpadd_s32(vget_low_s32(mmtmp0),vget_high_s32(mmtmp0)),
                              vpadd_s32(vget_low_s32(mmtmp1),vget_high_s32(mmtmp1)));
      mmtmp0 = vmull_s16(vrev32_s16(vmul_s16(((int16x4_t*)ul_ref128)[0],*(int16x4_t*)conjugate)), ((int16x4_t*)rxdataF128)[0]);
      mmtmp1 = vmull_s16(vrev32_s16(vmul_s16(((int16x4_t*)ul_ref128)[1],*(int16x4_t*)conjugate)), ((int16x4_t*)rxdataF128)[1]);
      mmtmp_im = vcombine_s32(vpadd_s32(vget_low_s32(mmtmp0),vget_high_s32(mmtmp0)),
                              vpadd_s32(vget_low_s32(mmtmp1),vget_high_s32(mmtmp1)));

      ul_ch128[0] = vcombine_s16(vmovn_s32(mmtmp_re),vmovn_s32(mmtmp_im));
      ul_ch128++;
      ul_ref128++;
      rxdataF128++;

      mmtmp0 = vmull_s16(((int16x4_t*)ul_ref128)[0],((int16x4_t*)rxdataF128)[0]);
      mmtmp1 = vmull_s16(((int16x4_t*)ul_ref128)[1],((int16x4_t*)rxdataF128)[1]);
      mmtmp_re = vcombine_s32(vpadd_s32(vget_low_s32(mmtmp0),vget_high_s32(mmtmp0)),
                              vpadd_s32(vget_low_s32(mmtmp1),vget_high_s32(mmtmp1)));
      mmtmp0 = vmull_s16(vrev32_s16(vmul_s16(((int16x4_t*)ul_ref128)[0],*(int16x4_t*)conjugate)), ((int16x4_t*)rxdataF128)[0]);
      mmtmp1 = vmull_s16(vrev32_s16(vmul_s16(((int16x4_t*)ul_ref128)[1],*(int16x4_t*)conjugate)), ((int16x4_t*)rxdataF128)[1]);
      mmtmp_im = vcombine_s32(vpadd_s32(vget_low_s32(mmtmp0),vget_high_s32(mmtmp0)),
                              vpadd_s32(vget_low_s32(mmtmp1),vget_high_s32(mmtmp1)));

      ul_ch128[0] = vcombine_s16(vmovn_s32(mmtmp_re),vmovn_s32(mmtmp_im));
      ul_ch128++;
      ul_ref128++;
      rxdataF128++;

230

231
#endif
232
233
234
        ul_ch128+=3;
        ul_ref128+=3;
        rxdataF128+=3;
235
236
237
      }

      alpha_ind = 0;
238
239
240
241

      if((cyclic_shift != 0)) {
        // Compensating for the phase shift introduced at the transmitte
#ifdef DEBUG_CH
Mongazon's avatar
Mongazon committed
242
        LOG_M("drs_est_pre.m","drsest_pre",ul_ch_estimates[0],300*12,1,1);
243
#endif
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263

        for(i=symbol_offset; i<symbol_offset+Msc_RS; i++) {
          ul_ch_estimates_re = ((int16_t*) ul_ch_estimates[aa])[i<<1];
          ul_ch_estimates_im = ((int16_t*) ul_ch_estimates[aa])[(i<<1)+1];
          //    ((int16_t*) ul_ch_estimates[aa])[i<<1] =  (i%2 == 1? 1:-1) * ul_ch_estimates_re;
          ((int16_t*) ul_ch_estimates[aa])[i<<1] =
            (int16_t) (((int32_t) (alpha_re[alpha_ind]) * (int32_t) (ul_ch_estimates_re) +
                        (int32_t) (alpha_im[alpha_ind]) * (int32_t) (ul_ch_estimates_im))>>15);

          //((int16_t*) ul_ch_estimates[aa])[(i<<1)+1] =  (i%2 == 1? 1:-1) * ul_ch_estimates_im;
          ((int16_t*) ul_ch_estimates[aa])[(i<<1)+1] =
            (int16_t) (((int32_t) (alpha_re[alpha_ind]) * (int32_t) (ul_ch_estimates_im) -
                        (int32_t) (alpha_im[alpha_ind]) * (int32_t) (ul_ch_estimates_re))>>15);

          alpha_ind+=cyclic_shift;

          if (alpha_ind>11)
            alpha_ind-=12;
        }

264
#ifdef DEBUG_CH
Mongazon's avatar
Mongazon committed
265
        LOG_M("drs_est_post.m","drsest_post",ul_ch_estimates[0],300*12,1,1);
266
#endif
267
268
269
      }

      // Convert to time domain for visualization
Cedric Roux's avatar
Cedric Roux committed
270
      memset(temp_in_ifft_0,0,frame_parms->ofdm_symbol_size*sizeof(int32_t));
271
272
      for(i=0; i<Msc_RS; i++)
        ((int32_t*)temp_in_ifft_0)[i] = ul_ch_estimates[aa][symbol_offset+i];
Cedric Roux's avatar
Cedric Roux committed
273

274
275
      switch(frame_parms->N_RB_DL) {
      case 6:
276
	idft128((int16_t*) temp_in_ifft_0,
277
278
279
280
	       (int16_t*) ul_ch_estimates_time[aa],
	       1);
	break;
      case 25:
281
	idft512((int16_t*) temp_in_ifft_0,
282
283
284
285
	       (int16_t*) ul_ch_estimates_time[aa],
	       1);
	break;
      case 50:
286
	idft1024((int16_t*) temp_in_ifft_0,
287
288
289
290
	       (int16_t*) ul_ch_estimates_time[aa],
	       1);
	break;
      case 100:
291
	idft2048((int16_t*) temp_in_ifft_0,
292
293
294
	       (int16_t*) ul_ch_estimates_time[aa],
	       1);
	break;
295
296
      }

Cedric Roux's avatar
Cedric Roux committed
297
#if T_TRACER
298
      if (aa == 0)
Cedric Roux's avatar
Cedric Roux committed
299
        T(T_ENB_PHY_UL_CHANNEL_ESTIMATE, T_INT(0), T_INT(eNB->ulsch[UE_id]->rnti),
Cedric Roux's avatar
Cedric Roux committed
300
          T_INT(proc->frame_rx), T_INT(subframe),
301
          T_INT(0), T_BUFFER(ul_ch_estimates_time[0], 512  * 4));
Cedric Roux's avatar
Cedric Roux committed
302
#endif
303

304
305
#ifdef DEBUG_CH

knopp's avatar
knopp committed
306
      if (aa==1) {
307
        if (Ns == 0) {
Mongazon's avatar
Mongazon committed
308
309
310
          LOG_M("rxdataF_ext.m","rxF_ext",&rxdataF_ext[aa][symbol_offset],512*2,2,1);
          LOG_M("tmpin_ifft.m","drs_in",temp_in_ifft_0,512,1,1);
          LOG_M("drs_est0.m","drs0",ul_ch_estimates_time[aa],512,1,1);
311
        } else
Mongazon's avatar
Mongazon committed
312
          LOG_M("drs_est1.m","drs1",ul_ch_estimates_time[aa],512,1,1);
313
      }
314

315
316
317
#endif


318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
    
    
	if (Ns&1) {//we are in the second slot of the sub-frame, so do the interpolation

	  ul_ch1 = &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*pilot_pos1];
	  ul_ch2 = &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*pilot_pos2];
	  
	  
	  // Estimation of phase difference between the 2 channel estimates
	  delta_phase = lte_ul_freq_offset_estimation(frame_parms,
						      ul_ch_estimates[aa],
						      N_rb_alloc);
	  // negative phase index indicates negative Im of ru
	  //    msg("delta_phase: %d\n",delta_phase);
	  
333
#ifdef DEBUG_CH
334
	  LOG_D(PHY,"lte_ul_channel_estimation: ul_ch1 = %p, ul_ch2 = %p, pilot_pos1=%d, pilot_pos2=%d\n",ul_ch1, ul_ch2, pilot_pos1,pilot_pos2);
335
#endif
336
337
338
339
340
341
342
343
	  
	  for (k=0; k<frame_parms->symbols_per_tti; k++) {
	    
	    // we scale alpha and beta by SCALE (instead of 0x7FFF) to avoid overflows
	    //	    alpha = (int16_t) (((int32_t) SCALE * (int32_t) (pilot_pos2-k))/(pilot_pos2-pilot_pos1));
	    //	    beta  = (int16_t) (((int32_t) SCALE * (int32_t) (k-pilot_pos1))/(pilot_pos2-pilot_pos1));

	    
344
#ifdef DEBUG_CH
345
	    LOG_D(PHY,"lte_ul_channel_estimation: k=%d, alpha = %d, beta = %d\n",k,alpha,beta);
346
#endif
347
          //symbol_offset_subframe = frame_parms->N_RB_UL*12*k;
348
	    
349
          // interpolate between estimates
350
351
352
353
354
355
356
357
358
359
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
	    if ((k != pilot_pos1) && (k != pilot_pos2))  {
	      //          multadd_complex_vector_real_scalar((int16_t*) ul_ch1,alpha,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],1,Msc_RS);
	      //          multadd_complex_vector_real_scalar((int16_t*) ul_ch2,beta ,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],0,Msc_RS);
	      
	      //          multadd_complex_vector_real_scalar((int16_t*) ul_ch1,SCALE,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],1,Msc_RS);
	      //          multadd_complex_vector_real_scalar((int16_t*) ul_ch2,SCALE,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],0,Msc_RS);
	      //          msg("phase = %d\n",ru[2*cmax(((delta_phase/7)*(k-3)),0)]);
	      
	      // the phase is linearly interpolated
	      current_phase1 = (delta_phase/7)*(k-pilot_pos1);
	      current_phase2 = (delta_phase/7)*(k-pilot_pos2);
	      //          msg("sym: %d, current_phase1: %d, current_phase2: %d\n",k,current_phase1,current_phase2);
	      // set the right quadrant
	      (current_phase1 > 0) ? (ru1 = ru_90) : (ru1 = ru_90c);
	      (current_phase2 > 0) ? (ru2 = ru_90) : (ru2 = ru_90c);
	      // take absolute value and clip
	      current_phase1 = cmin(abs(current_phase1),127);
	      current_phase2 = cmin(abs(current_phase2),127);
	      
	      //          msg("sym: %d, current_phase1: %d, ru: %d + j%d, current_phase2: %d, ru: %d + j%d\n",k,current_phase1,ru1[2*current_phase1],ru1[2*current_phase1+1],current_phase2,ru2[2*current_phase2],ru2[2*current_phase2+1]);
	      
	      // rotate channel estimates by estimated phase
	      rotate_cpx_vector((int16_t*) ul_ch1,
				&ru1[2*current_phase1],
				(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],
				Msc_RS,
				15);
	      
	      rotate_cpx_vector((int16_t*) ul_ch2,
				&ru2[2*current_phase2],
				(int16_t*) &tmp_estimates[0],
				Msc_RS,
				15);
	      
	      // Combine the two rotated estimates
	      multadd_complex_vector_real_scalar((int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],SCALE,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],1,Msc_RS);
	      multadd_complex_vector_real_scalar((int16_t*) &tmp_estimates[0],SCALE,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],0,Msc_RS);
	      
	      /*
		if ((k<pilot_pos1) || ((k>pilot_pos2))) {
		
391
                multadd_complex_vector_real_scalar((int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],SCALE>>1,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],1,Msc_RS);
392
		
393
                multadd_complex_vector_real_scalar((int16_t*) &tmp_estimates[0],SCALE>>1,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],0,Msc_RS);
394
395
396
		
		} else {
		
397
                multadd_complex_vector_real_scalar((int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],SCALE>>1,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],1,Msc_RS);
398
		
399
                multadd_complex_vector_real_scalar((int16_t*) &tmp_estimates[0],SCALE>>1,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],0,Msc_RS);
400
		
401
                //              multadd_complex_vector_real_scalar((int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],alpha,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],1,Msc_RS);
402
		
403
                //              multadd_complex_vector_real_scalar((int16_t*) &tmp_estimates[0],beta ,(int16_t*) &ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],0,Msc_RS);
404
405
406
407
408
409
410
		
		}
	      */
	      
	      //      memcpy(&ul_ch_estimates[aa][frame_parms->N_RB_UL*12*k],ul_ch1,Msc_RS*sizeof(int32_t));
	    }
	  } //for(k=...
411

412
413
414
415
	  // because of the scaling of alpha and beta we also need to scale the final channel estimate at the pilot positions
	  
	  //    multadd_complex_vector_real_scalar((int16_t*) ul_ch1,SCALE,(int16_t*) ul_ch1,1,Msc_RS);
	  //    multadd_complex_vector_real_scalar((int16_t*) ul_ch2,SCALE,(int16_t*) ul_ch2,1,Msc_RS);
416

417
	} //if (Ns&1)
418
419

    } //for(aa=...
420

421
422
423
  } //if(l==...


424

425
  return(0);
426
}
427

gauthier's avatar
gauthier committed
428
extern uint16_t transmission_offset_tdd[16];
429
//#define DEBUG_SRS
430

gauthier's avatar
gauthier committed
431
int32_t lte_srs_channel_estimation(LTE_DL_FRAME_PARMS *frame_parms,
432
433
                                   LTE_eNB_COMMON *common_vars,
                                   LTE_eNB_SRS *srs_vars,
434
                                   SOUNDINGRS_UL_CONFIG_DEDICATED *soundingrs_ul_config_dedicated,
435
                                   unsigned char subframe,
436
                                   unsigned char eNB_id)
437
{
438

439
  int aa;
440
  int N_symb,symbol;
441
  uint8_t nb_antennas_rx = frame_parms->nb_antennas_rx;
442
443
444
445
#ifdef DEBUG_SRS
  char fname[40], vname[40];
#endif

446
447
  //uint8_t Ssrs  = frame_parms->soundingrs_ul_config_common.srs_SubframeConfig;
  //uint8_t T_SFC = (Ssrs<=7 ? 5 : 10);
448
449

  N_symb = 2*7-frame_parms->Ncp;
450
  symbol = N_symb-1; //SRS is always in last symbol of subframe
451

452
  /*
453
     msg("SRS channel estimation eNB %d, subframs %d, %d %d %d %d %d\n",eNB_id,sub_frame_number,
454
455
456
457
458
459
460
     SRS_parms->Csrs,
     SRS_parms->Bsrs,
     SRS_parms->kTC,
     SRS_parms->n_RRC,
     SRS_parms->Ssrs);
  */

461
  //if ((1<<(sub_frame_number%T_SFC))&transmission_offset_tdd[Ssrs]) {
462

463
464
465
466
467
  if (generate_srs(frame_parms,
		   soundingrs_ul_config_dedicated,
		   &srs_vars->srs[eNB_id],
		   0x7FFF,
		   subframe)==-1) {
Mongazon's avatar
Mongazon committed
468
      LOG_E(PHY,"lte_srs_channel_estimation: Error in generate_srs\n");
469
470
471
      return(-1);
    }

472
    for (aa=0; aa<nb_antennas_rx; aa++) {
473
#ifdef DEBUG_SRS
474
475
476
477
      LOG_E(PHY,"SRS channel estimation eNB %d, subframs %d, aarx %d, %p, %p, %p\n",eNB_id,sub_frame_number,aa,
	    &common_vars->rxdataF[aa][frame_parms->ofdm_symbol_size*symbol],
	    srs_vars->srs,
	    srs_vars->srs_ch_estimates[aa]);
478
#endif
479
      
Mongazon's avatar
Mongazon committed
480
481
      //LOG_M("eNB_rxF.m","rxF",&common_vars->rxdataF[0][aa][2*frame_parms->ofdm_symbol_size*symbol],2*(frame_parms->ofdm_symbol_size),2,1);
      //LOG_M("eNB_srs.m","srs_eNB",common_vars->srs,(frame_parms->ofdm_symbol_size),1,1);
482

483

484
      mult_cpx_conj_vector((int16_t*) &common_vars->rxdataF[aa][2*frame_parms->ofdm_symbol_size*symbol],
485
			   (int16_t*) srs_vars->srs,
486
			   (int16_t*) srs_vars->srs_ch_estimates[aa],
487
488
489
			   frame_parms->ofdm_symbol_size,
			   15,
			   0);
490
491

#ifdef DEBUG_SRS
492
493
      sprintf(fname,"srs_ch_est%d.m",aa);
      sprintf(vname,"srs_est%d",aa);
Mongazon's avatar
Mongazon committed
494
      LOG_M(fname,vname,srs_vars->srs_ch_estimates[aa],frame_parms->ofdm_symbol_size,1,1);
495
496
#endif
    }
497

498
499
  /*
    else {
500
    for (aa=0;aa<nb_antennas_rx;aa++)
501
    bzero(srs_vars->srs_ch_estimates[eNB_id][aa],frame_parms->ofdm_symbol_size*sizeof(int));
502
503
504
505
506
    }
  */
  return(0);
}

gauthier's avatar
gauthier committed
507
int16_t lte_ul_freq_offset_estimation(LTE_DL_FRAME_PARMS *frame_parms,
508
509
510
511
                                      int32_t *ul_ch_estimates,
                                      uint16_t nb_rb)
{

512
#if defined(__x86_64__) || defined(__i386__)
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
568
569
570
571
572
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
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
  int k, rb;
  int a_idx = 64;
  uint8_t conj_flag = 0;
  uint8_t output_shift;
  int pilot_pos1 = 3 - frame_parms->Ncp;
  int pilot_pos2 = 10 - 2*frame_parms->Ncp;
  __m128i *ul_ch1 = (__m128i*)&ul_ch_estimates[pilot_pos1*frame_parms->N_RB_UL*12];
  __m128i *ul_ch2 = (__m128i*)&ul_ch_estimates[pilot_pos2*frame_parms->N_RB_UL*12];
  int32_t avg[2];
  int16_t Ravg[2];
  Ravg[0]=0;
  Ravg[1]=0;
  int16_t iv, rv, phase_idx;
  __m128i avg128U1, avg128U2, R[3], mmtmpD0,mmtmpD1,mmtmpD2,mmtmpD3;

  // round(tan((pi/4)*[1:1:N]/N)*pow2(15))
  int16_t alpha[128] = {201, 402, 603, 804, 1006, 1207, 1408, 1610, 1811, 2013, 2215, 2417, 2619, 2822, 3024, 3227, 3431, 3634, 3838, 4042, 4246, 4450, 4655, 4861, 5066, 5272, 5479, 5686, 5893, 6101, 6309, 6518, 6727, 6937, 7147, 7358, 7570, 7782, 7995, 8208, 8422, 8637, 8852, 9068, 9285, 9503, 9721, 9940, 10160, 10381, 10603, 10825, 11049, 11273, 11498, 11725, 11952, 12180, 12410, 12640, 12872, 13104, 13338, 13573, 13809, 14046, 14285, 14525, 14766, 15009, 15253, 15498, 15745, 15993, 16243, 16494, 16747, 17001, 17257, 17515, 17774, 18035, 18298, 18563, 18829, 19098, 19368, 19640, 19915, 20191, 20470, 20750, 21033, 21318, 21605, 21895, 22187, 22481, 22778, 23078, 23380, 23685, 23992, 24302, 24615, 24931, 25250, 25572, 25897, 26226, 26557, 26892, 27230, 27572, 27917, 28266, 28618, 28975, 29335, 29699, 30067, 30440, 30817, 31198, 31583, 31973, 32368, 32767};

  // compute log2_maxh (output_shift)
  avg128U1 = _mm_setzero_si128();
  avg128U2 = _mm_setzero_si128();

  for (rb=0; rb<nb_rb; rb++) {
    avg128U1 = _mm_add_epi32(avg128U1,_mm_madd_epi16(ul_ch1[0],ul_ch1[0]));
    avg128U1 = _mm_add_epi32(avg128U1,_mm_madd_epi16(ul_ch1[1],ul_ch1[1]));
    avg128U1 = _mm_add_epi32(avg128U1,_mm_madd_epi16(ul_ch1[2],ul_ch1[2]));

    avg128U2 = _mm_add_epi32(avg128U2,_mm_madd_epi16(ul_ch2[0],ul_ch2[0]));
    avg128U2 = _mm_add_epi32(avg128U2,_mm_madd_epi16(ul_ch2[1],ul_ch2[1]));
    avg128U2 = _mm_add_epi32(avg128U2,_mm_madd_epi16(ul_ch2[2],ul_ch2[2]));

    ul_ch1+=3;
    ul_ch2+=3;
  }

  avg[0] = (((int*)&avg128U1)[0] +
            ((int*)&avg128U1)[1] +
            ((int*)&avg128U1)[2] +
            ((int*)&avg128U1)[3])/(nb_rb*12);

  avg[1] = (((int*)&avg128U2)[0] +
            ((int*)&avg128U2)[1] +
            ((int*)&avg128U2)[2] +
            ((int*)&avg128U2)[3])/(nb_rb*12);

  //    msg("avg0 = %d, avg1 = %d\n",avg[0],avg[1]);
  avg[0] = cmax(avg[0],avg[1]);
  avg[1] = log2_approx(avg[0]);
  output_shift = cmax(0,avg[1]-10);
  //output_shift  = (log2_approx(avg[0])/2)+ log2_approx(frame_parms->nb_antennas_rx-1)+1;
  //    msg("avg= %d, shift = %d\n",avg[0],output_shift);

  ul_ch1 = (__m128i*)&ul_ch_estimates[pilot_pos1*frame_parms->N_RB_UL*12];
  ul_ch2 = (__m128i*)&ul_ch_estimates[pilot_pos2*frame_parms->N_RB_UL*12];

  // correlate and average the 2 channel estimates ul_ch1*ul_ch2
  for (rb=0; rb<nb_rb; rb++) {
    mmtmpD0 = _mm_madd_epi16(ul_ch1[0],ul_ch2[0]);
    mmtmpD1 = _mm_shufflelo_epi16(ul_ch1[0],_MM_SHUFFLE(2,3,0,1));
    mmtmpD1 = _mm_shufflehi_epi16(mmtmpD1,_MM_SHUFFLE(2,3,0,1));
    mmtmpD1 = _mm_sign_epi16(mmtmpD1,*(__m128i*)&conjugate);
    mmtmpD1 = _mm_madd_epi16(mmtmpD1,ul_ch2[0]);
    mmtmpD0 = _mm_srai_epi32(mmtmpD0,output_shift);
    mmtmpD1 = _mm_srai_epi32(mmtmpD1,output_shift);
    mmtmpD2 = _mm_unpacklo_epi32(mmtmpD0,mmtmpD1);
    mmtmpD3 = _mm_unpackhi_epi32(mmtmpD0,mmtmpD1);
    R[0] = _mm_packs_epi32(mmtmpD2,mmtmpD3);

    mmtmpD0 = _mm_madd_epi16(ul_ch1[1],ul_ch2[1]);
    mmtmpD1 = _mm_shufflelo_epi16(ul_ch1[1],_MM_SHUFFLE(2,3,0,1));
    mmtmpD1 = _mm_shufflehi_epi16(mmtmpD1,_MM_SHUFFLE(2,3,0,1));
    mmtmpD1 = _mm_sign_epi16(mmtmpD1,*(__m128i*)&conjugate);
    mmtmpD1 = _mm_madd_epi16(mmtmpD1,ul_ch2[1]);
    mmtmpD0 = _mm_srai_epi32(mmtmpD0,output_shift);
    mmtmpD1 = _mm_srai_epi32(mmtmpD1,output_shift);
    mmtmpD2 = _mm_unpacklo_epi32(mmtmpD0,mmtmpD1);
    mmtmpD3 = _mm_unpackhi_epi32(mmtmpD0,mmtmpD1);
    R[1] = _mm_packs_epi32(mmtmpD2,mmtmpD3);

    mmtmpD0 = _mm_madd_epi16(ul_ch1[2],ul_ch2[2]);
    mmtmpD1 = _mm_shufflelo_epi16(ul_ch1[2],_MM_SHUFFLE(2,3,0,1));
    mmtmpD1 = _mm_shufflehi_epi16(mmtmpD1,_MM_SHUFFLE(2,3,0,1));
    mmtmpD1 = _mm_sign_epi16(mmtmpD1,*(__m128i*)&conjugate);
    mmtmpD1 = _mm_madd_epi16(mmtmpD1,ul_ch2[2]);
    mmtmpD0 = _mm_srai_epi32(mmtmpD0,output_shift);
    mmtmpD1 = _mm_srai_epi32(mmtmpD1,output_shift);
    mmtmpD2 = _mm_unpacklo_epi32(mmtmpD0,mmtmpD1);
    mmtmpD3 = _mm_unpackhi_epi32(mmtmpD0,mmtmpD1);
    R[2] = _mm_packs_epi32(mmtmpD2,mmtmpD3);

    R[0] = _mm_add_epi16(_mm_srai_epi16(R[0],1),_mm_srai_epi16(R[1],1));
    R[0] = _mm_add_epi16(_mm_srai_epi16(R[0],1),_mm_srai_epi16(R[2],1));

    Ravg[0] += (((short*)&R)[0] +
                ((short*)&R)[2] +
                ((short*)&R)[4] +
                ((short*)&R)[6])/(nb_rb*4);

    Ravg[1] += (((short*)&R)[1] +
                ((short*)&R)[3] +
                ((short*)&R)[5] +
                ((short*)&R)[7])/(nb_rb*4);

    ul_ch1+=3;
    ul_ch2+=3;
  }

  // phase estimation on Ravg
  //   Ravg[0] = 56;
  //   Ravg[1] = 0;
  rv = Ravg[0];
  iv = Ravg[1];

  if (iv<0)
    iv = -Ravg[1];

  if (rv<iv) {
    rv = iv;
    iv = Ravg[0];
    conj_flag = 1;
  }

  //   msg("rv = %d, iv = %d\n",rv,iv);
  //   msg("max_avg = %d, log2_approx = %d, shift = %d\n",avg[0], avg[1], output_shift);

  for (k=0; k<6; k++) {
    (iv<(((int32_t)(alpha[a_idx]*rv))>>15)) ? (a_idx -= 32>>k) : (a_idx += 32>>k);
  }

  (conj_flag==1) ? (phase_idx = 127-(a_idx>>1)) : (phase_idx = (a_idx>>1));

  if (Ravg[1]<0)
    phase_idx = -phase_idx;

  return(phase_idx);
648
649
650
#elif defined(__arm__)
  return(0);
#endif
651
}