abstraction.c 7.91 KB
Newer Older
1 2 3 4 5
#include <math.h>
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
6
#include <errno.h>
7 8 9 10 11 12 13 14 15 16 17 18 19

#include "PHY/TOOLS/defs.h"
#include "defs.h"

// NEW code with lookup table for sin/cos based on delay profile (TO BE TESTED)

double **cos_lut=NULL,**sin_lut=NULL;


//#if 1



20
void init_freq_channel(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples) {
21 22 23 24


  double delta_f,freq;  // 90 kHz spacing
  double delay;
25 26
  int16_t f;
  uint8_t l;
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56


  cos_lut = (double **)malloc(n_samples*sizeof(double*));
  sin_lut = (double **)malloc(n_samples*sizeof(double*));
  


  delta_f = nb_rb*180000/(n_samples-1);

  for (f=-(n_samples>>1);f<(n_samples>>1);f++) {
    freq=delta_f*(double)f*1e-6;// due to the fact that delays is in mus

    cos_lut[f+(n_samples>>1)] = (double *)malloc((int)desc->nb_taps*sizeof(double));
    sin_lut[f+(n_samples>>1)] = (double *)malloc((int)desc->nb_taps*sizeof(double));


    for (l=0;l<(int)desc->nb_taps;l++) {
      if (desc->nb_taps==1) 
	delay = desc->delays[l];
      else
	delay = desc->delays[l]+NB_SAMPLES_CHANNEL_OFFSET/desc->BW;

      cos_lut[f+(n_samples>>1)][l] = cos(2*M_PI*freq*delay);
      sin_lut[f+(n_samples>>1)][l] = sin(2*M_PI*freq*delay);
      //printf("values cos:%d, sin:%d\n", cos_lut[f][l], sin_lut[f][l]);
      
    }
  }
}

57
void freq_channel(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples) {
58 59


60 61
  int16_t f;
  uint8_t aarx,aatx,l;
62 63 64 65 66 67 68 69 70 71
  double *clut,*slut;
  static int freq_channel_init=0;

 // printf("no of samples:%d,",n_samples);
 // printf("no of taps:%d,",(int)desc->nb_taps);

  if (freq_channel_init == 0) {
    init_freq_channel(desc,nb_rb,n_samples);
    freq_channel_init=1;
  }
72 73
   
  start_meas(&desc->interp_freq);
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
  for (f=-n_samples/2;f<n_samples/2;f++) {
	clut = cos_lut[n_samples/2+f];
        slut = sin_lut[n_samples/2+f];
        
      for (aarx=0;aarx<desc->nb_rx;aarx++) {
	for (aatx=0;aatx<desc->nb_tx;aatx++) {
	  desc->chF[aarx+(aatx*desc->nb_rx)][n_samples/2+f].x=0.0;
	  desc->chF[aarx+(aatx*desc->nb_rx)][n_samples/2+f].y=0.0;
	  for (l=0;l<(int)desc->nb_taps;l++) {
		
	    desc->chF[aarx+(aatx*desc->nb_rx)][n_samples/2+f].x+=(desc->a[l][aarx+(aatx*desc->nb_rx)].x*clut[l]+
								  desc->a[l][aarx+(aatx*desc->nb_rx)].y*slut[l]);
	    desc->chF[aarx+(aatx*desc->nb_rx)][n_samples/2+f].y+=(-desc->a[l][aarx+(aatx*desc->nb_rx)].x*slut[l]+
								  desc->a[l][aarx+(aatx*desc->nb_rx)].y*clut[l]);
	  }
	}
      }
  }
92
  stop_meas(&desc->interp_freq);
93 94 95 96 97 98 99
}

double compute_pbch_sinr(channel_desc_t *desc,
			 channel_desc_t *desc_i1, 
			 channel_desc_t *desc_i2,
			 double snr_dB,double snr_i1_dB,
			 double snr_i2_dB,
100
			 uint16_t nb_rb) {
101 102

  double avg_sinr,snr=pow(10.0,.1*snr_dB),snr_i1=pow(10.0,.1*snr_i1_dB),snr_i2=pow(10.0,.1*snr_i2_dB);
103 104
  uint16_t f;
  uint8_t aarx,aatx;
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
  double S;
  struct complex S_i1;
  struct complex S_i2;

  avg_sinr=0.0;
  //  printf("nb_rb %d\n",nb_rb);
  for (f=(nb_rb-6);f<(nb_rb+6);f++) {
    S = 0.0;
    S_i1.x =0.0;
    S_i1.y =0.0;
    S_i2.x =0.0;
    S_i2.y =0.0;
    for (aarx=0;aarx<desc->nb_rx;aarx++) {
      for (aatx=0;aatx<desc->nb_tx;aatx++) {
	S    += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc->chF[aarx+(aatx*desc->nb_rx)][f].x + 
	         desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc->chF[aarx+(aatx*desc->nb_rx)][f].y);
	//	printf("%d %d chF[%d] => (%f,%f)\n",aarx,aatx,f,desc->chF[aarx+(aatx*desc->nb_rx)][f].x,desc->chF[aarx+(aatx*desc->nb_rx)][f].y);
	       
	if (desc_i1) {
	  S_i1.x += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].x + 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].y);
	  S_i1.y += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].y - 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].x);
	}
	if (desc_i2) {
	  S_i2.x += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].x + 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].y);
	  S_i2.y += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].y - 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].x);
	}
      }
    } 
    //    printf("snr %f f %d : S %f, S_i1 %f, S_i2 %f\n",snr,f-nb_rb,S,snr_i1*sqrt(S_i1.x*S_i1.x + S_i1.y*S_i1.y),snr_i2*sqrt(S_i2.x*S_i2.x + S_i2.y*S_i2.y));
    avg_sinr += (snr*S/(desc->nb_tx+snr_i1*sqrt(S_i1.x*S_i1.x + S_i1.y*S_i1.y)+snr_i2*sqrt(S_i2.x*S_i2.x + S_i2.y*S_i2.y)));
  }
  //  printf("avg_sinr %f (%f,%f,%f)\n",avg_sinr/12.0,snr,snr_i1,snr_i2);
  return(10*log10(avg_sinr/12.0));
}
 

double compute_sinr(channel_desc_t *desc,
		    channel_desc_t *desc_i1, 
		    channel_desc_t *desc_i2,
		    double snr_dB,double snr_i1_dB,
		    double snr_i2_dB,
150
		    uint16_t nb_rb) {
151 152

  double avg_sinr,snr=pow(10.0,.1*snr_dB),snr_i1=pow(10.0,.1*snr_i1_dB),snr_i2=pow(10.0,.1*snr_i2_dB);
153 154
  uint16_t f;
  uint8_t aarx,aatx;
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 188 189 190 191
  double S;
  struct complex S_i1;
  struct complex S_i2;

  avg_sinr=0.0;
  //  printf("nb_rb %d\n",nb_rb);
  for (f=0;f<2*nb_rb;f++) {
    S = 0.0;
    S_i1.x =0.0;
    S_i1.y =0.0;
    S_i2.x =0.0;
    S_i2.y =0.0;
    for (aarx=0;aarx<desc->nb_rx;aarx++) {
      for (aatx=0;aatx<desc->nb_tx;aatx++) {
	S    += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc->chF[aarx+(aatx*desc->nb_rx)][f].x + 
	         desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc->chF[aarx+(aatx*desc->nb_rx)][f].y);
	if (desc_i1) {
	  S_i1.x += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].x + 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].y);
	  S_i1.y += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].y - 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i1->chF[aarx+(aatx*desc->nb_rx)][f].x);
	}
	if (desc_i2) {
	  S_i2.x += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].x + 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].y);
	  S_i2.y += (desc->chF[aarx+(aatx*desc->nb_rx)][f].x*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].y - 
		     desc->chF[aarx+(aatx*desc->nb_rx)][f].y*desc_i2->chF[aarx+(aatx*desc->nb_rx)][f].x);
	}
      }
    }
    //        printf("f %d : S %f, S_i1 %f, S_i2 %f\n",f-nb_rb,snr*S,snr_i1*sqrt(S_i1.x*S_i1.x + S_i1.y*S_i1.y),snr_i2*sqrt(S_i2.x*S_i2.x + S_i2.y*S_i2.y));
    avg_sinr += (snr*S/(desc->nb_tx+snr_i1*sqrt(S_i1.x*S_i1.x + S_i1.y*S_i1.y)+snr_i2*sqrt(S_i2.x*S_i2.x + S_i2.y*S_i2.y)));
  }
  //  printf("avg_sinr %f (%f,%f,%f)\n",avg_sinr/12.0,snr,snr_i1,snr_i2);
  return(10*log10(avg_sinr/(nb_rb*2)));
}

192
int pbch_polynomial_degree;
193
double pbch_awgn_polynomial[7];
194 195 196

void load_pbch_desc(FILE *pbch_file_fd) {

197
  int i, ret;
198 199
  char dummy[25];

200 201 202 203 204
  ret = fscanf(pbch_file_fd,"%d",&pbch_polynomial_degree);
  if (ret < 0) {
    printf("fscanf failed: %s\n", strerror(errno));
    exit(EXIT_FAILURE);
  }
205 206 207 208 209 210 211 212
  if (pbch_polynomial_degree>6) {
    printf("Illegal degree for pbch interpolation polynomial %d\n",pbch_polynomial_degree);
    exit(-1);
  }

  printf("PBCH polynomial : ");

  for (i=0;i<=pbch_polynomial_degree;i++) {
213 214 215 216 217
    ret = fscanf(pbch_file_fd,"%s",dummy);
    if (ret < 0) {
      printf("fscanf failed: %s\n", strerror(errno));
      exit(EXIT_FAILURE);
    }
218 219
    pbch_awgn_polynomial[i] = strtod(dummy,NULL);
    printf("%f ",pbch_awgn_polynomial[i]);
220 221 222 223 224 225 226
  }
  printf("\n");
} 

double pbch_bler(double sinr) {

  int i;
227
  double log10_bler=pbch_awgn_polynomial[pbch_polynomial_degree];
228
  double sinrpow=sinr;
229
  double bler=0.0;
230 231
  //  printf("log10bler %f\n",log10_bler);
  if (sinr<-7.9)
232
    bler= 1.0;
233
  else if (sinr>=0.0)
234 235 236 237 238 239 240 241 242
    bler=0.0001;
  else  {
    for (i=1;i<=pbch_polynomial_degree;i++) {
      //    printf("sinrpow %f\n",sinrpow);
      log10_bler += (pbch_awgn_polynomial[pbch_polynomial_degree-i]*sinrpow);
      sinrpow *= sinr;
      //    printf("log10bler %f\n",log10_bler);
    }
    bler = pow(10.0,log10_bler);
243
  }
244 245 246
  printf ("sinr %f bler %f\n",sinr,bler);
  return(bler);

247
}