Skip to content

O
OpenXG-RAN
Commit 866e1962 authored by lfarizav's avatar lfarizav
Browse files
Options

Ziggurat normal random number generator implemented for frequency domain simulations

parent 38f263cc
Hide whitespace changes
Inline Side-by-side
Showing with 89 additions and 74 deletions
openair1/SIMULATION/RF/rf.c
View file @ 866e1962
...@@ -442,8 +442,8 @@ clock_t start=clock();*/ ...@@ -442,8 +442,8 @@ clock_t start=clock();*/
else else
{ {
//rx128_gain_lin=mm_mul_set1_ps(rx_gain_lin); //rx128_gain_lin=mm_mul_set1_ps(rx_gain_lin);
gauss0_sqrt_NOW=sqrt_NOW*ziggurat(); gauss0_sqrt_NOW=sqrt_NOW*ziggurat(0.0,1.0);
gauss1_sqrt_NOW=sqrt_NOW*ziggurat(); gauss1_sqrt_NOW=sqrt_NOW*ziggurat(0.0,1.0);
rx128_re = _mm_loadu_pd(&r_re[a][2*i]);//r_re[a][i],r_re[a][i+1] rx128_re = _mm_loadu_pd(&r_re[a][2*i]);//r_re[a][i],r_re[a][i+1]
rx128_im = _mm_loadu_pd(&r_im[a][2*i]);//r_im[a][i],r_im[a][i+1] rx128_im = _mm_loadu_pd(&r_im[a][2*i]);//r_im[a][i],r_im[a][i+1]
rx128_gain_lin = _mm_set1_pd(rx_gain_lin); rx128_gain_lin = _mm_set1_pd(rx_gain_lin);
......
openair1/SIMULATION/TOOLS/defs.h
View file @ 866e1962
...@@ -436,13 +436,13 @@ the value \f$\mathrm{sgn}(u)i\f$. The search requires at most \f$Nbits-1\f$ com ...@@ -436,13 +436,13 @@ the value \f$\mathrm{sgn}(u)i\f$. The search requires at most \f$Nbits-1\f$ com
*/ */
int gauss(unsigned int *gauss_LUT,unsigned char Nbits); int gauss(unsigned int *gauss_LUT,unsigned char Nbits);
double zigguratdouble(unsigned int r, double sigma);
double gaussdouble(double,double); double gaussdouble(double,double);
void randominit(unsigned int seed_init); void randominit(unsigned int seed_init);
void setup_nor(); void table_nor(unsigned long seed);
double nfix(void);
double uniformrandom(void); double uniformrandom(void);
void uniformrandomSSE(__m128d *d1,__m128d *d2); void uniformrandomSSE(__m128d *d1,__m128d *d2);
double ziggurat(); double ziggurat(double mean, double variance);
int freq_channel(channel_desc_t *desc,uint16_t nb_rb, int16_t n_samples); int freq_channel(channel_desc_t *desc,uint16_t nb_rb, int16_t n_samples);
int freq_channel_prach(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples,int16_t prach_fmt,int16_t n_ra_prb); int freq_channel_prach(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples,int16_t prach_fmt,int16_t n_ra_prb);
int init_freq_channel(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples); int init_freq_channel(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples);
......
openair1/SIMULATION/TOOLS/rangen_double.c
View file @ 866e1962
...@@ -29,9 +29,6 @@ ...@@ -29,9 +29,6 @@
#include "defs.h" #include "defs.h"
static unsigned int seed, iy, ir[98]; static unsigned int seed, iy, ir[98];
static uint32_t kn[128];static double wn[128],fn[128];
static unsigned long iz,jz,jsr=123456789;
static long hz;
/* /*
@defgroup _uniformdouble @defgroup _uniformdouble
...@@ -106,16 +103,67 @@ void randominit(unsigned seed_init) ...@@ -106,16 +103,67 @@ void randominit(unsigned seed_init)
#endif #endif
#endif #endif
/*!\brief Uniform linear congruential random number generator on \f$[0,1)\f$. Returns a double-precision floating-point number.*/
double uniformrandom(void)
{
int j;
j=1 + 97.0*iy/mod;
iy=ir[j];
seed = a*seed; /* mod 2**32 */
ir[j] = seed;
return( (double) iy/mod );
}
/*void uniformrandomSSE(__m128d *d1,__m128d *d2)
{
int i,j;
__m128d u_rand128;
__m128i j128;
int j_array[]={(int) (1 + 9.103678167e-08*iy_array[0]),(int) (1 + 9.103678167e-08*iy_array[1]),(int) (1 + 9.103678167e-08*iy_array[2]),(int) (1 + 9.103678167e-08*iy_array[3])};
//j128=_mm_setr_epi32(1 + 391.0*iy_array[0]/mod,1 + 391.0*iy_array[1]/mod,1 + 391.0*iy_array[2]/mod,1 + 391.0*iy_array[3]/mod);
//j128=_mm_setr_epi32(391,391,391,391);
//j128=_mm_mul_epu32(j128,_mm_setr_epi32(iy_array[0],iy_array[1],iy_array[2],iy_array[3]));
//j128=_mm_slli_epi32(j128,32);
//j=1 + 97.0*iy/mod;
iy_array[0]=ir[j_array[0]];
iy_array[1]=ir[j_array[1]];
iy_array[2]=ir[j_array[2]];
iy_array[3]=ir[j_array[3]];
seed_array[0] = a*seed_array[0];
seed_array[1] = a*seed_array[1];
seed_array[2] = a*seed_array[2];
seed_array[3] = a*seed_array[3];
ir[j_array[0]] = seed_array[0];
ir[j_array[1]] = seed_array[1];
ir[j_array[2]] = seed_array[2];
ir[j_array[3]] = seed_array[3];
*d1=_mm_setr_pd (2*((double) iy_array[0]/mod)-1, 2*((double) iy_array[1])-1);
*d2=_mm_setr_pd (2*((double) iy_array[2]/mod)-1, 2*((double) iy_array[3])-1);
//return ((double) iy/mod );
return;
}*/
/*!\brief Ziggurat random number generator based on rejection sampling. It returns a pseudorandom normally distributed double number between -4.5 and 4.5*/
//Procedure to create tables for normal distribution kn,wn and fn //Procedure to create tables for normal distribution kn,wn and fn
#define SHR3 (jz=jsr, jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5),jz+jsr) #define SHR3 (jz=jsr, jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5),jz+jsr)
#define UNI (0.5+(signed) SHR3 * 0.2328306e-9) #define UNI (0.5+(signed) SHR3 * 0.2328306e-9)
#define NOR (hz=SHR3, iz=hz&127, (abs(hz)<kn[iz])? hz*wn[iz] : nfix()) #define NOR (hz=SHR3,iz=(hz&127),(abs(hz)<kn[iz])? hz*wn[iz] : nfix())
double nfix() static double wn[128],fn[128];
static unsigned long iz,jz,jsr=123456789,kn[128];
static int32_t hz;
//#define NOR (hz=SHR3, printf("hz %d\n",hz),sign=(hz&128)>>7,printf("sign %s\n",(sign)?"-":"+"),iz=hz&127,printf("iz %d\n",iz), (abs(hz)<kn[iz])? (sign)?(-1)*hz*wn[iz]:hz*wn[iz] : (sign)?(-1)*nfix():nfix())
double nfix(void)
{ {
const double r = 3.442620; static double x, y; const double r = 3.442620;
static double x, y;
long l;
double d;
for (;;) for (;;)
{ {
x=hz*wn[iz]; x=hz * wn[iz];
//printf("iz %d,x %e,l %d, hz %d,d %e wn %e\n",iz, x, l, hz, d, wn[iz]);
if (iz==0) if (iz==0)
{ {
do do
...@@ -124,19 +172,28 @@ double nfix() ...@@ -124,19 +172,28 @@ double nfix()
y = - log (UNI); y = - log (UNI);
} }
while (y+y < x*x); while (y+y < x*x);
/*printf("x*x %e,y+y %e\n",x*x,y+y);
printf("return 1: %e\n",(hz>0)? r+x : -r-x);*/
return (hz>0)? r+x : -r-x; return (hz>0)? r+x : -r-x;
} }
if (fn[iz]+UNI*(fn[iz-1]-fn[iz])<exp(-0.5*x*x)) if (fn[iz]+UNI*(fn[iz-1]-fn[iz])<exp(-0.5*x*x)){
//printf("return 2: %e\n",x);
return x; return x;
}
hz = SHR3; hz = SHR3;
iz = hz&127; iz = hz&127;
if (abs(hz) < kn[iz]) if (abs(hz) < kn[iz]){
return (hz*wn[iz]); /*printf("return 3: %e\n",(double)hz*wn[iz]);
printf("return 3: iz %d, hz %d, wn[%d] %e, hz*wz[%d] %e\n",iz,hz,iz,wn[iz],iz,wn[iz]*(double)hz);*/
return ((hz)*wn[iz]);
}
} }
} }
void setup_nor() void table_nor(unsigned long seed)
{ {
jsr=seed;
printf("seed is %d\n",seed);
double dn = 3.442619855899; double dn = 3.442619855899;
int i; int i;
const double m1 = 2147483648.0; const double m1 = 2147483648.0;
...@@ -146,21 +203,21 @@ void setup_nor() ...@@ -146,21 +203,21 @@ void setup_nor()
q = vn/exp(-0.5*dn*dn); q = vn/exp(-0.5*dn*dn);
kn[0] = (int)((dn/q)*m1); kn[0] = ((dn/q)*m1);
kn[1] = 0; kn[1] = 0;
wn[0] = (double) ( q / m1 ); wn[0] = ( q / m1 );
wn[127] = (double) ( dn / m1 ); wn[127] = ( dn / m1 );
fn[0] = 1.0; fn[0] = 1.0;
fn[127] = (double) ( exp ( - 0.5 * dn * dn ) ); fn[127] = ( exp ( - 0.5 * dn * dn ) );
for ( i = 126; 1 <= i; i-- ) for ( i = 126; 1 <= i; i-- )
{ {
dn = sqrt (-2.0 * log ( vn/dn + exp(-0.5*dn*dn))); dn = sqrt (-2.0 * log ( vn/dn + exp(-0.5*dn*dn)));
kn[i+1] = (int) ((dn / tn)*m1); kn[i+1] = ((dn / tn)*m1);
tn = dn; tn = dn;
fn[i] = (double) (exp (-0.5*dn*dn)); fn[i] = (exp (-0.5*dn*dn));
wn[i] = (double) (dn / m1); wn[i] = (dn / m1);
} }
for ( i = 0; i <= 127; i++ ){ for ( i = 0; i <= 127; i++ ){
printf("i %d: kn %d, fn %e, wn %e\n",i,kn[i],fn[i],wn[i]); printf("i %d: kn %d, fn %e, wn %e\n",i,kn[i],fn[i],wn[i]);
...@@ -168,63 +225,21 @@ void setup_nor() ...@@ -168,63 +225,21 @@ void setup_nor()
return; return;
} }
double ziggurat(double mean, double variance)
/*!\brief Uniform linear congruential random number generator on \f$[0,1)\f$. Returns a double-precision floating-point number.*/
double uniformrandom(void)
{
int j;
j=1 + 97.0*iy/mod;
iy=ir[j];
seed = a*seed; /* mod 2**32 */
ir[j] = seed;
return( (double) iy/mod );
}
/*void uniformrandomSSE(__m128d *d1,__m128d *d2)
{
int i,j;
__m128d u_rand128;
__m128i j128;
int j_array[]={(int) (1 + 9.103678167e-08*iy_array[0]),(int) (1 + 9.103678167e-08*iy_array[1]),(int) (1 + 9.103678167e-08*iy_array[2]),(int) (1 + 9.103678167e-08*iy_array[3])};
//j128=_mm_setr_epi32(1 + 391.0*iy_array[0]/mod,1 + 391.0*iy_array[1]/mod,1 + 391.0*iy_array[2]/mod,1 + 391.0*iy_array[3]/mod);
//j128=_mm_setr_epi32(391,391,391,391);
//j128=_mm_mul_epu32(j128,_mm_setr_epi32(iy_array[0],iy_array[1],iy_array[2],iy_array[3]));
//j128=_mm_slli_epi32(j128,32);
//j=1 + 97.0*iy/mod;
iy_array[0]=ir[j_array[0]];
iy_array[1]=ir[j_array[1]];
iy_array[2]=ir[j_array[2]];
iy_array[3]=ir[j_array[3]];
seed_array[0] = a*seed_array[0];
seed_array[1] = a*seed_array[1];
seed_array[2] = a*seed_array[2];
seed_array[3] = a*seed_array[3];
ir[j_array[0]] = seed_array[0];
ir[j_array[1]] = seed_array[1];
ir[j_array[2]] = seed_array[2];
ir[j_array[3]] = seed_array[3];
*d1=_mm_setr_pd (2*((double) iy_array[0]/mod)-1, 2*((double) iy_array[1])-1);
*d2=_mm_setr_pd (2*((double) iy_array[2]/mod)-1, 2*((double) iy_array[3])-1);
//return ((double) iy/mod );
return;
}*/
/*!\brief Ziggurat random number generator based on rejection sampling. It returns a pseudorandom normally distributed double number between -4.5 and 4.5*/
double ziggurat()
{ {
//double nor=NOR;
//printf("NOR %e\n",nor);
return NOR; return NOR;
} }
/* /*
@defgroup _gaussdouble Gaussian random number generator based on modified Box-Muller transformation. @defgroup _gaussdouble Gaussian random number generator based on modified Box-Muller transformation.
@ingroup numerical @ingroup numerical
*/ */
/*!\brief Gaussian random number generator based on modified Box-Muller transformation.Returns a double-precision floating-point number. */ /*!\brief Gaussian random number generator based on modified Box-Muller transformation.Returns a double-precision floating-point number. */
#define random_SSE //#define random_SSE
#ifdef random_SSE #ifdef random_SSE
double gaussdouble(double mean, double variance)//It is necessary to improve the function. However if we enable SSE the gain in time it is not too much. double gaussdouble(double mean, double variance)//It is necessary to improve the function.
{ {
static int iset=0; static int iset=0;
static double gset; static double gset;
...@@ -293,11 +308,11 @@ double gaussdouble(double mean, double variance) ...@@ -293,11 +308,11 @@ double gaussdouble(double mean, double variance)
do { do {
/*count++; /*count++;
clock_t start=clock();*/ clock_t start=clock();*/
v1 = 2.0*uniformrandom()-1.0; v1 = 2.0*UNI-1.0;
/*clock_t stop=clock(); /*clock_t stop=clock();
printf("UE_freq_channel time is %f s, AVERAGE time is %f s, count %d, sum %e\n",(float) (stop-start)/CLOCKS_PER_SEC,(float) (sum+stop-start)/(count*CLOCKS_PER_SEC),count,sum+stop-start); printf("UE_freq_channel time is %f s, AVERAGE time is %f s, count %d, sum %e\n",(float) (stop-start)/CLOCKS_PER_SEC,(float) (sum+stop-start)/(count*CLOCKS_PER_SEC),count,sum+stop-start);
sum=(sum+stop-start);*/ sum=(sum+stop-start);*/
v2 = 2.0*uniformrandom()-1.0; v2 = 2.0*UNI-1.0;
r = v1*v1+v2*v2; r = v1*v1+v2*v2;
//printf("Inside do: r %e\n",r); //printf("Inside do: r %e\n",r);
} while (r >= 1.0); } while (r >= 1.0);
......
targets/SIMU/USER/oaisim_functions.c
View file @ 866e1962
...@@ -990,12 +990,12 @@ void init_seed(uint8_t set_seed) ...@@ -990,12 +990,12 @@ void init_seed(uint8_t set_seed)
if(set_seed) { if(set_seed) {
randominit (oai_emulation.info.seed); randominit (oai_emulation.info.seed);
setup_nor();//Setup for the normal probability distribution function table_nor(oai_emulation.info.seed);//Setup for the normal probability distribution function
set_taus_seed (oai_emulation.info.seed); set_taus_seed (oai_emulation.info.seed);
} else { } else {
randominit (0); randominit (0);
setup_nor();//Setup for the normal probability distribution function table_nor(1234567891);//Setup for the normal probability distribution function. Always use a non-zero unsigned long argument.
set_taus_seed (0); set_taus_seed (0);
} }
} }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment