1 /*
2 * Copyright 2018 Oticon A/S
3 *
4 * SPDX-License-Identifier: Apache-2.0
5 */
6 #include <stdlib.h>
7 #include <math.h>
8 #include <complex.h>
9 #include "bs_tracing.h"
10 #include "bs_types.h"
11
12 /**
13 * Initialize the random generators
14 */
bs_random_init(unsigned int seed)15 int bs_random_init(unsigned int seed){
16 srandom(seed);
17 srand(seed);
18 return 0;
19 }
20
21 /**
22 * Free any memory this module may have allocated
23 */
bs_random_free()24 void bs_random_free(){
25
26 }
27
28 /**
29 * Generate a Gausian random number N(0,1): with mean 0 and variance 1
30 * In general to get a N(ave,std) do:
31 * Random_Gaus()*std + ave;
32 *
33 * It uses the Marsaglia polar method:
34 * http://en.wikipedia.org/wiki/Marsaglia_polar_method
35 * Based on the "MIT like" licensed code
36 * ftp://ftp.taygeta.com/pub/c/boxmuller.c:
37 * "(c) Copyright 1994, Everett F. Carter Jr.
38 * Permission is granted by the author to use
39 * this software for any application provided this
40 * copyright notice is preserved."
41 *
42 * but with div by 0 protection
43 */
bs_random_Gaus()44 double bs_random_Gaus(){
45 double x,w;
46 static double y;
47 static int use_last = 0;
48 double toreturn;
49
50 if (use_last == 0) {
51 do {
52 x = 2*((double)random() / RAND_MAX)-1;
53 y = 2*((double)random() / RAND_MAX)-1;
54
55 w = x * x + y * y;
56 } while (w >= 1 || w == 0);
57
58 w = sqrt(-2.0 * log(w) / w);
59 toreturn = x * w;
60 y = y * w;
61 } else {
62 toreturn = y;
63 }
64
65 use_last = 1 - use_last;
66
67 return toreturn;
68 }
69
70 /**
71 * Like Random_Gaus but return a complex double instead
72 * Returns a sample of a N(0,1) + 1i*N(0,1) with uncorrelated I and Q parts
73 */
bs_random_Gaus_c()74 double complex bs_random_Gaus_c(){
75 double x, y, w;
76 double complex toreturn;
77
78 do {
79 x = 2*((double)random() / RAND_MAX)-1;
80 y = 2*((double)random() / RAND_MAX)-1;
81
82 w = x * x + y * y;
83 } while(w >= 1 || w == 0);
84 w = sqrt(-2.0 * log(w) / w);
85
86 toreturn = x * w + ( y * w )*I;
87
88 return toreturn;
89 }
90
91 /**
92 * Fill in a buffer of size <size> with complex random gausian numbers
93 * each a sample of a N(0,1) + 1i*N(0,1) with uncorrelated I and Q parts
94 */
bs_random_Gaus_c_buffer(double complex * buffer,uint size)95 void bs_random_Gaus_c_buffer(double complex* buffer, uint size){
96 double x, y, w;
97 uint i;
98
99 //TOOPT: this takes a sizeable amount of time for the 2G4_channel_Indoorv1,
100 // it could make sense to optimize it to avoid using log
101 // roughly time is used as follows: log: 383/856, sqrt = 12/856, random 133*2/856
102 // an option would be: http://en.wikipedia.org/wiki/Ziggurat_algorithm
103
104 for (i = 0 ; i < size; i++){
105 do {
106 x = 2*((double)random() / RAND_MAX)-1;
107 y = 2*((double)random() / RAND_MAX)-1;
108
109 w = x * x + y * y;
110 } while(w >= 1 || w == 0);
111 w = sqrt(-2.0 * log(w) / w);
112
113 buffer[i] = x * w + ( y * w )*I;
114 }
115 }
116
117 /**
118 * Return a Generalized Pareto random number (double)
119 * Where the generalized pareto parameters are (like in MATLAB)
120 * k: shape (-inf..inf) (xi for some authors)
121 * sigma: scale (0..inf)
122 * theta: location (-inf..inf) (mu for some authors)
123 */
bs_random_GPRND(double k,double sigma,double theta)124 double bs_random_GPRND(double k, double sigma, double theta){
125 double u = ( ((double)random()) / ((double)RAND_MAX)); //U(0,1) : uniformly distributed random number from 0 to 1
126 double r;
127 if ( fabs(k) < 0.005 ){ //effectively there is no difference in the shape of the distribution from k = 0
128 r = -log(u);
129 } else {
130 r = expm1(-k*log(u)) / k; // ( u^(-k) - 1 ) / k
131 }
132 r = theta + sigma*r;
133 return r;
134 }
135
136 /**
137 * Draws a bit {0,1}
138 */
bs_random_bit()139 char bs_random_bit(){
140 return ( random() & 1 );
141 }
142
143 /**
144 * Draws a number from a U[0,1]
145 */
bs_random_uniform()146 double bs_random_uniform(){
147 return ( ((double)random()) / ((double)RAND_MAX));
148 }
149
150 /**
151 * Draws a number from a U[a,b] (both double)
152 */
bs_random_uniformR(double a,double b)153 double bs_random_uniformR(double a, double b){
154 double u = ( ((double)random()) / ((double)RAND_MAX));
155 u *=(b-a);
156 u +=a;
157 return u;
158 }
159
160 /**
161 * Draws a number from a U[a,b] (both integer)
162 * b > a
163 */
bs_random_uniformRi(int a,int b)164 int bs_random_uniformRi(int a, int b){
165 if ( b <= a )
166 return a; //safety feature
167
168 int i_u;
169 while (1) {
170 double u = ( ((double)random()) / ((double)RAND_MAX)); //[0..1]
171 u *=(b+1-a); //[0.. (b+1-a)]
172 u +=a; //[a..b+1]
173 i_u = floor(u); //[a,..,b+ups]
174 if (i_u <= b) //although extremely unlikely, we could get RAND_MAX out of the random generator and otherwise produce b+1
175 break;
176 }
177 return i_u;
178 }
179
bs_random_uint32()180 uint32_t bs_random_uint32(){
181 return (uint32_t)random();
182 }
183
