signal_energy.c 5.88 KB
Newer Older
ghaddab's avatar
ghaddab committed
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
/*******************************************************************************
    OpenAirInterface 
    Copyright(c) 1999 - 2014 Eurecom

    OpenAirInterface is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.


    OpenAirInterface is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with OpenAirInterface.The full GNU General Public License is 
   included in this distribution in the file called "COPYING". If not, 
   see <http://www.gnu.org/licenses/>.

  Contact Information
  OpenAirInterface Admin: openair_admin@eurecom.fr
  OpenAirInterface Tech : openair_tech@eurecom.fr
  OpenAirInterface Dev  : openair4g-devel@eurecom.fr
  
ghaddab's avatar
ghaddab committed
26
  Address      : Eurecom, Campus SophiaTech, 450 Route des Chappes, CS 50193 - 06904 Biot Sophia Antipolis cedex, FRANCE
ghaddab's avatar
ghaddab committed
27 28

 *******************************************************************************/
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
#include "defs.h"

#ifndef EXPRESSMIMO_TARGET
#include "mmintrin.h"
#endif //EXPRESSMIMO_TARGET

// Compute Energy of a complex signal vector, removing the DC component!
// input  : points to vector
// length : length of vector in complex samples

#define shift 4
#define shift_DC 0


#ifndef EXPRESSMIMO_TARGET
gauthier's avatar
gauthier committed
44
int32_t signal_energy(int32_t *input,uint32_t length) {
45

gauthier's avatar
gauthier committed
46 47
  int32_t i;
  int32_t temp,temp2;
48 49 50 51
  register __m64 mm0,mm1,mm2,mm3;
  __m64 *in = (__m64 *)input;

#ifdef MAIN
gauthier's avatar
gauthier committed
52
  int16_t *printb;
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
#endif

  mm0 = _m_pxor(mm0,mm0);
  mm3 = _m_pxor(mm3,mm3);

  for (i=0;i<length>>1;i++) {
    
    mm1 = in[i]; 
    mm2 = mm1;
    mm1 = _m_pmaddwd(mm1,mm1);
    mm1 = _m_psradi(mm1,shift);// shift any 32 bits blocs of the word by the value shift
    mm0 = _m_paddd(mm0,mm1);// add the two 64 bits words 4 bytes by 4 bytes
    //    temp2 = mm0;
    //    printf("%d %d\n",((int *)&temp2)[0],((int *)&temp2)[1]);


gauthier's avatar
gauthier committed
69
    //    printb = (int16_t *)&mm2;
70 71 72 73 74
    //    printf("mm2 %d : %d %d %d %d\n",i,printb[0],printb[1],printb[2],printb[3]);

    mm2 = _m_psrawi(mm2,shift_DC);
    mm3 = _m_paddw(mm3,mm2);// add the two 64 bits words 2 bytes by 2 bytes

gauthier's avatar
gauthier committed
75
    //    printb = (int16_t *)&mm3;
76 77 78 79 80 81
    //    printf("mm3 %d : %d %d %d %d\n",i,printb[0],printb[1],printb[2],printb[3]);

  }

  /*
#ifdef MAIN
gauthier's avatar
gauthier committed
82
  printb = (int16_t *)&mm3;
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
  printf("%d %d %d %d\n",printb[0],printb[1],printb[2],printb[3]);
#endif
  */
  mm1 = mm0;

  mm0 = _m_psrlqi(mm0,32);

  mm0 = _m_paddd(mm0,mm1);

  temp = _m_to_int(mm0);

  temp/=length;
  temp<<=shift;   // this is the average of x^2

  // now remove the DC component
  

  mm2 = _m_psrlqi(mm3,32);
  mm2 = _m_paddw(mm2,mm3);

  mm2 = _m_pmaddwd(mm2,mm2);

  temp2 = _m_to_int(mm2);

  temp2/=(length*length);

  temp2<<=(2*shift_DC);
#ifdef MAIN
  printf("E x^2 = %d\n",temp);  
#endif
  temp -= temp2;
#ifdef MAIN
  printf("(E x)^2=%d\n",temp2);
#endif
  _mm_empty();
  _m_empty();



  return((temp>0)?temp:1);
}

gauthier's avatar
gauthier committed
125
int32_t signal_energy_nodc(int32_t *input,uint32_t length) {
126

gauthier's avatar
gauthier committed
127 128
  int32_t i;
  int32_t temp;
129 130 131 132
  register __m64 mm0,mm1,mm2,mm3;
  __m64 *in = (__m64 *)input;

#ifdef MAIN
gauthier's avatar
gauthier committed
133
  int16_t *printb;
134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
#endif

  mm0 = _m_pxor(mm0,mm0);
  mm3 = _m_pxor(mm3,mm3);

  for (i=0;i<length>>1;i++) {
    
    mm1 = in[i]; 
    mm2 = mm1;
    mm1 = _m_pmaddwd(mm1,mm1);// SIMD complex multiplication
    mm1 = _m_psradi(mm1,shift);
    mm0 = _m_paddd(mm0,mm1);
    //    temp2 = mm0;
    //    printf("%d %d\n",((int *)&in[i])[0],((int *)&in[i])[1]);


gauthier's avatar
gauthier committed
150
    //    printb = (int16_t *)&mm2;
151 152 153 154 155 156 157
    //    printf("mm2 %d : %d %d %d %d\n",i,printb[0],printb[1],printb[2],printb[3]);


  }

  /*
#ifdef MAIN
gauthier's avatar
gauthier committed
158
  printb = (int16_t *)&mm3;
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
  printf("%d %d %d %d\n",printb[0],printb[1],printb[2],printb[3]);
#endif
  */
  mm1 = mm0;

  mm0 = _m_psrlqi(mm0,32);

  mm0 = _m_paddd(mm0,mm1);

  temp = _m_to_int(mm0);

  temp/=length;
  temp<<=shift;   // this is the average of x^2

#ifdef MAIN
  printf("E x^2 = %d\n",temp);  
#endif
  _mm_empty();
  _m_empty();



  return((temp>0)?temp:1);
}

gauthier's avatar
gauthier committed
184
double signal_energy_fp(double **s_re,double **s_im,uint32_t nb_antennas,uint32_t length,uint32_t offset) {
185

gauthier's avatar
gauthier committed
186
  int32_t aa,i;
187 188 189 190 191 192 193 194 195 196
  double V=0.0;

  for (i=0;i<length;i++) {
    for (aa=0;aa<nb_antennas;aa++) {
      V= V + (s_re[aa][i+offset]*s_re[aa][i+offset]) + (s_im[aa][i+offset]*s_im[aa][i+offset]); 
    }
  }
  return(V/length/nb_antennas);
}

gauthier's avatar
gauthier committed
197
double signal_energy_fp2(struct complex *s,uint32_t length) {
198

gauthier's avatar
gauthier committed
199
  int32_t i;
200 201 202 203 204 205 206 207 208 209 210
  double V=0.0;

  for (i=0;i<length;i++) {
    //    printf("signal_energy_fp2 : %f,%f => %f\n",s[i].x,s[i].y,V);
    //      V= V + (s[i].y*s[i].x) + (s[i].y*s[i].x); 
      V= V + (s[i].x*s[i].x) + (s[i].y*s[i].y); 
    }
  return(V/length);
}
#else

gauthier's avatar
gauthier committed
211
int32_t signal_energy(int32_t *input,uint32_t length) {
212 213 214 215 216 217 218 219 220 221 222 223 224 225
}

#endif

#ifdef MAIN
#define LENGTH 256
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
main(int argc,char **argv) {

  int input[LENGTH];
  int energy=0,dc_r=0,dc_i=0;
gauthier's avatar
gauthier committed
226
  int16_t s=1,i;
227 228 229 230 231 232 233 234
  int amp;

  amp = atoi(argv[1]);// arguments to integer
  if (argc>1)
    printf("Amp = %d\n",amp);

  for (i=0;i<LENGTH;i++) {
    s = -s;
gauthier's avatar
gauthier committed
235 236 237 238 239
    ((int16_t*)input)[2*i]     = 31 + (int16_t)(amp*sin(2*M_PI*i/LENGTH));
    ((int16_t*)input)[1+(2*i)] = 30 + (int16_t)(amp*cos(2*M_PI*i/LENGTH));
    energy += (((int16_t*)input)[2*i]*((int16_t*)input)[2*i]) + (((int16_t*)input)[1+(2*i)]*((int16_t*)input)[1+(2*i)]);
    dc_r += ((int16_t*)input)[2*i];
    dc_i += ((int16_t*)input)[1+(2*i)];
240 241 242 243 244 245 246 247 248 249 250 251


  }
  energy/=LENGTH;
  dc_r/=LENGTH;
  dc_i/=LENGTH;

  printf("signal_energy = %d dB(%d,%d,%d,%d)\n",dB_fixed(signal_energy(input,LENGTH)),signal_energy(input,LENGTH),energy-(dc_r*dc_r+dc_i*dc_i),energy,(dc_r*dc_r+dc_i*dc_i));
  printf("dc = (%d,%d)\n",dc_r,dc_i);
}
#endif