- Fix: avoid infinite loop on mulitcast read thread if no data to read

Pre-ci test passed

git-svn-id: http://svn.eurecom.fr/openair4G/trunk@4000 818b1a75-f10b-46b9-bf7c-635c3b92a50f
This commit is contained in:
Cedric Roux
2013-06-25 09:05:02 +00:00
commit 56c66c09aa
4111 changed files with 2133103 additions and 0 deletions
+232
View File
@@ -0,0 +1,232 @@
#include <math.h>
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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
void init_freq_channel(channel_desc_t *desc,u16 nb_rb,s16 n_samples) {
double delta_f,freq; // 90 kHz spacing
double delay;
s16 f;
u8 l;
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]);
}
}
}
void freq_channel(channel_desc_t *desc,u16 nb_rb,s16 n_samples) {
s16 f;
u8 aarx,aatx,l;
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;
}
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]);
}
}
}
}
}
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,
u16 nb_rb) {
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);
u16 f;
u8 aarx,aatx;
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,
u16 nb_rb) {
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);
u16 f;
u8 aarx,aatx;
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)));
}
u8 pbch_polynomial_degree;
double a[7];
void load_pbch_desc(FILE *pbch_file_fd) {
int i;
char dummy[25];
fscanf(pbch_file_fd,"%d",&pbch_polynomial_degree);
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++) {
fscanf(pbch_file_fd,"%s",dummy);
a[i] = strtod(dummy,NULL);
printf("%f ",a[i]);
}
printf("\n");
}
double pbch_bler(double sinr) {
int i;
double log10_bler=a[pbch_polynomial_degree];
double sinrpow=sinr;
// printf("log10bler %f\n",log10_bler);
if (sinr<-7.9)
return(1.0);
else if (sinr>=0.0)
return(1e-4);
for (i=1;i<=pbch_polynomial_degree;i++) {
// printf("sinrpow %f\n",sinrpow);
log10_bler += (a[pbch_polynomial_degree-i]*sinrpow);
sinrpow *= sinr;
// printf("log10bler %f\n",log10_bler);
}
return(pow(10.0,log10_bler));
}
+292
View File
@@ -0,0 +1,292 @@
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "defs.h"
#define PI 3.1416
#define Am 20
#define MCS_COUNT 24
#define MCL 70 /*minimum coupling loss (MCL) in dB*/
#define theta_3dB (65*PI/180)
enum sector {SEC1, SEC2, SEC3};
scenario_desc_t scenario;
void get_chan_desc(node_desc_t* node_tx, node_desc_t* node_rx, channel_desc_t *ch_desc, scenario_desc_t* scenario) {
double dist;
dist = sqrt(pow((node_tx->x - node_rx->x), 2) + pow((node_tx->y - node_rx->y), 2));
/* conversion of distance into KM 3gpp (36-942)*/
ch_desc->path_loss_dB = (128.1 + 37.6 * log10(dist/1000));
}
int main() {
int enb_count = 16;
int ue_count = 50;
double sect_angle[3]={0,2*PI/3,4*PI/3};
double gain_max;
double theta;
double min_path_loss = 0;
int att_enb_index;
node_desc_t *enb_data[enb_count];
node_desc_t *ue_data[ue_count];
channel_desc_t *ul_channel[ue_count][enb_count];
channel_desc_t *dl_channel[ue_count][enb_count];
int count;
int mcs;
int ue_index, enb_index;
int return_value;
int nb_rb = 25; //No. of resource blocks
double sinr[enb_count][2*nb_rb];
double sinr_eff[ue_count][MCS_COUNT];
double bler[ue_count][MCS_COUNT];
double gain_sec[3];
double thermal_noise;
double interference;
double coupling;
FILE *fp;
char buffer[100];
char *sinr_bler;
double tlu_sinr;
double tlu_bler;
int line_num;
char *file_name[]= {"bler_1.csv", "bler_2.csv", "bler_3.csv", "bler_4.csv", "bler_5.csv", "bler_6.csv", "bler_7.csv", "bler_8.csv",
"bler_9.csv", "bler_10.csv", "bler_11.csv", "bler_12.csv", "bler_13.csv", "bler_14.csv", "bler_15.csv", "bler_16.csv",
"bler_17.csv", "bler_18.csv", "bler_19.csv", "bler_20.csv", "bler_21.csv", "bler_22.csv"};
double beta[MCS_COUNT] = {0, 0, 0, 0, 0.9459960937499999, 1.2912109374999994, 1.0133789062499998, 1.000390625,
1.02392578125, 1.8595703124999998, 2.424389648437498, 2.3946533203124982, 2.5790039062499988,
2.4084960937499984, 2.782617187499999, 2.7868652343749996, 3.92099609375, 4.0392578125,
4.56109619140625, 5.03338623046875, 5.810888671875, 6.449108886718749};
double enb_position[][2] = {{1100,1100},{1100,2100},{1100,3100},{1100,4100},
{2100,1100},{2100,2100},{2100,3100},{2100,4100},
{3100,1100},{3100,2100},{3100,3100},{3100,4100},
{4100,1100},{4100,2100},{4100,3100},{4100,4100}};
double ue_position[][2] = {{3340,4740},{1500,620},{1780,4220},{1300,3540},{780,3100},
{1140,540},{1340,3660},{860,1220},{2700,2140},{3860,3060},
{3740,1060},{1700,3060},{2180,1620},{4420,1060},{1300,3340},
{3700,3180},{3780,540},{1700,4380},{4140,4740},{820,4380},
{3300,1540},{2100,1780},{1780,2260},{1940,2620},{1580,1700},
{1460,1940},{940,1340},{2100,3540},{1260,4340},{2940,4060},
{3980,940},{540,2220},{3060,2140},{4620,3940},{4260,2820},
{3860,3500},{4140,4140},{3900,3500},{1500,2140},{2620,3820},
{3420,2820},{1580,3940},{660,2100},{2740,1180},{2500,2500},
{3580,3580},{3740,3140},{3020,3020},{4340,4140},{980,4300}};
randominit(0);
/////////////////////////////////////////////////////////////////////////////////////////////////
int tabl_len=0;
double local_table[MCS_COUNT][9][9];
for (mcs = 5; mcs <= MCS_COUNT; mcs++) {
fp = fopen(file_name[mcs - 1],"r");
if (fp == NULL) {
printf("ERROR: Unable to open the file\n");
}
else {
fgets(buffer, 100, fp);
tabl_len=0;
while (!feof(fp)) {
sinr_bler = strtok(buffer, ";");
local_table[mcs-1][0][tabl_len] = atof(sinr_bler);
sinr_bler = strtok(NULL,";");
local_table[mcs-1][1][tabl_len] = atof(sinr_bler);
tabl_len++;
fgets(buffer, 100, fp);
}
fclose(fp);
}
printf("\n table for mcs %d\n",mcs);
for (tabl_len=0;tabl_len<9;tabl_len++)
printf("%lf %lf \n ",local_table[mcs-1][0][tabl_len],local_table[mcs-1][1][tabl_len]);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
for (enb_index = 0; enb_index < enb_count; enb_index++)
enb_data[enb_index] = (node_desc_t *)(malloc(sizeof(node_desc_t)));
for (ue_index = 0; ue_index < ue_count; ue_index++)
ue_data[ue_index] = (node_desc_t *)(malloc(sizeof(node_desc_t)));
for (enb_index = 0; enb_index < enb_count; enb_index++) {
enb_data[enb_index]->x = enb_position[enb_index][0];
enb_data[enb_index]->y = enb_position[enb_index][1];
enb_data[enb_index]->tx_power_dBm = 40;
enb_data[enb_index]->ant_gain_dBi = 15;
enb_data[enb_index]->rx_noise_level = 5; //value in db
enb_data[enb_index]->n_sectors = 3;
}
for (ue_index = 0; ue_index < ue_count; ue_index++) {
ue_data[ue_index]->x = ue_position[ue_index][0];
ue_data[ue_index]->y = ue_position[ue_index][1];
ue_data[ue_index]->phi_rad = 2 * PI;
ue_data[ue_index]->tx_power_dBm = 20;
ue_data[ue_index]->ant_gain_dBi = 0;
ue_data[ue_index]->rx_noise_level = 9; //value in db
}
for (ue_index = 0; ue_index < ue_count; ue_index++) {
min_path_loss = 10000;
for (enb_index = 0; enb_index < enb_count; enb_index++) {
ul_channel[ue_index][enb_index] = new_channel_desc_scm(1, 1, SCM_C, 7.68, 0, 0, 0);
dl_channel[ue_index][enb_index] = new_channel_desc_scm(1, 1, SCM_C, 7.68, 0, 0, 0);
//printf("ue %d enb %d\n", ue_index, enb_index);
/* Calculating the angle in the range -pi to pi from the slope */
//(ue_data[ue_index])->alpha_rad[enb_index] = (double)(atan2((ue_data[ue_index]->x - enb_data[enb_index]->x), (ue_data[ue_index]->y - enb_data[enb_index]->y)));
ue_data[ue_index]->alpha_rad[enb_index] = atan2((ue_data[ue_index]->x - enb_data[enb_index]->x), (ue_data[ue_index]->y - enb_data[enb_index]->y));
//printf("angle is tan %lf\n", ue_data[ue_index]->alpha_rad[enb_index]);
if ((ue_data[ue_index]->alpha_rad[enb_index]) < 0) {
ue_data[ue_index]->alpha_rad[enb_index] = 2*PI + ue_data[ue_index]->alpha_rad[enb_index];
//printf("angle in radians is %lf\n", ue_data[ue_index]->alpha_rad[enb_index]);
}
for(count = 0; count < enb_data[enb_index]->n_sectors; count++) {
theta = sect_angle[count] - ue_data[ue_index]->alpha_rad[enb_index];
gain_sec[count] = -(Am < (12 * pow((theta/theta_3dB),2)) ? Am : (12 * pow((theta/theta_3dB),2)));
}
/* gain = -min(Am , 12 * (theta/theta_3dB)^2) */
gain_max = (gain_sec[SEC1] > gain_sec[SEC2]) ? ((gain_sec[SEC1] > gain_sec[SEC3]) ? gain_sec[SEC1]:gain_sec[SEC3]) :
((gain_sec[SEC2] > gain_sec[SEC3]) ? gain_sec[SEC2]:gain_sec[SEC3]);
get_chan_desc(enb_data[enb_index], ue_data[ue_index], ul_channel[ue_index][enb_index], &scenario);
get_chan_desc(enb_data[enb_index], ue_data[ue_index], dl_channel[ue_index][enb_index], &scenario);
//printf("Path loss for link between ue %d and enb %d is %lf and gain is %lf \n", ue_index, enb_index, dl_channel[ue_index][enb_index]->path_loss_dB, gain_max);
if (dl_channel[ue_index][enb_index]->path_loss_dB < min_path_loss) {
min_path_loss = dl_channel[ue_index][enb_index]->path_loss_dB;
att_enb_index = enb_index;
}
//return_value = random_channel(ul_channel[ue_index][enb_index]);
return_value = random_channel(dl_channel[ue_index][enb_index]);
/* Thermal noise is calculated using 10log10(K*T*B) K = Boltzmann´s constant T = room temperature B = bandwidth */
/* Taken as constant for the time being since the BW is not changing */
thermal_noise = -105; //value in dBm
if (0 == return_value) {
//freq_channel(ul_channel[ue_index][enb_index], nb_rb);
freq_channel(dl_channel[ue_index][enb_index], nb_rb);
coupling = MCL > (dl_channel[ue_index][enb_index]->path_loss_dB-(enb_data[enb_index]->ant_gain_dBi + gain_max)) ?
MCL : (dl_channel[ue_index][enb_index]->path_loss_dB-(enb_data[enb_index]->ant_gain_dBi + gain_max));
//printf ("coupling factor is %lf\n", coupling);
for (count = 0; count < (2 * nb_rb); count++) {
sinr[enb_index][count] = enb_data[enb_index]->tx_power_dBm
- coupling
- (thermal_noise + ue_data[ue_index]->rx_noise_level)
+ 10 * log10 (pow(dl_channel[ue_index][enb_index]->chF[0][count].r, 2)
+ pow(dl_channel[ue_index][enb_index]->chF[0][count].i, 2));
//printf("Dl_link SNR for res. block %d is %lf\n", count, sinr[enb_index][count]);
}
}
}
for (count = 0; count < 2 * nb_rb; count++) {
interference = 0;
for (enb_index = 0; enb_index < enb_count; enb_index++) {
if (att_enb_index != enb_index) {
interference += pow(10, 0.1 * sinr[enb_index][count]);
}
}
sinr[att_enb_index][count] -= 10*log10(1 + interference);
//printf("***Dl_link SINR for res. block %d is %lf\n", count, sinr[att_enb_index][count]);
for (mcs = 5; mcs <= MCS_COUNT; mcs++) {
sinr_eff[ue_index][mcs-1] += exp(-(pow(10, (sinr[att_enb_index][count])/10))/beta[mcs-1]);
//printf("Effective snr %lf\n",sinr_eff[ue_index][mcs-1]);
//sinr_eff[ue_index][mcs] += exp(-(sinr[att_enb_index][count])/beta[mcs]);
}
}
for (mcs = 5; mcs <= MCS_COUNT; mcs++) {
//printf("mcs value %d \n",mcs);
//printf("beta value %lf \n",-beta[mcs-1]);
//printf("snr_eff value %lf \n",log(sinr_eff[ue_index][mcs-1]));
sinr_eff[ue_index][mcs-1] = -beta[mcs-1] *log((sinr_eff[ue_index][mcs-1])/(2*nb_rb));//
//printf("snr_eff value %lf \n",sinr_eff[ue_index][mcs-1]);
sinr_eff[ue_index][mcs-1] = 10 * log10(sinr_eff[ue_index][mcs-1]);
sinr_eff[ue_index][mcs-1] *= 10;
sinr_eff[ue_index][mcs-1] = floor(sinr_eff[ue_index][mcs-1]);
if ((int)sinr_eff[ue_index][mcs-1]%2) {
sinr_eff[ue_index][mcs-1] += 1;
}
sinr_eff[ue_index][mcs-1] /= 10;
//printf("Effective snr %lf \n",sinr_eff[ue_index][mcs-1]);
bler[ue_index][mcs-1] = 0;
/*line_num = 0;
fp = fopen(file_name[mcs - 1],"r");
if (fp == NULL) {
printf("ERROR: Unable to open the file\n");
}
else {
fgets(buffer, 100, fp);
while (!feof(fp)) {
line_num++;
sinr_bler = strtok(buffer, ";");
tlu_sinr = atof(sinr_bler);
sinr_bler = strtok(NULL,";");
tlu_bler = atof(sinr_bler);
if (1 == line_num) {
if (sinr_eff[ue_index][mcs-1] < tlu_sinr) {
bler[ue_index][mcs-1] = 1;
break;
}
}
if (sinr_eff[ue_index][mcs-1] == tlu_sinr) {
bler[ue_index][mcs-1] = tlu_bler;
}
fgets(buffer, 100, fp);
}
fclose(fp);*/
for (tabl_len=0;tabl_len<9;tabl_len++) {
if(tabl_len==0)
if (sinr_eff[ue_index][mcs-1] < local_table[mcs-1][0][tabl_len]) {
bler[ue_index][mcs-1] = 1;
break;
}
if (sinr_eff[ue_index][mcs-1] == local_table[mcs-1][0][tabl_len]) {
bler[ue_index][mcs-1] = local_table[mcs-1][1][tabl_len];
}
}
//printf("\n###Dl_link UE %d attached to eNB %d \n MCS %d effective SNR %lf BLER %lf", ue_index, att_enb_index, mcs,sinr_eff[ue_index][mcs-1],bler[ue_index][mcs-1]);
}
//printf("\n\n");
printf("\n Ue_ix enb_ix mcs5 mcs6 mcs7 mcs8 mcs9 mcs10 mcs11 mcs12 mcs13\
mcs14 mcs15 mcs16 mcs17 mcs18 mcs19 mcs20 mcs21 mcs22\n");
printf("SINR %4d %4d %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f\
%+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f\n",
ue_index, att_enb_index, sinr_eff[ue_index][4], sinr_eff[ue_index][5], sinr_eff[ue_index][6], sinr_eff[ue_index][7],
sinr_eff[ue_index][8], sinr_eff[ue_index][9], sinr_eff[ue_index][10], sinr_eff[ue_index][11], sinr_eff[ue_index][12],
sinr_eff[ue_index][13], sinr_eff[ue_index][14], sinr_eff[ue_index][15], sinr_eff[ue_index][16], sinr_eff[ue_index][17],
sinr_eff[ue_index][18], sinr_eff[ue_index][19], sinr_eff[ue_index][20], sinr_eff[ue_index][21], sinr_eff[ue_index][22]);
printf("BLER %4d %4d %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f\
%+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f %+4.2f\n",
ue_index, att_enb_index, bler[ue_index][4], bler[ue_index][5], bler[ue_index][6], bler[ue_index][7],
bler[ue_index][8], bler[ue_index][9], bler[ue_index][10], bler[ue_index][11], bler[ue_index][12],
bler[ue_index][13], bler[ue_index][14], bler[ue_index][15], bler[ue_index][16], bler[ue_index][17],
bler[ue_index][18], bler[ue_index][19], bler[ue_index][20], bler[ue_index][21], bler[ue_index][22]);
}
}
+67
View File
@@ -0,0 +1,67 @@
% this script generates the .h file for the correlation matrices of the
% SCM-C for the 4x4 case (with polarization) as well as the 2x2 case
% (without polarization)
alpha_beta = ...
[-0.4616 + 0.5439i 0.0225 - 0.0595i;
0.2806 + 0.6476i 0.0088 + 0.0602i;
-0.1136 - 0.6818i 0.0307 + 0.0555i;
0.6944 + 0.5043i -0.0244 - 0.0028i;
0.4072 + 0.5626i 0.0828 - 0.2378i;
-0.7753 + 0.1776i 0.4194 - 0.2429i];
alpha = alpha_beta(:,1);
beta = alpha_beta(:,2);
Gamma = ...
[0.5953 0.4047 0 0;
0.4047 0.5953 0 0;
0 0 0.5953 -0.4047;
0 0 -0.4047 0.5953;
0.6134 0.3866 0 0;
0.3866 0.6134 0 0;
0 0 0.6134 -0.3866;
0 0 -0.3866 0.6134;
0.6090 0.3910 0 0;
0.3910 0.6090 0 0;
0 0 0.6090 -0.3910;
0 0 -0.3910 0.6090;
0.6430 0.3570 0 0;
0.3570 0.6430 0 0;
0 0 0.6430 -0.3570;
0 0 -0.3570 0.6430;
0.6935 0.3065 0 0;
0.3065 0.6935 0 0;
0 0 0.6935 -0.3065;
0 0 -0.3065 0.6935;
0.7535 0.2465 0 0;
0.2465 0.7535 0 0;
0 0 0.7535 -0.2465;
0 0 -0.2465 0.7535];
Gamma = reshape(Gamma.',4,4,[]);
%%
for i=1:size(Gamma,3)
A = [1 conj(alpha(i)); alpha(i) 1];
B = [1 conj(beta(i)); beta(i) 1];
R44(:,:,i) = kron(kron(A,Gamma(:,:,i)),B);
R22(:,:,i) = kron(A,B);
R22_sqrt(:,:,i) = sqrtm(R22(:,:,i));
end
%%
R22_sqrt_int = zeros(8,4,6);
R22_sqrt_int(1:2:end,:,:) = real(R22_sqrt);
R22_sqrt_int(2:2:end,:,:) = imag(R22_sqrt);
%%
fid = fopen('scm_corrmat.h','w');
fprintf(fid,'double R22_sqrt[][] = {\n');
for i=1:size(Gamma,3)
fprintf(fid,'{');
fprintf(fid,'%f, ',R22_sqrt_int(:,:,i));
fprintf(fid,'\b\b},\n');
end
fprintf(fid,'};\n');
fclose(fid)
+333
View File
@@ -0,0 +1,333 @@
#ifndef __SIMULATION_TOOLS_DEFS_H__
#define __SIMULATION_TOOLS_DEFS_H__
#include "PHY/TOOLS/defs.h"
/** @defgroup _numerical_ Useful Numerical Functions
*@{
The present clause specifies several numerical functions for testing of digital communication systems.
-# Generation of Uniform Random Bits
-# Generation of Quantized Gaussian Random Variables
-# Generation of Floating-point Gaussian Random Variables
-# Generic Multipath Channel Generator
* @defgroup _channel_ Multipath channel generator
* @ingroup _numerical_
* @{
*/
#define NB_SAMPLES_CHANNEL_OFFSET 4
typedef struct {
///Number of tx antennas
u8 nb_tx;
///Number of rx antennas
u8 nb_rx;
///number of taps
u8 nb_taps;
///linear amplitudes of taps
double *amps;
///Delays of the taps in mus. length(delays)=nb_taps. Has to be between 0 and Td.
double *delays;
///length of impulse response. should be set to 11+2*bw*t_max
u8 channel_length;
///channel state vector. size(state) = nb_taps * (n_tx * n_rx);
struct complex **a;
///interpolated (sample-spaced) channel impulse response. size(ch) = (n_tx * n_rx) * channel_length. ATTENTION: the dimensions of ch are the transposed ones of a. This is to allow the use of BLAS when applying the correlation matrices to the state.
struct complex **ch;
///Sampled frequency response (90 kHz resolution)
struct complex **chF;
///Maximum path delay in mus.
double Td;
///Channel bandwidth in MHz.
double BW;
///Ricean factor of first tap wrt other taps (0..1, where 0 means AWGN and 1 means Rayleigh channel).
double ricean_factor;
///Angle of arrival of wavefront (in radians). For Ricean channel only. This assumes that both RX and TX have linear antenna arrays with lambda/2 antenna spacing. Furhter it is assumed that the arrays are parallel to each other and that they are far enough apart so that we can safely assume plane wave propagation.
double aoa;
///If set to 1, aoa is randomized according to a uniform random distribution
s8 random_aoa;
///in Hz. if >0 generate a channel with a Clarke's Doppler profile with a maximum Doppler bandwidth max_Doppler. CURRENTLY NOT IMPLEMENTED!
double max_Doppler;
///Square root of the full correlation matrix size(R_tx) = nb_taps * (n_tx * n_rx) * (n_tx * n_rx).
struct complex **R_sqrt;
///path loss including shadow fading in dB
double path_loss_dB;
///additional delay of channel in samples.
s32 channel_offset;
///This parameter (0...1) allows for simple 1st order temporal variation. 0 means a new channel every call, 1 means keep channel constant all the time
double forgetting_factor;
///needs to be set to 1 for the first call, 0 otherwise.
u8 first_run;
/// initial phase for frequency offset simulation
double ip;
/// number of paths taken by transmit signal
u16 nb_paths;
} channel_desc_t;
typedef struct {
/// Number of sectors (set to 1 in case of an omnidirectional antenna)
u8 n_sectors;
/// Antenna orientation for each sector (for non-omnidirectional antennas) in radians wrt north
double alpha_rad[3];
/// Antenna 3dB beam width (in radians)
double phi_rad;
/// Antenna gain (dBi)
double ant_gain_dBi;
/// Tx power (dBm)
double tx_power_dBm;
/// Rx noise level (dB)
double rx_noise_level;
///x coordinate (cartesian, in m)
double x;
///y coordinate (cartesian, in m)
double y;
///z coordinate (antenna height, in m)
double z;
/// direction of travel in radians wrt north
double direction_rad;
/// speed of node (m/s)
double speed;
} node_desc_t;
typedef enum {
rural=0,
urban,
indoor
} scenario_t;
typedef struct {
/// Scenario classifcation
scenario_t scenario;
/// Carrier frequency in Hz
double carrier_frequency;
/// Bandwidth (in Hz)
double bandwidth;
/// path loss at 0m distance in dB
double path_loss_0;
/// path loss exponent
double path_loss_exponent;
/// shadow fading standard deviation [dB] (assuming log-normal shadow fading with 0 mean)
double shadow_fading_std;
/// correlation distance of shadow fading
double shadow_fading_correlation_distance;
/// Shadowing correlation between cells
double shadow_fading_correlation_cells;
/// Shadowing correlation between sectors
double shadow_fading_correlation_sectors;
/// Rice factor???
/// Walls (penetration loss)
/// Nodes in the scenario
node_desc_t* nodes;
} scenario_desc_t;
typedef enum {
custom=0,
SCM_A,
SCM_B,
SCM_C,
SCM_D,
EPA,
EVA,
ETU,
MBSFN,
Rayleigh8,
Rayleigh1,
Rayleigh1_800,
Rayleigh1_corr,
Rayleigh1_anticorr,
Rice8,
Rice1,
Rice1_corr,
Rice1_anticorr,
AWGN
} SCM_t;
/**
\brief This routine initializes a new channel descriptor
\param nb_tx Number of TX antennas
\param nb_rx Number of RX antennas
\param nb_taps Number of taps
\param channel_length Length of the interpolated channel impulse response
\param amps Linear amplitudes of the taps (length(amps)=channel_length). The values should sum up to 1.
\param delays Delays of the taps. If delays==NULL the taps are assumed to be spaced equidistantly between 0 and t_max.
\param R_sqrt Channel correlation matrix. If R_sqrt==NULL, no channel correlation is applied.
\param Td Maximum path delay in mus.
\param BW Channel bandwidth in MHz.
\param ricean_factor Ricean factor applied to all taps.
\param aoa Anlge of arrival
\param forgetting_factor This parameter (0...1) allows for simple 1st order temporal variation
\param max_Doppler This is the maximum Doppler frequency for Jakes' Model
\param channel_offset This is a time delay to apply to channel
\param path_loss_dB This is the path loss in dB
\param random_aoa If set to 1, AoA of ricean component is randomized
*/
channel_desc_t *new_channel_desc(u8 nb_tx,u8 nb_rx, u8 nb_taps, u8 channel_length, double *amps, double* delays, struct complex** R_sqrt, double Td, double BW, double ricean_factor, double aoa, double forgetting_factor, double max_Doppler, s32 channel_offset, double path_loss_dB,u8 random_aoa);
channel_desc_t *new_channel_desc_scm(u8 nb_tx,
u8 nb_rx,
SCM_t channel_model,
double BW,
double forgetting_factor,
s32 channel_offset,
double path_loss_dB);
/** \fn void random_channel(channel_desc_t *desc)
\brief This routine generates a random channel response (time domain) according to a tapped delay line model.
\param desc Pointer to the channel descriptor
*/
int random_channel(channel_desc_t *desc);
/**\fn void multipath_channel(channel_desc_t *desc,
double **tx_sig_re,
double **tx_sig_im,
double **rx_sig_re,
double **rx_sig_im,
u32 length,
u8 keep_channel)
\brief This function generates and applys a random frequency selective random channel model.
@param desc Pointer to channel descriptor
@param tx_sig_re input signal (real component)
@param tx_sig_im input signal (imaginary component)
@param rx_sig_re output signal (real component)
@param rx_sig_im output signal (imaginary component)
@param length Length of input signal
@param keep_channel Set to 1 to keep channel constant for null-B/F
*/
void multipath_channel(channel_desc_t *desc,
double **tx_sig_re,
double **tx_sig_im,
double **rx_sig_re,
double **rx_sig_im,
u32 length,
u8 keep_channel);
/*
\fn 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,
u16 nb_rb)
\brief This function computes the average SINR over all frequency resources of the PBCH. It is used for PHY abstraction of the PBCH BLER
@param desc Pointer to channel descriptor of eNB
@param desc Pointer to channel descriptor of interfering eNB 1
@param desc Pointer to channel descriptor of interfering eNB 2
@param snr_dB SNR of eNB
@param snr_i1_dB SNR of interfering eNB 1
@param snr_i2_dB SNR of interfering eNB 2
@param nb_rb Number of RBs in system
*/
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,
u16 nb_rb);
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,
u16 nb_rb);
double pbch_bler(double sinr);
/**@}*/
/**
* @defgroup _taus_ Tausworthe Uniform Random Variable Generator
* @ingroup _numerical_
* @{
\fn inline unsigned int taus()
\brief Tausworthe Uniform Random Generator. This is based on the hardware implementation described in
Lee et al, "A Hardware Gaussian Noise Generator Usign the Box-Muller Method and its Error Analysis," IEEE Trans. on Computers, 2006.
*/
inline unsigned int taus(void);
/**
\fn void set_taus_seed(unsigned int seed_init)
\brief Sets the seed for the Tausworthe generator.
@param seed_init 0 means generate based on CPU time, otherwise provide the seed
*/
void set_taus_seed(unsigned int seed_init);
/**@} */
/** @defgroup _gauss_ Generation of Quantized Gaussian Random Variables
* @ingroup _numerical_
* @{
This set of routines are used to generate quantized (i.e. fixed-point) Gaussian random noise efficiently.
The use of these routines allows for rapid computer simulation of digital communication systems. The method
is based on a lookup-table of the quantized normal probability distribution. The routines assume that the
continuous-valued Gaussian random-variable,\f$x\f$ is quantized
to \f$N\f$ bits over the interval \f$[-L\sigma,L\sigma)\f$ where \f$N\f$ and \f$L\f$ control the precision
and range of the quantization. The
random variable, \f$l\in\{-2^{N-1},-2^{N-1}+1,\cdots,0,1,\cdots,2^{N-1}-1\}\f$ corresponds to the event,
\f$E_l =
\begin{cases}
x\in\left[-\infty,-L\sigma\right) & l=-2^{N-1}, \\
x\in\left[\frac{lL\sigma}{2^{N-1}},\frac{(l+1)L\sigma}{2^{N-1}}\right) & <l>-2^{N-1}, \\
x\in\left[L\sigma,\infty\right) & l>-2^{N-1},
\end{cases}\f$
which occurs with probability
\f$\Pr(E_l) =
\begin{cases}
\mathrm{erfc}(L) & l=-2^{N-1}, \\
\mathrm{erfc}(L) & l>-2^{N-1}, \\
\mathrm{erf}\left(\frac{lL}{2^{N-1}}\right) \mathrm{erfc}\left(\frac{(l-1)L}{2^{N-1}}\right)& l>-2^{N-1}.
\end{cases}\f$
*/
/** \fn unsigned int *generate_gauss_LUT(unsigned char Nbits,unsigned char L)
\brief This routine generates a Gaussian pdf lookup table (LUT). The table has \f$2^{\mathrm{Nbits}-1}\f$ entries which represent
the right half of the pdf. The data stored in position \f$i\f$ is actually the scaled cumulative probability distribution,
\f$2^{31}\mathrm{erf}\left(\frac{iL}{2^{N-1}}\right)\f$. This represents the average number of times that the random variable
falls in the interval \f$\left[0,\frac{i}{2^{N-1}}\right)\f$. This format allows for rapid conversion of uniform 32-bit
random variables to \f$N\f$-bit Gaussian random variables using binary search.
@see gauss
@param Nbits Number of bits for the output variable
@param L Number of standard deviations in range
*/
unsigned int *generate_gauss_LUT(unsigned char Nbits,unsigned char L);
/** \fn int gauss(unsigned int *gauss_LUT,unsigned char Nbits);
\brief This routine returns a zero-mean unit-variance Gaussian random variable.
Given a 32-bit uniform random variable,
\f$\mathrm{u}\f$ (from \ref _taus_, we first extract the sign and then search in the monotonically increasing Gaussian LUT for
the two entries \f$(i,i+1)\f$ for which
\f$ 2^{31}\mathrm{erf}\left(\frac{i}{2^{Nbits-1}}\right) < |u| \leq 2^{31}\mathrm{erf}\left(\frac{i+1}{2^{Nbits-1}}\right) \f$ and assign
the value \f$\mathrm{sgn}(u)i\f$. The search requires at most \f$Nbits-1\f$ comparisons.
@see generate_gauss_LUT
@see taus
@param gauss_LUT pointer to lookup-table
@param Nbits number of bits for output variable ( between 1 and 16)
*/
int gauss(unsigned int *gauss_LUT,unsigned char Nbits);
double gaussdouble(double,double);
void randominit(unsigned int seed_init);
double uniformrandom(void);
void freq_channel(channel_desc_t *desc,u16 nb_rb, s16 n_samples);
void init_freq_channel(channel_desc_t *desc,u16 nb_rb,s16 n_samples);
u8 multipath_channel_nosigconv(channel_desc_t *desc);
void multipath_tv_channel(channel_desc_t *desc,
double **tx_sig_re,
double **tx_sig_im,
double **rx_sig_re,
double **rx_sig_im,
u16 length,
u8 keep_channel);
/**@} */
/**@} */
#endif
+138
View File
@@ -0,0 +1,138 @@
#include <assert.h>
#include <stdlib.h>
#include <math.h>
#include "defs.h"
unsigned int *generate_gauss_LUT(unsigned char Nbits,
unsigned char L
)
{
unsigned int *LUT_ptr,i;
LUT_ptr = (unsigned int *)malloc((1<<(Nbits-1))*sizeof(int));
assert(LUT_ptr);
for (i=0;i<(1<<(Nbits-1));i++) {
LUT_ptr[i] = (unsigned int)((double)((unsigned int)(1<<31))*erf(i*L/(double)(1<<(Nbits-1))));
#ifdef LUTDEBUG
printf("pos %d : LUT_ptr[%d]=%x (%f)\n",i,i,LUT_ptr[i],(double)(erf(i*L/(double)(1<<(Nbits-1)))));
#endif //LUTDEBUG
}
return(LUT_ptr);
}
int gauss(unsigned int *gauss_LUT,
unsigned char Nbits
) {
unsigned int search_pos,step_size,u,tmp,tmpm1,tmpp1,s;
// Get a 32-bit uniform random-variable
u = taus();
#ifdef DEBUG
printf("u = %u\n",u);
#endif //DEBUG
// if it is larger than 2^31 (here negative), save the sign and rescale down to 31-bits.
s = u & 0x80000000;
u &= 0x7fffffff;
#ifdef DEBUG
printf("u = %x,s=%d\n",u,s);
#endif //DEBUG
search_pos = (1<<(Nbits-2)); // starting position of the binary search
step_size = search_pos;
do {
step_size >>= 1;
tmp = gauss_LUT[search_pos];
tmpm1 = gauss_LUT[search_pos-1];
tmpp1 = gauss_LUT[search_pos+1];
#ifdef DEBUG
printf("search_pos %d, step_size %d: t %x tm %x,tp %x\n",search_pos,step_size,tmp,tmpm1,tmpp1);
#endif //DEBUG
if (u <= tmp)
if (u >tmpm1)
return s==0 ? (search_pos-1) : 1-search_pos;
else
search_pos -= step_size;
else
if (u <= tmpp1)
return s==0 ? search_pos : - search_pos;
else
search_pos += step_size;
}
while (step_size > 0);
// If it gets here we're beyond the positive edge so return max
return s==0 ? (1<<(Nbits-1))-1 : 1-((1<<(Nbits-1)));
}
#ifdef GAUSSMAIN
#define Nhistbits 8
void main(int argc,char **argv) {
unsigned int *gauss_LUT_ptr,i;
unsigned int hist[(1<<Nhistbits)];
int gvar,maxg=0,ming=9999999,maxnum=0,L,Ntrials,Nbits;
double meang=0.0,varg=0.0;
if (argc < 4) {
printf("Not enough arguments: %s Nbits L Ntrials\n",argv[0]);
exit(-1);
}
Nbits = atoi(argv[1]);
L = atoi(argv[2]);
Ntrials = atoi(argv[3]);
set_taus_seed();
// Generate Gaussian LUT 12-bit quantization over 5 standard deviations
gauss_LUT_ptr = generate_gauss_LUT(Nbits,L);
for (i=0;i<(1<<Nhistbits);i++)
hist[i] = 0;
for (i=0;i<Ntrials;i++) {
gvar = gauss(gauss_LUT_ptr,Nbits);
if (gvar == ((1<<(Nbits-1))-1))
maxnum++;
maxg = gvar > maxg ? gvar : maxg;
ming = gvar < ming ? gvar : ming;
meang += (double)gvar/Ntrials;
varg += (double)gvar*gvar/Ntrials;
gvar += (1<<(Nbits-1))-1;
hist[(gvar/(1<<(Nbits-Nhistbits)))]++;
// printf("%d\n",gauss(gauss_LUT_ptr,Nbits));
}
printf("Tail probability = %e(%x)\n",2*erfc((double)L*gauss_LUT_ptr[(1<<(Nbits-1))-1]/(unsigned int)(1<<31)),gauss_LUT_ptr[(1<<(Nbits-1))-1]);
printf("max %d, min %d, mean %f, stddev %f, Pr(maxnum)=%e(%d)\n",maxg,ming,meang,sqrt(varg),(double)maxnum/Ntrials,maxnum);
// for (i=0;i<(1<<Nhistbits);i++)
// printf("%d : %u\n",i,hist[i]);
free(gauss_LUT_ptr);
}
#endif //GAUSSMAIN
@@ -0,0 +1,201 @@
//#include <string.h>
#include <math.h>
//#include <unistd.h>
//#include <stdlib.h>
#include "PHY/TOOLS/defs.h"
//#define TEST_DEBUG
short max_array(const short *p, unsigned int size)
{
short max = *p;
unsigned int i;
for (i = 1; i < size; ++i)
if (max < p[i])
max = p[i];
return max;
}
short min_array(const short *p, unsigned int size)
{
short min = *p;
unsigned int i;
for (i = 1; i < size; ++i)
if (min > p[i])
min = p[i];
return min;
}
// Note that 'arraySize = 8*((3*8*6144)+12)';
void dlsch_LLR_quant(const short *llr, const unsigned int arraySize, const short Mlevel, short *llr_quant)
{
unsigned int i, j;
short max_llr, min_llr;
short llr_interval;
float quant_step;
float *transLevelArray;
min_llr = min_array(llr, arraySize);
max_llr = max_array(llr, arraySize);
llr_interval = (max_llr - min_llr);
quant_step = (float)llr_interval / Mlevel;
if ((Mlevel%2) != 0){
printf("Mlevel should be mutiple of 2...\n");
exit(-1);
}
transLevelArray = (float *)malloc((Mlevel+1)*sizeof(float));
if (!transLevelArray){
printf("Cannot allocate memory for transLevelArray...!\n");
exit(-1);
}
for(j=0; j < Mlevel+1; j++){
transLevelArray[j] = min_llr + j*quant_step;
}
for (i=0; i < arraySize; i++){
for(j=0; j < Mlevel; j++){
if ((transLevelArray[j] <= llr[i]) && (llr[i] <= transLevelArray[j+1])){
llr_quant[i] = (short)(transLevelArray[j] + quant_step/2); // mid-points are selected;
break;
}
if (transLevelArray[j+1] <= llr[i]){ // the last term;
llr_quant[i] = (short)(transLevelArray[j] + quant_step/2); // mid-points are selected;
}
}
}
#ifdef TEST_DEBUG
printf("min_llr: %d : max_llr: %d \n", min_llr, max_llr);
printf("llr_interval: %d : quant_step: %f \n", llr_interval, quant_step);
printf("transLevelArray = [");
for(j=0; j < Mlevel+1; j++)
printf("%f ", transLevelArray[j]);
printf("] \n\n");
#endif // TEST_DEBUG
}
/*
void dlsch_LLR_quant2(const short *llr, const unsigned int arraySize, const short bitsPerLLR, short *llr_quant)
{
unsigned int i;
short max_llr, min_llr;
short peak_llr;
unsigned char shift;
min_llr = min_array(llr, arraySize);
max_llr = max_array(llr, arraySize);
peak_llr = max(abs(max_llr), abs(min_llr));
shift = log2_approx(peak_llr);
if (bitsPerLLR < 1){
printf("Bits per LLR should be greater than 0\n");
exit(-1);
}
for (i=0; i < 6072; i++){
llr_quant[i] = llr[i] >> max(0, shift-bitsPerLLR); // quantize the LLRs;
llr_quant[i] <<= max(0, shift-bitsPerLLR); // reconstruct the LLRs;
printf("[%d] llr: %d : llr_quant %d \n",i, llr[i], llr_quant[i]);
}
exit(-1);
}
void clearbits(short *llr, const unsigned short startbit, const unsigned short nbits)
{
unsigned short i;
for (i=0; i < nbits; i++){
*llr &= ~(1 << (startbit + i));
}
}
void dlsch_LLR_quant3(const short *llr, const unsigned int arraySize, const short bitsPerLLR, short *llr_quant)
{
unsigned int i;
unsigned char shift;
unsigned int abs_mean_llr = 0;
//mean of the LLRs;
for(i=0; i < 6072; i++){
abs_mean_llr += abs(llr[i]);
}
abs_mean_llr /= 6072;
shift = log2_approx(abs_mean_llr);
printf("abs_mean_llr=%d, shift = %d\n", abs_mean_llr, shift);
exit(-1);
for (i=0; i < 6072; i++){
llr_quant[i] = llr[i] >> max(0, shift-bitsPerLLR); // quantize the LLRs;
printf("[%d] llr: %d : ", i, llr[i]);
//mask the llr_quant values such that their values do not exceed 2^(bitsPerLLR);
clearbits(&llr_quant[i], bitsPerLLR, sizeof(short)-bitsPerLLR-1);
printf("llr_quant_mask: %d : ", llr_quant[i]);
llr_quant[i] <<= max(0, shift-bitsPerLLR); // reconstruct the LLRs;
printf("llr_quant: %d \n", llr_quant[i]);
}
}
*/
void dlsch_MRC_relay_LLR(const short *llr_quant, const unsigned int arraySize, short *llr_quant_sum){
unsigned int i;
for (i=0; i < arraySize; i++){
llr_quant_sum[i] += llr_quant[i];
}
}
#ifdef TEST_DEBUG
void test_llr_quant()
{
short test[10] = {-1800, -5446, 345, 243, 130, -2111, 433, 4210, -10, -134};
short channel_output[10];
unsigned int i;
dlsch_LLR_quant(test, 10, 2, channel_output);
for (i = 0; i < 10; i++){
printf("llr: %d : llr_quant %d \n",test[i], channel_output[i]);
}
}
void main() {
test_llr_quant();
}
#endif // TEST_DEBUG
@@ -0,0 +1,199 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "defs.h"
#include "SIMULATION/RF/defs.h"
//#define DEBUG_CH
/*
int init_multipath_channel(channel_desc_t *channel) {
int i;
if (channel->nb_taps<=0) {
msg("init_multipath_channel: nb_taps must be > 0\n");
return(-1);
}
if (channel->channel_length<=0) {
msg("init_multipath_channel: channel_length must be > 0\n");
return(-1);
}
if (channel->nb_tx<=0) {
msg("init_multipath_channel: nb_tx must be > 0\n");
return(-1);
}
if (channel->nb_rx<=0) {
msg("init_multipath_channel: nb_rx must be > 0\n");
return(-1);
}
channel->amps = (double*) malloc(channel->nb_taps*sizeof(double));
if (!channel->amps) {
msg("init_multipath_channel: cannot allocate amps\n");
return(-1);
}
channel->delays = (double*) malloc(channel->nb_taps*sizeof(double));
if (!channel->delays) {
msg("init_multipath_channel: cannot allocate delays\n");
return(-1);
}
channel->state = (struct complex **) malloc(channel->nb_tx*channel->nb_rx*sizeof(struct complex *));
if (!channel->state) {
msg("init_multipath_channel: cannot allocate state\n");
return(-1);
}
for (i=0; i<channel->nb_tx*channel->nb_rx; i++) {
channel->state[i] = (struct complex *) malloc(channel->nb_taps*sizeof(struct complex));
if (!channel->state[i]) {
msg("init_multipath_channel: cannot allocate state\n");
return(-1);
}
}
channel->ch = (struct complex **) malloc(channel->nb_tx*channel->nb_rx*sizeof(struct complex *));
if (!channel->ch) {
msg("init_multipath_channel: cannot allocate ch\n");
return(-1);
}
for (i=0; i<channel->nb_tx*channel->nb_rx; i++) {
channel->ch[i] = (struct complex *) malloc(channel->channel_length*sizeof(struct complex));
if (!channel->ch[i]) {
msg("init_multipath_channel: cannot allocate ch\n");
return(-1);
}
}
channel->R_tx_sqrt = (struct complex **) malloc(channel->nb_tx*channel->nb_tx*sizeof(struct complex *));
if (!channel->Rx_tx_sqrt) {
msg("init_multipath_channel: cannot allocate R_tx_sqrt\n");
return(-1);
}
for (i=0; i<channel->nb_tx*channel->nb_tx; i++) {
channel->R_tx_sqrt[i] = (struct complex *) malloc(channel->channel_length*sizeof(struct complex));
if (!channel->R_tx_sqrt[i]) {
msg("init_multipath_channel: cannot allocate R_tx_sqrt\n");
return(-1);
}
}
channel->R_rx_sqrt = (struct complex **) malloc(channel->nb_rx*channel->nb_rx*sizeof(struct complex *));
if (!channel->R_rx_sqrt) {
msg("init_multipath_channel: cannot allocate R_rx_sqrt\n");
return(-1);
}
for (i=0; i<channel->nb_rx*channel->nb_rx; i++) {
channel->R_rx_sqrt[i] = (struct complex *) malloc(channel->channel_length*sizeof(struct complex));
if (!channel->R_rx_sqrt[i]) {
msg("init_multipath_channel: cannot allocate R_rx_sqrt\n");
return(-1);
}
}
return(0);
}
int free_multipath_channel(channel_desc_t *channel) {
int i;
free(channel->amps);
free(channel->delays);
for (i=0; i<channel->nb_tx*channel->nb_rx; i++) {
free(channel->state[i] = (struct complex *) malloc(channel->nb_taps*sizeof(struct complex));
if (!channel->state[i]) {
msg("init_multipath_channel: cannot allocate state\n");
return(-1);
}
}
free(channel->state);
channel->ch = (struct complex **) malloc(channel->nb_tx*channel->nb_rx*sizeof(struct complex *));
if (!channel->ch) {
msg("init_multipath_channel: cannot allocate ch\n");
return(-1);
}
for (i=0; i<channel->nb_tx*channel->nb_rx; i++) {
channel->ch[i] = (struct complex *) malloc(channel->channel_length*sizeof(struct complex));
if (!channel->ch[i]) {
msg("init_multipath_channel: cannot allocate ch\n");
return(-1);
}
}
channel->R_tx_sqrt = (struct complex *) malloc(channel->nb_tx*channel->nb_tx*sizeof(struct complex));
if (!channel->Rx_tx_sqrt) {
msg("init_multipath_channel: cannot allocate Rx_tx_sqrt\n");
return(-1);
}
channel->R_rx_sqrt = (struct complex *) malloc(channel->nb_rx*channel->nb_rx*sizeof(struct complex));
if (!channel->R_rx_sqrt) {
msg("init_multipath_channel: cannot allocate R_rx_sqrt\n");
return(-1);
}
return(0);
}
*/
u8 multipath_channel_nosigconv(channel_desc_t *desc){
random_channel(desc);
return(1);
}
//#define DEBUG_CH
void multipath_channel(channel_desc_t *desc,
double **tx_sig_re,
double **tx_sig_im,
double **rx_sig_re,
double **rx_sig_im,
u32 length,
u8 keep_channel) {
int i,ii,j,l;
struct complex rx_tmp,tx;
double path_loss = pow(10,desc->path_loss_dB/20);
int dd;
dd = abs(desc->channel_offset);
#ifdef DEBUG_CH
printf("[CHANNEL] keep = %d : path_loss = %g (%f), nb_rx %d, nb_tx %d, dd %d, len %d \n",keep_channel,path_loss,desc->path_loss_dB,desc->nb_rx,desc->nb_tx,dd,desc->channel_length);
#endif
if (keep_channel) {
// do nothing - keep channel
} else {
random_channel(desc);
}
#ifdef DEBUG_CH
for (l = 0;l<(int)desc->channel_length;l++) {
printf("%p (%f,%f) ",desc->ch[0],desc->ch[0][l].x,desc->ch[0][l].y);
}
printf("\n");
#endif
for (i=0;i<((int)length-dd);i++) {
for (ii=0;ii<desc->nb_rx;ii++) {
rx_tmp.x = 0;
rx_tmp.y = 0;
for (j=0;j<desc->nb_tx;j++) {
for (l = 0;l<(int)desc->channel_length;l++) {
if ((i>=0) && (i-l)>=0) {
tx.x = tx_sig_re[j][i-l];
tx.y = tx_sig_im[j][i-l];
}
else {
tx.x =0;
tx.y =0;
}
rx_tmp.x += (tx.x * desc->ch[ii+(j*desc->nb_rx)][l].x) - (tx.y * desc->ch[ii+(j*desc->nb_rx)][l].y);
rx_tmp.y += (tx.y * desc->ch[ii+(j*desc->nb_rx)][l].x) + (tx.x * desc->ch[ii+(j*desc->nb_rx)][l].y);
} //l
} // j
rx_sig_re[ii][i+dd] = rx_tmp.x*path_loss;
rx_sig_im[ii][i+dd] = rx_tmp.y*path_loss;
/*
if ((ii==0)&&((i%32)==0)) {
printf("%p %p %f,%f => %e,%e\n",rx_sig_re[ii],rx_sig_im[ii],rx_tmp.x,rx_tmp.y,rx_sig_re[ii][i-dd],rx_sig_im[ii][i-dd]);
}
*/
//rx_sig_re[ii][i] = sqrt(.5)*(tx_sig_re[0][i] + tx_sig_re[1][i]);
//rx_sig_im[ii][i] = sqrt(.5)*(tx_sig_im[0][i] + tx_sig_im[1][i]);
} // ii
} // i
}
@@ -0,0 +1,186 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "defs.h"
#include "SIMULATION/RF/defs.h"
#include <complex.h>
void tv_channel(channel_desc_t *desc,double complex ***H,u16 length);
double frand_a_b(double a, double b);
void tv_conv(double complex **h, double complex *x, double complex *y, u16 nb_samples, u8 nb_taps, int delay);
void multipath_tv_channel(channel_desc_t *desc,
double **tx_sig_re,
double **tx_sig_im,
double **rx_sig_re,
double **rx_sig_im,
u16 length,
u8 keep_channel) {
double complex **tx,**rx,***H_t,*rx_temp;//, *tv_H_t;
double path_loss = pow(10,desc->path_loss_dB/20);
int i,j,k,dd;
dd = abs(desc->channel_offset);
#ifdef DEBUG_CH
printf("[TV CHANNEL] keep = %d : path_loss = %g (%f), nb_rx %d, nb_tx %d, dd %d, len %d max_doppler %d\n",keep_channel,path_loss,desc->path_loss_dB,desc->nb_rx,desc->nb_tx,dd,desc->channel_length,desc->max_Doppler);
#endif
tx = (double complex **)malloc(desc->nb_tx*sizeof(double complex));
rx = (double complex **)malloc(desc->nb_rx*sizeof(double complex));
H_t= (double complex ***) malloc(desc->nb_tx*desc->nb_rx*sizeof(double complex **));
// tv_H_t = (double complex *) malloc(length*sizeof(double complex));
rx_temp= (double complex *) calloc(length,sizeof(double complex));
for(i=0;i<desc->nb_tx;i++){
tx[i] = (double complex *)calloc(length,sizeof(double complex));
}
for(i=0;i<desc->nb_rx;i++){
rx[i] = (double complex *)calloc(length,sizeof(double complex));
}
for(i=0;i<desc->nb_tx*desc->nb_rx;i++){
H_t[i] = (double complex **) malloc(length*sizeof(double complex *));
for(j=0;j<length;j++) {
H_t[i][j] = (double complex *) calloc (desc->nb_taps,sizeof(double complex));
}
}
for (i=0;i<desc->nb_tx;i++){
for(j=0;j<length;j++) {
tx[i][j] = tx_sig_re[i][j] +I*tx_sig_im[i][j];
}
}
if (keep_channel) {
// do nothing - keep channel
} else {
tv_channel(desc,H_t,length);
}
for(i=0;i<desc->nb_rx;i++){
for(j=0;j<desc->nb_tx;j++){
tv_conv(H_t[i+(j*desc->nb_rx)],tx[j],rx_temp,length,desc->nb_taps,dd);
for(k=0;k<length;k++){
rx[i][k] += rx_temp[k];
}
}
}
for(i=0;i<desc->nb_rx;i++){
for(j=0;j<length;j++){
rx_sig_re[i][j] = creal(rx[i][j])*path_loss;
rx_sig_im[i][j] = cimag(rx[i][j])*path_loss;
}
}
/* for(k=0;k<length;k++) {
tv_H_t[k] = H_t[0][k][0];
}*/
for(i=0;i<desc->nb_tx;i++){
free(tx[i]);
}
free(tx);
for(i=0;i<desc->nb_rx;i++) {
free(rx[i]);
}
free(rx);
for(i=0;i<desc->nb_rx*desc->nb_tx;i++){
for(j=0;j<length;j++){
free(H_t[i][j]);
}
free(H_t[i]);
}
free(H_t);
free(rx_temp);
}
//TODO: make phi_rad a parameter of this function
void tv_channel(channel_desc_t *desc,double complex ***H,u16 length){
int i,j,p,l,k;
double *alpha,*phi_rad,pi=acos(-1),*w_Hz;
alpha = (double *)calloc(desc->nb_paths,sizeof(double));
phi_rad = (double *)calloc(desc->nb_paths,sizeof(double));
w_Hz = (double *)calloc(desc->nb_paths,sizeof(double));
for(i=0;i<desc->nb_paths;i++) {
w_Hz[i]=desc->max_Doppler*cos(frand_a_b(0,2*pi));
phi_rad[i]=frand_a_b(0,2*pi);
}
if(desc->ricean_factor == 1){
for(i=0;i<desc->nb_paths;i++) {
alpha[i]=1/sqrt(desc->nb_paths);
}
}
else {
alpha[0]=sqrt(desc->ricean_factor/(desc->ricean_factor+1));
for(i=1;i<desc->nb_paths;i++) {
alpha[i] = (1/sqrt(desc->nb_paths-1))*(sqrt(1/(desc->ricean_factor+1)));
}
}
/*
// This is the code when we only consider a SISO case
for(i=0;i<length;i++)
{
for(j=0;j<desc->nb_taps;j++)
{
for(p=0;p<desc->nb_paths;p++)
{
H[i][j] += sqrt(desc->amps[j]/2)*alpha[p]*cexp(-I*(2*pi*w_Hz[p]*i*(1/(desc->BW*1e6))+phi_rad[p]));
}
}
}
for(j=0;j<desc->nb_paths;j++)
{
phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->BW)+phi_rad[j],2*pi);
}
*/
// if MIMO
for (i=0;i<desc->nb_rx;i++){
for(j=0;j<desc->nb_tx;j++){
for(k=0;k<length;k++){
for(l=0;l<desc->nb_taps;l++){
H[i+(j*desc->nb_rx)][k][l] = 0;
for(p=0;p<desc->nb_paths;p++){
H[i+(j*desc->nb_rx)][k][l] += sqrt(desc->amps[l]/2)*alpha[p]*cexp(I*(2*pi*w_Hz[p]*k*(1/(desc->BW*1e6))+phi_rad[p]));
}
}
}
for(j=0;j<desc->nb_paths;j++){
phi_rad[j] = fmod(2*pi*w_Hz[j]*(length-1)*(1/desc->BW)+phi_rad[j],2*pi);
}
}
}
free(alpha);
free(w_Hz);
free(phi_rad);
}
// time varying convolution
void tv_conv(double complex **h, double complex *x, double complex *y, u16 nb_samples, u8 nb_taps, int dd) {
int i,j;
for(i=0;i<((int)nb_samples-dd);i++) {
for(j=0;j<nb_taps;j++) {
if(i>j)
y[i+dd] += creal(h[i][j])*creal(x[i-j])-cimag(h[i][j])*cimag(x[i-j]) + I*(creal(h[i][j])*cimag(x[i-j])+cimag(h[i][j])*creal(x[i-j]));
}
}
}
double frand_a_b(double a, double b){
return (rand()/(double)RAND_MAX)*(b-a)+a;
}
+878
View File
@@ -0,0 +1,878 @@
#include <math.h>
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "PHY/TOOLS/defs.h"
#include "defs.h"
#include "scm_corrmat.h"
#include "UTIL/LOG/log.h"
//#define DEBUG_CH
channel_desc_t *new_channel_desc(u8 nb_tx,
u8 nb_rx,
u8 nb_taps,
u8 channel_length,
double *amps,
double *delays,
struct complex** R_sqrt,
double Td,
double BW,
double ricean_factor,
double aoa,
double forgetting_factor,
double max_Doppler,
s32 channel_offset,
double path_loss_dB,
u8 random_aoa) {
channel_desc_t *chan_desc = (channel_desc_t *)malloc(sizeof(channel_desc_t));
u16 i,j;
double delta_tau;
LOG_I(OCM,"[CHANNEL] Getting new channel descriptor, nb_tx %d, nb_rx %d, nb_taps %d, channel_length %d\n",
nb_tx,nb_rx,nb_taps,channel_length);
chan_desc->nb_tx = nb_tx;
chan_desc->nb_rx = nb_rx;
chan_desc->nb_taps = nb_taps;
chan_desc->channel_length = channel_length;
chan_desc->amps = amps;
LOG_D(OCM,"[CHANNEL] Doing delays ...\n");
if (delays==NULL) {
chan_desc->delays = (double*) malloc(nb_taps*sizeof(double));
delta_tau = Td/nb_taps;
for (i=0; i<nb_taps; i++)
chan_desc->delays[i] = ((double)i)*delta_tau;
}
else
chan_desc->delays = delays;
chan_desc->Td = Td;
chan_desc->BW = BW;
chan_desc->ricean_factor = ricean_factor;
chan_desc->aoa = aoa;
chan_desc->random_aoa = random_aoa;
chan_desc->forgetting_factor = forgetting_factor;
chan_desc->channel_offset = channel_offset;
chan_desc->path_loss_dB = path_loss_dB;
chan_desc->first_run = 1;
chan_desc->ip = 0.0;
chan_desc->max_Doppler = max_Doppler;
chan_desc->ch = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->chF = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->a = (struct complex**) malloc(nb_taps*sizeof(struct complex*));
LOG_D(OCM,"[CHANNEL] Filling ch \n");
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->ch[i] = (struct complex*) malloc(channel_length * sizeof(struct complex));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex)); // allocate for up to 100 RBs, 12 samples per RB
LOG_D(OCM,"[CHANNEL] Filling a (nb_taps %d)\n",nb_taps);
for (i = 0; i<nb_taps; i++) {
LOG_D(OCM,"tap %d (%p,%d)\n",i,&chan_desc->a[i],nb_tx*nb_rx * sizeof(struct complex));
chan_desc->a[i] = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
}
LOG_D(OCM,"[CHANNEL] Doing R_sqrt ...\n");
if (R_sqrt == NULL) {
chan_desc->R_sqrt = (struct complex**) calloc(nb_taps,sizeof(struct complex*));
for (i = 0; i<nb_taps; i++) {
chan_desc->R_sqrt[i] = (struct complex*) calloc(nb_tx*nb_rx*nb_tx*nb_rx,sizeof(struct complex));
for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
chan_desc->R_sqrt[i][j].x = 1.0;
chan_desc->R_sqrt[i][j].y = 0.0;
}
}
}
else {
chan_desc->R_sqrt = R_sqrt;
}
for (i = 0; i<nb_taps; i++) {
for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
LOG_D(OCM,"Rsqrt[%d][%d] %f %f\n",i,j,chan_desc->R_sqrt[i][j].x,chan_desc->R_sqrt[i][j].y);
}
}
LOG_D(OCM,"[CHANNEL] RF %f\n",chan_desc->ricean_factor);
for (i=0;i<chan_desc->nb_taps;i++)
LOG_D(OCM,"[CHANNEL] tap %d: amp %f, delay %f\n",i,chan_desc->amps[i],chan_desc->delays[i]);
chan_desc->nb_paths=10;
return(chan_desc);
}
double mbsfn_delays[] = {0,.03,.15,.31,.37,1.09,12.490,12.52,12.64,12.80,12.86,13.58,27.49,27.52,27.64,27.80,27.86,28.58};
double mbsfn_amps_dB[] = {0,-1.5,-1.4,-3.6,-0.6,-7.0,-10,-11.5,-11.4,-13.6,-10.6,-17.0,-20,-21.5,-21.4,-23.6,-20.6,-27};
double scm_c_delays[] = {0, 0.0125, 0.0250, 0.3625, 0.3750, 0.3875, 0.2500, 0.2625, 0.2750, 1.0375, 1.0500, 1.0625, 2.7250, 2.7375, 2.7500, 4.6000, 4.6125, 4.6250};
double scm_c_amps_dB[] = {0.00, -2.22, -3.98, -1.86, -4.08, -5.84, -1.08, -3.30, -5.06, -9.08, -11.30, -13.06, -15.14, -17.36, -19.12, -20.64, -22.85, -24.62};
double epa_delays[] = { 0,.03,.07,.09,.11,.19,.41};
double epa_amps_dB[] = {0.0,-1.0,-2.0,-3.0,-8.0,-17.2,-20.8};
double eva_delays[] = { 0,.03,.15,.31,.37,.71,1.09,1.73,2.51};
double eva_amps_dB[] = {0.0,-1.5,-1.4,-3.6,-0.6,-9.1,-7.0,-12.0,-16.9};
double etu_delays[] = { 0,.05,.12,.2,.23,.5,1.6,2.3,5.0};
double etu_amps_dB[] = {-1.0,-1.0,-1.0,0.0,0.0,0.0,-3.0,-5.0,-7.0};
double default_amps_lin[] = {0.3868472 , 0.3094778 , 0.1547389 , 0.0773694 , 0.0386847 , 0.0193424 , 0.0096712 , 0.0038685};
double default_amp_lin[] = {1};
//correlation matrix for a 2x2 channel with full Tx correlation
struct complex R_sqrt_22_corr_tap[16] = {{0.70711,0}, {0.0, 0.0}, {0.70711,0}, {0.0, 0.0},
{0.0, 0.0}, {0.70711,0}, {0.0, 0.0}, {0.70711,0},
{0.70711,0}, {0.0, 0.0}, {0.70711,0}, {0.0, 0.0},
{0.0, 0.0}, {0.70711,0}, {0.0, 0.0}, {0.70711,0}};
struct complex *R_sqrt_22_corr[1] = {R_sqrt_22_corr_tap};
//correlation matrix for a fully correlated 2x1 channel (h1==h2)
struct complex R_sqrt_21_corr_tap[4] = {{0.70711,0}, {0.70711,0}, {0.70711,0}, {0.70711,0}};
struct complex *R_sqrt_21_corr[1] = {R_sqrt_21_corr_tap};
//correlation matrix for a 2x2 channel with full Tx anti-correlation
struct complex R_sqrt_22_anticorr_tap[16] = {{0.70711,0}, {0.0, 0.0}, {-0.70711,0}, {0.0, 0.0},
{0.0, 0.0}, {0.70711,0}, {0.0, 0.0}, {-0.70711,0},
{-0.70711,0}, {0.0, 0.0}, {0.70711,0}, {0.0, 0.0},
{0.0, 0.0}, {-0.70711,0}, {0.0, 0.0}, {0.70711,0}};
struct complex *R_sqrt_22_anticorr[1] = {R_sqrt_22_anticorr_tap};
//correlation matrix for a fully anti-correlated 2x1 channel (h1==-h2)
struct complex R_sqrt_21_anticorr_tap[4] = {{0.70711,0}, {-0.70711,0}, {-0.70711,0}, {0.70711,0}};
struct complex *R_sqrt_21_anticorr[1] = {R_sqrt_21_anticorr_tap};
struct complex **R_sqrt_ptr2;
channel_desc_t *new_channel_desc_scm(u8 nb_tx,
u8 nb_rx,
SCM_t channel_model,
double BW,
double forgetting_factor,
s32 channel_offset,
double path_loss_dB) {
channel_desc_t *chan_desc = (channel_desc_t *)malloc(sizeof(channel_desc_t));
u16 i,j;
double sum_amps;
double aoa,ricean_factor,Td,maxDoppler;
int channel_length,nb_taps;
chan_desc->nb_tx = nb_tx;
chan_desc->nb_rx = nb_rx;
chan_desc->BW = BW;
chan_desc->forgetting_factor = forgetting_factor;
chan_desc->channel_offset = channel_offset;
chan_desc->path_loss_dB = path_loss_dB;
chan_desc->first_run = 1;
chan_desc->ip = 0.0;
LOG_I(OCM,"Channel Model (inside of new_channel_desc_scm)=%d\n\n", channel_model);
switch (channel_model) {
case SCM_A:
LOG_W(OCM,"channel model not yet supported\n");
free(chan_desc);
return(NULL);
case SCM_B:
LOG_W(OCM,"channel model not yet supported\n");
free(chan_desc);
return(NULL);
case SCM_C:
chan_desc->nb_taps = 18;
chan_desc->Td = 4.625;
chan_desc->channel_length = (int) (2*chan_desc->BW*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->BW*chan_desc->Td));
sum_amps = 0;
chan_desc->amps = (double*) malloc(chan_desc->nb_taps*sizeof(double));
for (i = 0; i<chan_desc->nb_taps; i++) {
chan_desc->amps[i] = pow(10,.1*scm_c_amps_dB[i]);
sum_amps += chan_desc->amps[i];
}
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->amps[i] /= sum_amps;
chan_desc->delays = scm_c_delays;
chan_desc->ricean_factor = 1;
chan_desc->aoa = 0;
chan_desc->random_aoa = 0;
chan_desc->ch = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->chF = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->a = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->a[i] = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
if (nb_tx==2 && nb_rx==2) {
for (i = 0; i<6; i++)
chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
}
else if (nb_tx==2 && nb_rx==1) {
for (i = 0; i<6; i++)
chan_desc->R_sqrt[i] = (struct complex*) &R21_sqrt[i][0];
}
else if (nb_tx==1 && nb_rx==2) {
for (i = 0; i<6; i++)
chan_desc->R_sqrt[i] = (struct complex*) &R12_sqrt[i][0];
}
else {
for (i = 0; i<6; i++) {
chan_desc->R_sqrt[i] = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
chan_desc->R_sqrt[i][j].x = 1.0;
chan_desc->R_sqrt[i][j].y = 0.0;
}
LOG_W(OCM,"correlation matrix not implemented for nb_tx==%d and nb_rx==%d, using identity\n", nb_tx, nb_rx);
}
}
break;
case SCM_D:
LOG_W(OCM,"channel model not yet supported\n");
free(chan_desc);
return(NULL);
case EPA:
chan_desc->nb_taps = 7;
chan_desc->Td = .410;
chan_desc->channel_length = (int) (2*chan_desc->BW*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->BW*chan_desc->Td));
sum_amps = 0;
chan_desc->amps = (double*) malloc(chan_desc->nb_taps*sizeof(double));
for (i = 0; i<chan_desc->nb_taps; i++) {
chan_desc->amps[i] = pow(10,.1*epa_amps_dB[i]);
sum_amps += chan_desc->amps[i];
}
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->amps[i] /= sum_amps;
chan_desc->delays = epa_delays;
chan_desc->ricean_factor = 1;
chan_desc->aoa = 0;
chan_desc->random_aoa = 0;
chan_desc->ch = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->chF = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->a = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->a[i] = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
if (nb_tx==2 && nb_rx==2) {
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
for (i = 0; i<6; i++)
chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
}
else {
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
for (i = 0; i<6; i++) {
chan_desc->R_sqrt[i] = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
chan_desc->R_sqrt[i][j].x = 1.0;
chan_desc->R_sqrt[i][j].y = 0.0;
}
LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
}
}
break;
case EVA:
chan_desc->nb_taps = 9;
chan_desc->Td = 2.51;
chan_desc->channel_length = (int) (2*chan_desc->BW*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->BW*chan_desc->Td));
sum_amps = 0;
chan_desc->amps = (double*) malloc(chan_desc->nb_taps*sizeof(double));
for (i = 0; i<chan_desc->nb_taps; i++) {
chan_desc->amps[i] = pow(10,.1*eva_amps_dB[i]);
sum_amps += chan_desc->amps[i];
}
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->amps[i] /= sum_amps;
chan_desc->delays = eva_delays;
chan_desc->ricean_factor = 1;
chan_desc->aoa = 0;
chan_desc->random_aoa = 0;
chan_desc->ch = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->chF = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->a = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->a[i] = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
if (nb_tx==2 && nb_rx==2) {
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
for (i = 0; i<6; i++)
chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
}
else {
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
for (i = 0; i<6; i++) {
chan_desc->R_sqrt[i] = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
chan_desc->R_sqrt[i][j].x = 1.0;
chan_desc->R_sqrt[i][j].y = 0.0;
}
LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
}
}
break;
case ETU:
chan_desc->nb_taps = 9;
chan_desc->Td = 5.0;
chan_desc->channel_length = (int) (2*chan_desc->BW*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->BW*chan_desc->Td));
sum_amps = 0;
chan_desc->amps = (double*) malloc(chan_desc->nb_taps*sizeof(double));
for (i = 0; i<chan_desc->nb_taps; i++) {
chan_desc->amps[i] = pow(10,.1*etu_amps_dB[i]);
sum_amps += chan_desc->amps[i];
}
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->amps[i] /= sum_amps;
chan_desc->delays = etu_delays;
chan_desc->ricean_factor = 1;
chan_desc->aoa = 0;
chan_desc->random_aoa = 0;
chan_desc->ch = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->chF = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->a = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->a[i] = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
if (nb_tx==2 && nb_rx==2) {
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
for (i = 0; i<6; i++)
chan_desc->R_sqrt[i] = (struct complex*) &R22_sqrt[i][0];
}
else {
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
for (i = 0; i<6; i++) {
chan_desc->R_sqrt[i] = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
chan_desc->R_sqrt[i][j].x = 1.0;
chan_desc->R_sqrt[i][j].y = 0.0;
}
LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
}
}
break;
case MBSFN:
chan_desc->nb_taps = 18;
chan_desc->Td = 28.58;
chan_desc->channel_length = (int) (2*chan_desc->BW*chan_desc->Td + 1 + 2/(M_PI*M_PI)*log(4*M_PI*chan_desc->BW*chan_desc->Td));
sum_amps = 0;
chan_desc->amps = (double*) malloc(chan_desc->nb_taps*sizeof(double));
for (i = 0; i<chan_desc->nb_taps; i++) {
chan_desc->amps[i] = pow(10,.1*mbsfn_amps_dB[i]);
sum_amps += chan_desc->amps[i];
}
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->amps[i] /= sum_amps;
chan_desc->delays = mbsfn_delays;
chan_desc->ricean_factor = 1;
chan_desc->aoa = 0;
chan_desc->random_aoa = 0;
chan_desc->ch = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->chF = (struct complex**) malloc(nb_tx*nb_rx*sizeof(struct complex*));
chan_desc->a = (struct complex**) malloc(chan_desc->nb_taps*sizeof(struct complex*));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->ch[i] = (struct complex*) malloc(chan_desc->channel_length * sizeof(struct complex));
for (i = 0; i<nb_tx*nb_rx; i++)
chan_desc->chF[i] = (struct complex*) malloc(1200 * sizeof(struct complex));
for (i = 0; i<chan_desc->nb_taps; i++)
chan_desc->a[i] = (struct complex*) malloc(nb_tx*nb_rx * sizeof(struct complex));
chan_desc->R_sqrt = (struct complex**) malloc(6*sizeof(struct complex**));
for (i = 0; i<6; i++) {
chan_desc->R_sqrt[i] = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
chan_desc->R_sqrt[i][j].x = 1.0;
chan_desc->R_sqrt[i][j].y = 0.0;
}
LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
}
break;
case Rayleigh8:
nb_taps = 8;
Td = 0.8;
channel_length = (int)11+2*BW*Td;
ricean_factor = 1;
aoa = .03;
maxDoppler = 0;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amps_lin,
NULL,
NULL,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
0);
break;
case Rice8:
nb_taps = 8;
Td = 0.8;
channel_length = (int)11+2*BW*Td;
ricean_factor = 0.1;
aoa = .03;
maxDoppler = 0;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amps_lin,
NULL,
NULL,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
1);
break;
case Rayleigh1:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 1;
aoa = .03;
maxDoppler = 0;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
NULL,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
0);
break;
case Rayleigh1_800:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 1;
aoa = .03;
maxDoppler = 800;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
NULL,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
0);
break;
case Rayleigh1_corr:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 1;
aoa = .03;
maxDoppler = 0;
if ((nb_tx==2) && (nb_rx==1)) {
R_sqrt_ptr2 = R_sqrt_21_corr;
}
else if ((nb_tx==2) && (nb_rx==2)) {
R_sqrt_ptr2 = R_sqrt_22_corr;
}
else
R_sqrt_ptr2 = NULL;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
R_sqrt_ptr2,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
0);
break;
case Rayleigh1_anticorr:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 1;
aoa = .03;
maxDoppler = 0;
if ((nb_tx==2) && (nb_rx==1)) {
R_sqrt_ptr2 = R_sqrt_21_anticorr;
}
else if ((nb_tx==2) && (nb_rx==2)) {
R_sqrt_ptr2 = R_sqrt_22_anticorr;
}
else
R_sqrt_ptr2 = NULL;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
R_sqrt_ptr2,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
0);
break;
case Rice1:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 0.1;
aoa = .03;
maxDoppler = 0;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
NULL,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
0);
break;
case AWGN:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 0.0;
aoa = 0.0;
maxDoppler = 0;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
NULL,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
0);
printf("AWGN: ricean_factor %f\n",chan_desc->ricean_factor);
break;
case Rice1_corr:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 0.1;
aoa = .03;
maxDoppler = 0;
if ((nb_tx==2) && (nb_rx==1)) {
R_sqrt_ptr2 = R_sqrt_21_corr;
}
else if ((nb_tx==2) && (nb_rx==2)) {
R_sqrt_ptr2 = R_sqrt_22_corr;
}
else
R_sqrt_ptr2 = NULL;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
R_sqrt_ptr2,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
1);
break;
case Rice1_anticorr:
nb_taps = 1;
Td = 0;
channel_length = 1;
ricean_factor = 0.1;
aoa = .03;
maxDoppler = 0;
if ((nb_tx==2) && (nb_rx==1)) {
R_sqrt_ptr2 = R_sqrt_21_anticorr;
}
else if ((nb_tx==2) && (nb_rx==2)) {
R_sqrt_ptr2 = R_sqrt_22_anticorr;
}
else
R_sqrt_ptr2 = NULL;
chan_desc = new_channel_desc(nb_tx,
nb_rx,
nb_taps,
channel_length,
default_amp_lin,
NULL,
R_sqrt_ptr2,
Td,
BW,
ricean_factor,
aoa,
forgetting_factor,
maxDoppler,
channel_offset,
path_loss_dB,
1);
break;
default:
LOG_W(OCM,"channel model not yet supported\n");
free(chan_desc);
return(NULL);
}
LOG_D(OCM,"[CHANNEL] RF %f\n",chan_desc->ricean_factor);
for (i=0;i<chan_desc->nb_taps;i++)
LOG_D(OCM,"[CHANNEL] tap %d: amp %f, delay %f\n",i,chan_desc->amps[i],chan_desc->delays[i]);
chan_desc->nb_paths = 10;
return(chan_desc);
}
int random_channel(channel_desc_t *desc) {
double s;
int i,k,l,aarx,aatx;
struct complex anew[NB_ANTENNAS_TX*NB_ANTENNAS_RX],acorr[NB_ANTENNAS_TX*NB_ANTENNAS_RX];
struct complex phase, alpha, beta;
if ((desc->nb_tx>NB_ANTENNAS_TX) || (desc->nb_rx > NB_ANTENNAS_RX)) {
msg("random_channel.c: Error: temporary buffer for channel not big enough\n");
return(-1);
}
for (i=0;i<(int)desc->nb_taps;i++) {
for (aarx=0;aarx<desc->nb_rx;aarx++) {
for (aatx=0;aatx<desc->nb_tx;aatx++) {
anew[aarx+(aatx*desc->nb_rx)].x = sqrt(desc->ricean_factor*desc->amps[i]/2) * gaussdouble(0.0,1.0);
anew[aarx+(aatx*desc->nb_rx)].y = sqrt(desc->ricean_factor*desc->amps[i]/2) * gaussdouble(0.0,1.0);
if ((i==0) && (desc->ricean_factor != 1.0)) {
if (desc->random_aoa==1) {
desc->aoa = uniformrandom()*2*M_PI;
}
// this assumes that both RX and TX have linear antenna arrays with lambda/2 antenna spacing.
// Furhter it is assumed that the arrays are parallel to each other and that they are far enough apart so
// that we can safely assume plane wave propagation.
phase.x = cos(M_PI*((aarx-aatx)*sin(desc->aoa)));
phase.y = sin(M_PI*((aarx-aatx)*sin(desc->aoa)));
anew[aarx+(aatx*desc->nb_rx)].x += phase.x * sqrt(1.0-desc->ricean_factor);
anew[aarx+(aatx*desc->nb_rx)].y += phase.y * sqrt(1.0-desc->ricean_factor);
}
#ifdef DEBUG_CH
printf("(%d,%d,%d) %f->(%f,%f) (%f,%f) phase (%f,%f)\n",aarx,aatx,i,desc->amps[i],anew[aarx+(aatx*desc->nb_rx)].x,anew[aarx+(aatx*desc->nb_rx)].y,desc->aoa,desc->ricean_factor,phase.x,phase.y);
#endif
} //aatx
} //aarx
/*
// for debugging set a=anew;
for (aarx=0;aarx<desc->nb_rx;aarx++) {
for (aatx=0;aatx<desc->nb_tx;aatx++) {
desc->a[i][aarx+(aatx*desc->nb_rx)].x = anew[aarx+(aatx*desc->nb_rx)].x;
desc->a[i][aarx+(aatx*desc->nb_rx)].y = anew[aarx+(aatx*desc->nb_rx)].y;
printf("anew(%d,%d) = %f+1j*%f\n",aatx,aarx,anew[aarx+(aatx*desc->nb_rx)].x, anew[aarx+(aatx*desc->nb_rx)].y);
}
}
*/
//apply correlation matrix
//compute acorr = R_sqrt[i] * anew
alpha.x = 1.0;
alpha.y = 0.0;
beta.x = 0.0;
beta.y = 0.0;
cblas_zgemv(CblasRowMajor, CblasNoTrans, desc->nb_tx*desc->nb_rx, desc->nb_tx*desc->nb_rx,
(void*) &alpha, (void*) desc->R_sqrt[i/3], desc->nb_rx*desc->nb_tx,
(void*) anew, 1, (void*) &beta, (void*) acorr, 1);
/*
for (aarx=0;aarx<desc->nb_rx;aarx++) {
for (aatx=0;aatx<desc->nb_tx;aatx++) {
desc->a[i][aarx+(aatx*desc->nb_rx)].x = acorr[aarx+(aatx*desc->nb_rx)].x;
desc->a[i][aarx+(aatx*desc->nb_rx)].y = acorr[aarx+(aatx*desc->nb_rx)].y;
printf("tap %d, acorr1(%d,%d) = %f+1j*%f\n",i,aatx,aarx,acorr[aarx+(aatx*desc->nb_rx)].x, acorr[aarx+(aatx*desc->nb_rx)].y);
}
}
*/
if (desc->first_run==1){
cblas_zcopy(desc->nb_tx*desc->nb_rx, (void*) acorr, 1, (void*) desc->a[i], 1);
}
else {
// a = alpha*acorr+beta*a
// a = beta*a
// a = a+alpha*acorr
alpha.x = sqrt(1-desc->forgetting_factor);
alpha.y = 0;
beta.x = sqrt(desc->forgetting_factor);
beta.y = 0;
cblas_zscal(desc->nb_tx*desc->nb_rx, (void*) &beta, (void*) desc->a[i], 1);
cblas_zaxpy(desc->nb_tx*desc->nb_rx, (void*) &alpha, (void*) acorr, 1, (void*) desc->a[i], 1);
// desc->a[i][aarx+(aatx*desc->nb_rx)].x = (sqrt(desc->forgetting_factor)*desc->a[i][aarx+(aatx*desc->nb_rx)].x) + sqrt(1-desc->forgetting_factor)*anew.x;
// desc->a[i][aarx+(aatx*desc->nb_rx)].y = (sqrt(desc->forgetting_factor)*desc->a[i][aarx+(aatx*desc->nb_rx)].y) + sqrt(1-desc->forgetting_factor)*anew.y;
}
/*
for (aarx=0;aarx<desc->nb_rx;aarx++) {
for (aatx=0;aatx<desc->nb_tx;aatx++) {
//desc->a[i][aarx+(aatx*desc->nb_rx)].x = acorr[aarx+(aatx*desc->nb_rx)].x;
//desc->a[i][aarx+(aatx*desc->nb_rx)].y = acorr[aarx+(aatx*desc->nb_rx)].y;
printf("tap %d, a(%d,%d) = %f+1j*%f\n",i,aatx,aarx,desc->a[i][aarx+(aatx*desc->nb_rx)].x, desc->a[i][aarx+(aatx*desc->nb_rx)].y);
}
}
*/
} //nb_taps
//memset((void *)desc->ch[aarx+(aatx*desc->nb_rx)],0,(int)(desc->channel_length)*sizeof(struct complex));
for (aarx=0;aarx<desc->nb_rx;aarx++) {
for (aatx=0;aatx<desc->nb_tx;aatx++) {
if (desc->channel_length == 1) {
desc->ch[aarx+(aatx*desc->nb_rx)][0].x = desc->a[0][aarx+(aatx*desc->nb_rx)].x;
desc->ch[aarx+(aatx*desc->nb_rx)][0].y = desc->a[0][aarx+(aatx*desc->nb_rx)].y;
}
else {
for (k=0;k<(int)desc->channel_length;k++) {
desc->ch[aarx+(aatx*desc->nb_rx)][k].x = 0.0;
desc->ch[aarx+(aatx*desc->nb_rx)][k].y = 0.0;
for (l=0;l<desc->nb_taps;l++) {
if ((k - (desc->delays[l]*desc->BW) - NB_SAMPLES_CHANNEL_OFFSET) == 0)
s = 1.0;
else
s = sin(M_PI*(k - (desc->delays[l]*desc->BW) - NB_SAMPLES_CHANNEL_OFFSET))/
(M_PI*(k - (desc->delays[l]*desc->BW) - NB_SAMPLES_CHANNEL_OFFSET));
desc->ch[aarx+(aatx*desc->nb_rx)][k].x += s*desc->a[l][aarx+(aatx*desc->nb_rx)].x;
desc->ch[aarx+(aatx*desc->nb_rx)][k].y += s*desc->a[l][aarx+(aatx*desc->nb_rx)].y;
// printf("l %d : desc->ch.x %f\n",l,desc->a[l][aarx+(aatx*desc->nb_rx)].x);
} //nb_taps
#ifdef DEBUG_CH
printf("(%d,%d,%d)->(%f,%f)\n",k,aarx,aatx,desc->ch[aarx+(aatx*desc->nb_rx)][k].x,desc->ch[aarx+(aatx*desc->nb_rx)][k].y);
#endif
}
} //channel_length
} //aatx
} //aarx
if (desc->first_run==1)
desc->first_run = 0;
return (0);
}
#ifdef RANDOM_CHANNEL_MAIN
#define BW 5.0
#define Td 2.0
main(int argc,char **argv) {
double amps[8] = {.8,.2,.1,.04,.02,.01,.005};
struct complex ch[(int)(1+2*BW*Td)],phase;
int i;
randominit();
phase.x = 1.0;
phase.y = 0;
random_channel(amps,Td, 8,BW,ch,(double)1.0,&phase);
/*
for (i=0;i<(11+2*BW*Td);i++){
printf("%f + sqrt(-1)*%f\n",ch[i].x,ch[i].y);
}
*/
}
#endif
+104
View File
@@ -0,0 +1,104 @@
#ifdef USER_MODE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#endif
#include "defs.h"
static unsigned int seed, iy, ir[98];
/*
@defgroup _uniformdouble
@ingroup numerical Uniform linear congruential random number generator.
*/
/*!\brief Initialization routine for Uniform/Gaussian random number generators. */
#define a 1664525lu
#define mod 4294967296.0 /* is 2**32 */
void randominit(unsigned seed_init)
{
int i;
// this need to be integrated with the existing rng, like taus: navid
msg("Initializing random number generator, seed %x\n",seed_init);
if (seed_init == 0) {
srand((unsigned)time(NULL));
seed = (unsigned int) rand();
}
else {
seed = seed_init;
}
if (seed % 2 == 0) seed += 1; /* seed and mod are relative prime */
for (i=1; i<=97; i++)
{
seed = a*seed; /* mod 2**32 */
ir[i]= seed; /* initialize the shuffle table */
}
iy=1;
}
/*!\brief Uniform linear congruential random number generator on \f$[0,1)\f$. Returns a double-precision floating-point number.*/
double uniformrandom(void)
{
#define a 1664525lu
#define mod 4294967296.0 /* is 2**32 */
int j;
j=1 + 97.0*iy/mod;
iy=ir[j];
seed = a*seed; /* mod 2**32 */
ir[j] = seed;
return( (double) iy/mod );
}
/*
@defgroup _gaussdouble Gaussian random number generator based on modified Box-Muller transformation.
@ingroup numerical
*/
/*!\brief Gaussian random number generator based on modified Box-Muller transformation.Returns a double-precision floating-point number. */
double gaussdouble(double mean, double variance)
{
static int iset=0;
static double gset;
double fac,r,v1,v2;
if (iset == 0){
do{
v1 = 2.0*uniformrandom()-1.0;
v2 = 2.0*uniformrandom()-1.0;
r = v1*v1+v2*v2;
} while (r >= 1.0);
fac = sqrt(-2.0*log(r)/r);
gset= v1*fac;
iset=1;
return(sqrt(variance)*v2*fac + mean);
} else {
iset=0;
return(sqrt(variance)*gset + mean);
}
}
#ifdef MAIN
main(int argc,char **argv) {
int i;
randominit();
for (i=0;i<10;i++) {
printf("%f\n",gaussdouble(0.0,1.0));
}
}
#endif
+45
View File
@@ -0,0 +1,45 @@
alpha = [-0.4616 + 0.5439i, 0.2806 + 0.6476i, -0.1136 - 0.6818i, 0.6944 + 0.5043i, 0.4072 + 0.5626i, -0.7753 + 0.1776i];
beta = [ 0.0225 - 0.0595i, 0.0088 + 0.0602i, 0.0307 + 0.0555i, -0.0244 - 0.0028i, 0.0828 - 0.2378i, 0.4194 - 0.2429i];
n_tx = 2;
n_rx = 1;
A = zeros(length(alpha),n_tx,n_tx);
B = zeros(length(alpha),n_rx,n_rx);
C = zeros(length(alpha),n_tx*n_rx, n_tx*n_rx);
D = zeros(length(alpha),n_tx*n_rx, n_tx*n_rx);
fprintf(1,'double R%d%d_sqrt[%d][%d] = {\n', n_tx, n_rx, length(alpha), 2*(n_tx*n_rx)^2);
for i=1:length(alpha)
if (n_tx==1)
A(i,:,:) = 1;
elseif (n_tx==2)
A(i,:,:) = [1 conj(alpha(i)); alpha(i) 1];
else
error('n_tx must be <= 2')
end
if (n_rx==1)
B(i,:,:) = 1;
elseif (n_rx==2)
B(i,:,:) = [1 conj(beta(i)); beta(i) 1];
else
error('n_tx must be <= 2')
end
C(i,:,:) = kron(squeeze(A(i,:,:)),squeeze(B(i,:,:)));
D(i,:,:) = sqrtm(squeeze(C(i,:,:)));
fprintf(1,'{');
for j=1:((n_tx*n_rx)^2-1)
fprintf(1,'%f, %f,', real(D(i,j)), imag(D(i,j)));
end
j=(n_tx*n_rx)^2;
fprintf(1,'%f, %f},\n', real(D(i,j)), imag(D(i,j)));
end
fprintf(1,'};\n');
+24
View File
@@ -0,0 +1,24 @@
double R22_sqrt[6][32] = {
{0.921700, -0.000000, 0.010380, -0.027448, -0.250153, 0.294754, 0.005961, 0.010769, 0.010380, 0.027448, 0.921700, 0.000000, -0.011595, -0.004130, -0.250153, 0.294754, -0.250153, -0.294754, -0.011595, 0.004130, 0.921700, 0.000000, 0.010380, -0.027448, 0.005961, -0.010769, -0.250153, -0.294754, 0.010380, 0.027448, 0.921700, 0.000000},
{0.923810, 0.000000, 0.004069, 0.027832, 0.151730, 0.350180, -0.009882, 0.006114, 0.004069, -0.027832, 0.923810, 0.000000, 0.011218, -0.003029, 0.151730, 0.350180, 0.151730, -0.350180, 0.011218, 0.003029, 0.923810, -0.000000, 0.004069, 0.027832, -0.009882, -0.006114, 0.151730, -0.350180, 0.004069, -0.027832, 0.923810, 0.000000},
{0.927613, 0.000000, 0.014253, 0.025767, -0.061171, -0.367133, 0.009258, -0.007340, 0.014253, -0.025767, 0.927613, -0.000000, -0.011138, -0.003942, -0.061171, -0.367133, -0.061171, 0.367133, -0.011138, 0.003942, 0.927613, 0.000000, 0.014253, 0.025767, 0.009258, 0.007340, -0.061171, 0.367133, 0.014253, -0.025767, 0.927613, 0.000000},
{0.869794, -0.000000, -0.010613, -0.001218, 0.399115, 0.289852, -0.004464, -0.004096, -0.010613, 0.001218, 0.869794, -0.000000, -0.005276, -0.002978, 0.399115, 0.289852, 0.399115, -0.289852, -0.005276, 0.002978, 0.869794, -0.000000, -0.010613, -0.001218, -0.004464, 0.004096, 0.399115, -0.289852, -0.010613, 0.001218, 0.869794, 0.000000},
{0.919726, -0.000000, 0.038700, -0.111146, 0.217804, 0.300925, 0.045531, -0.013659, 0.038700, 0.111146, 0.919726, 0.000000, -0.027201, 0.038983, 0.217804, 0.300925, 0.217804, -0.300925, -0.027201, -0.038983, 0.919726, 0.000000, 0.038700, -0.111146, 0.045531, 0.013659, 0.217804, -0.300925, 0.038700, 0.111146, 0.919726, 0.000000},
{0.867608, -0.000000, 0.194097, -0.112414, -0.418811, 0.095938, -0.081264, 0.075727, 0.194097, 0.112414, 0.867608, -0.000000, -0.106125, -0.032801, -0.418811, 0.095938, -0.418811, -0.095938, -0.106125, 0.032801, 0.867608, 0.000000, 0.194097, -0.112414, -0.081264, -0.075727, -0.418811, -0.095938, 0.194097, 0.112414, 0.867608, 0.000000},
};
double R21_sqrt[6][8] = {
{0.922167, 0.000000,-0.250280, 0.294903,-0.250280, -0.294903,0.922167, 0.000000},
{0.924238, 0.000000,0.151801, 0.350342,0.151801, -0.350342,0.924238, 0.000000},
{0.928080, 0.000000,-0.061202, -0.367318,-0.061202, 0.367318,0.928080, 0.000000},
{0.869860, 0.000000,0.399145, 0.289874,0.399145, -0.289874,0.869860, 0.000000},
{0.927225, 0.000000,0.219580, 0.303378,0.219580, -0.303378,0.927225, 0.000000},
{0.896133, 0.000000,-0.432581, 0.099092,-0.432581, -0.099092,0.896133, 0.000000},
};
double R12_sqrt[6][8] = {
{0.999494, 0.000000,0.011256, -0.029765,0.011256, 0.029765,0.999494, 0.000000},
{0.999537, 0.000000,0.004402, 0.030114,0.004402, -0.030114,0.999537, 0.000000},
{0.999497, 0.000000,0.015358, 0.027764,0.015358, -0.027764,0.999497, 0.000000},
{0.999925, -0.000000,-0.012201, -0.001400,-0.012201, 0.001400,0.999925, 0.000000},
{0.991912, 0.000000,0.041738, -0.119870,0.041738, 0.119870,0.991912, 0.000000},
{0.968169, 0.000000,0.216594, -0.125443,0.216594, 0.125443,0.968169, 0.000000},
};
+81
View File
@@ -0,0 +1,81 @@
#ifdef USER_MODE
#include <time.h>
#include <stdlib.h>
#else
#include <asm/io.h>
#include <asm/rtai.h>
#endif
#ifdef RTAI_ENABLED
#include <rtai.h>
#include <rtai_sched.h>
#define time(x) (unsigned int)(rt_get_time_ns())
#endif
unsigned int s0, s1, s2, b;
//----------------------------------------------
//
/*!
*/
//
inline unsigned int taus(void) {
b = (((s0 << 13) ^ s0) >> 19);
s0 = (((s0 & 0xFFFFFFFE) << 12)^ b);
b = (((s1 << 2) ^ s1) >> 25);
s1 = (((s1 & 0xFFFFFFF8) << 4)^ b);
b = (((s2 << 3) ^ s2) >> 11);
s2 = (((s2 & 0xFFFFFFF0) << 17)^ b);
return s0 ^ s1 ^ s2;
}
void set_taus_seed(unsigned int seed_init) {
#ifdef USER_MODE
struct drand48_data buffer;
unsigned long result = 0;
#endif
if (seed_init == 0) {
s0 = (unsigned int)time(NULL);
s1 = (unsigned int)time(NULL);
s2 = (unsigned int)time(NULL);
}
else {
#ifndef USER_MODE
s0 = (unsigned int)0x1e23d852;
s1 = (unsigned int)0x81f38a1c;
s2 = (unsigned int)0xfe1a133e;
#else
/* Use reentrant version of rand48 to ensure that no conflicts with other generators occur */
srand48_r((long int)seed_init, &buffer);
mrand48_r(&buffer, (long int *)&result);
s0 = result;
mrand48_r(&buffer, (long int *)&result);
s1 = result;
mrand48_r(&buffer, (long int *)&result);
s2 = result;
#endif
}
}
#ifdef MAIN
main() {
unsigned int i,rand;
set_taus_seed();
for (i=0;i<10;i++){
rand = taus();
printf("%u\n",rand);
}
}
#endif //MAIN