它是否存在具有平坦功率谱的随机发生器？ [英] Does it exist a random generator with flat power spectrum?
本文介绍了它是否存在具有平坦功率谱的随机发生器？的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我对一些rand（）函数性能进行了以下测试。随机0,2,3使用c函数rand（）
我通过制作循环表来设计random4（）
我通过制作数据阵列的fft，通过计算功率谱来进行测试以寻找最佳的测试。我还计算了每秒的随机微积分操作。
结果通过使用fft我可以看到所有方法都相似但是灾难因为rms误差大约为2.9-3而平均值= 5.5所以频谱的平均误差为54％平均误差
另一个结果算法random0,2,3每秒计算45-46百万个随机数，random4计算每秒256兆随机数（在发布模式下使用3Ghz处理器，在调试模式下仅20兆数/秒）
random3避免使用rand（）时它达到RAND_MAX
I made following tests of some rand() functions performance. The random0,2,3 uses the c function rand()
I designed the random4() by making a loop up table
I performed a test looking for the best one by calculation of the the power spectrum by making fft of array of data. I calculated also the random calculus operations per second.
As result by using the fft I could see that all the methods are similar but a disaster because rms error is about 2.9-3 and the average value=5.5 so spectrum has an average error of 54% from the flat one
Another result algorithms random0,2,3 calculates 45-46 million random numbers per second and the random4 calculates 256 mega random numbers per second (in release mode by using a 3Ghz processor, in debug mode only 20 mega numbers/second)
The random3 avoids use rand() when it reach RAND_MAX
#include <stdio.h>
#include <stdlib.h>//rand() used
#include <limits.h>//used by random0()
#include <math.h> //to make fft test
#include <time.h>//to make timer test and seeds the randoms
#define SWAP(t,a,b) ( (t) = (a), (a) = (b), (b) = (t) ) //use: SWAP(dtemp, d1, d2);
//Definitions, constants and functions used by random4 test:
#define RANS_MAX 131072
double rans_m[RANS_MAX];
double rans();//generates random number 0 <= rans() < 1
void srans(int seed); //Must be run one time before using rans()
//Four random functions:
long random0(long max);//good
long random2(long max) { return rand() % max; }//better than random0 22% and quicker
long random3(long max) { long r; do{r=rand();}while (r==RAND_MAX);return (long) ((r*max)/RAND_MAX); }//It is avoided to use rand() when it reach RAND_MAX
long random4(long max) { return (long) (max*rans());}//It is the best of all but needs 177ms to execute srans() before to fill rams_m[] array
//Constants and functions used for testing
#define NUM 10000000L
typedef long (*FUNCT)(long max);
const long MAX=128;//rand values maximum=MAX-1
long sum[MAX+1];
//tests the random functions by generation of numbers between 0-19 NUM times
void test(FUNCT f_rand);
//Constants and functions used to calculate power spectrum
#define PI 3.1415926535897932384626433832795
#define TWOPI 6.283185307179586476925286766559
#define PI2 6.283185307179586476925286766559//input: data_in[0..number_of_samples]
//output:power[0..NFFT/2] and frecuencia[0..NFFT/2]. CAUTION!!: fix power[] and frecuencia[] dimension to 2 x number_of_samples to avoid stack overflow!!!!
//NFFT=real size of power[] and frecuencia[] that is greater than number_of_samples in the form 2^N
//sample_rate=samples per second of the acquisition system
void power_spectrum(double data_in[], long number_of_samples,long sample_rate,double power[],double frecuencia[],long &NFFT);
long main()
{
printf("\nTimings will be reduced by compile in RELEASE instead DEBUG mode");
clock_t ini,fin;ini=clock();
srans(time (NULL));//initializes rands_m[] with seed=1000;
fin=clock();double milliseconds=1000.0*(fin-ini)/CLOCKS_PER_SEC;
srand(time (NULL));//initializes rand()
printf("\nTime needed for srans()=%lg milliseconds",milliseconds);
printf("\n\nTesting random0()");test(random0);
printf("\n\nTesting random2()");test(random2);
printf("\n\nTesting random3()");test(random3);
printf("\n\nTesting random4()");test(random4);
printf("\n\n===FIN===");getchar();getchar();
return 1;
}
void test(FUNCT f_rand)
{
long i,max=MAX;
for (i=0;i<MAX+1;i++) sum[i]=0;
clock_t ini,fin;ini=clock();
for (i=0;i<NUM;i++)
{
sum[(*f_rand)(max)]++;
}
fin=clock();double seconds=1.0*(fin-ini)/CLOCKS_PER_SEC;
long min=LONG_MAX,max1=-1;
for (i=0;i<MAX;i++)
{
if (min>sum[i]) min=sum[i];
if (max1<sum[i]) max1=sum[i];
}
printf("\nMegasamples/second=%lg",1.0e-6*NUM/seconds);
printf("\nMinimum repeatings=%li",min);
printf("\nMaximum repeatings=%li",max1);
printf("\nMaximum-minimums =%li (less is better)",max1-min);
if (sum[MAX]!=0) printf("\nFAILS maximum value reached one time every %lg", (double) (1.0*NUM/sum[MAX]));
double power[MAX],data_in[MAX],frecuencia[MAX];
long NFFT;
for (i=0;i<MAX;i++) data_in[i]=(*f_rand)(max);
power_spectrum(data_in, MAX, MAX,power,frecuencia,NFFT);
double mid=0.0,rms=0.0,x;
for (i=0;i<MAX/2;i++) mid+=power[i];
mid=2.0*mid/MAX;
for (i=0;i<MAX/2;i++) { x=mid-power[i];rms+=x*x;}
rms=sqrt(2.0*rms/MAX);
printf("\nRMS of FFT power =%lg (less is better)",rms);
}
long random0(long max)
{
static long sum=0;
long x;
x=(rand()+sum) % max;
sum+=RAND_MAX;
if (sum>(LONG_MAX-RAND_MAX))
sum=0;
return x;
}
void srans(int seed=0)
{
long i,k,l=0;double temp;
for (i=0;i<RANS_MAX;i++) rans_m[i]=1.0*i/RANS_MAX;
for (i=0;i<(8L*1024L*1024L);i++)
{
k=i % RANS_MAX;l=(k+rand()) % RANS_MAX;
SWAP(temp,rans_m[k],rans_m[l]);
//temp=rans_m[k];rans_m[k]=rans_m[l];rans_m[l]=temp;
}
for (i=0;i<(seed % RANS_MAX);i++) //seed is applied, but not longer than 128k times
rans();
}
double rans()
{
static int index=-1;
index=(index+1)% RANS_MAX; //index must be < RANS_MAX
return rans_m[index];
}
//input: data_in[0..number_of_samples]
//output:power[0..NFFT/2] and frecuencia[0..NFFT/2]. CAUTION!!: fix power[] and frecuencia[] dimension to 2 x number_of_samples to avoid stack overflow!!!!
//NFFT=real size of power[] and frecuencia[] that is greater than number_of_samples in the form 2^N
//sample_rate=samples per second of the acquisition system
void power_spectrum(double data_in[], long number_of_samples,long sample_rate,double power[],double frecuencia[],long &NFFT)
{
long i;
double *data;
double *X;
NFFT = (long) pow(2.0, ceil(log((double)number_of_samples)/log(2.0)));
/* allocate memory for NFFT complex numbers (note the +1) */
X = (double *) malloc((2*NFFT+1) * sizeof(double));
/* Storing x(n) in a complex array to make it work with four1.
This is needed even though x(n) is purely real in this case. */
for(i=0; i<number_of_samples; i++)
{
X[2*i+1] = data_in[i];
X[2*i+2] = 0.0;
}
/* pad the remainder of the array with zeros (0 + 0 j) */
for(i=number_of_samples; i<NFFT; i++)
{
X[2*i+1] = 0.0;
X[2*i+2] = 0.0;
}
data=X;
long n, mmax, m, j, istep;
double wtemp, wr, wpr, wpi, wi, theta;
double tempr, tempi;
n = NFFT << 1;
j = 1;
for (i = 1; i < n; i += 2)
{
if (j > i)
{
tempr = data[j]; data[j] = data[i]; data[i] = tempr;
tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;
}
m = n >> 1;
while (m >= 2 && j > m)
{
j -= m;
m >>= 1;
}
j += m;
}
mmax = 2;
while (n > mmax)
{
istep = 2*mmax;
theta = TWOPI/(1.0*mmax);
wtemp = sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for (m = 1; m < mmax; m += 2)
{
for (i = m; i <= n; i += istep)
{
j =i + mmax;
tempr = wr*data[j] - wi*data[j+1];
tempi = wr*data[j+1] + wi*data[j];
data[j] = data[i] - tempr;
data[j+1] = data[i+1] - tempi;
data[i] += tempr;
data[i+1] += tempi;
}
wr = (wtemp = wr)*wpr - wi*wpi + wr;
wi = wi*wpr + wtemp*wpi + wi;
}
mmax = istep;
}
//I have now in data[] the real and imaginary parts of the FFT
//I calculate now Power and frequency:
double kk=0.5*number_of_samples;//It must be 0.5*NFFT, but we introducted samples only in the 0-number_of_samples segments, others were filled with zeroesd
for(i=1; i<=NFFT/2; i++)
{
frecuencia[i-1] =1.0*(i-1)*sample_rate/NFFT;
power[i-1] =data[2*i+1]*data[2*i+1]+ data[2*i+2]*data[2*i+2];
power[i-1] =sqrt(power[i-1])/kk;
}
free(X);
}
推荐答案
我测试了drand48（参见函数random5）和std :: mt19937 gen（seed（）） （参见函数random6）但不幸的是，性能类似于之前使用的random0和random3：
I have tested also drand48 (see function random5) and std::mt19937 gen(seed()) (see function random6) but unfortunately the performance is similar to the random0 and random3 used before:
#include <stdio.h>
#include <random> //Used in random6
#include <stdlib.h>//rand() used
#include <limits.h>//used by random0()
#include <math.h> //to make fft test
#include <time.h>//to make timer test and seeds the randoms
#define SWAP(t,a,b) ( (t) = (a), (a) = (b), (b) = (t) ) //use: SWAP(dtemp, d1, d2);
const long MAX = 128;//rand values maximum=MAX-1
//Definitions, constants and functions used by random4 test:
#define RANS_MAX 131072
double rans_m[RANS_MAX];
double rans();//generates random number 0 <= rans() < 1
void srans(long seed); //Must be run one time before using rans()
//Definitions for drand48 and random5
#define C_DRAND 16807
#define A_DRAND 2147483647.0
double yz_DRAND = 0.1;
void srand48(long seed) { yz_DRAND = (double)seed; }
double drand48();
//Definitions for random6 (C++11 only)
std::random_device seed;
std::mt19937 gen(seed());
std::uniform_int_distribution<int> dist(0, MAX-1);//Ojo, si no falla
//Four random functions: (worst=random3 and 5)
long random0(long max);//good
long random2(long max) { return rand() % max; }//better than random0 22% and quicker
long random3(long max) { long r; do { r = rand(); } while (r == RAND_MAX); return (long)((r*max) / RAND_MAX); }//Bad! It is avoided to use rand() when it reach RAND_MAX
long random4(long max) { return (long)(max*rans()); }//It is the best of all but needs 177ms to execute srans() before to fill rams_m[] array
long random5(long max) { return (long)(max*drand48()); }//As bad as random3
long random6(long max) { return dist(gen); }
//Constants and functions used for testing
#define NUM 10000000L
typedef long(*FUNCT)(long max);
long sum[MAX + 1];
//tests the random functions by generation of numbers between 0-19 NUM times
void test(FUNCT f_rand);
//Constants and functions used to calculate power spectrum
#define PI 3.1415926535897932384626433832795
#define TWOPI 6.283185307179586476925286766559
#define PI2 6.283185307179586476925286766559//input: data_in[0..number_of_samples]
//output:power[0..NFFT/2] and frecuencia[0..NFFT/2]. CAUTION!!: fix power[] and frecuencia[] dimension to 2 x number_of_samples to avoid stack overflow!!!!
//NFFT=real size of power[] and frecuencia[] that is greater than number_of_samples in the form 2^N
//sample_rate=samples per second of the acquisition system
void power_spectrum(double data_in[], long number_of_samples, long sample_rate, double power[], double frecuencia[], long &NFFT);
long main()
{
printf("\nTimings will be reduced by compile in RELEASE instead DEBUG mode");x
clock_t ini, fin; ini = clock();
srans(time(NULL));//initializes rands_m[] with seed=1000;
fin = clock(); double milliseconds = 1000.0*(fin - ini) / CLOCKS_PER_SEC;
srand(time(NULL));//initializes rand()
srand48(time(NULL));//initializes drand48()
printf("\nTime needed for srans()=%lg milliseconds", milliseconds);
printf("\n\nTesting random0()"); test(random0);
printf("\n\nTesting random2()"); test(random2);
printf("\n\nTesting random3()"); test(random3);
printf("\n\nTesting random4()"); test(random4);
printf("\n\nTesting random5()"); test(random5);
printf("\n\nTesting random6()"); test(random6);
printf("\n\n===FIN==="); getchar(); getchar();
return 1;
}
void test(FUNCT f_rand)
{
long i, max = MAX;
for (i = 0; i<max mode="hold" />
clock_t ini, fin; ini = clock();
for (i = 0; i<num;> {
sum[(*f_rand)(max)]++;
}
fin = clock(); double seconds = 1.0*(fin - ini) / CLOCKS_PER_SEC;
long min = LONG_MAX, max1 = -1;
for (i = 0; i<max;> {
if (min>sum[i]) min = sum[i];
if (max1<sum[i])> }
printf("\nMegasamples/second=%lg", 1.0e-6*NUM / seconds);
printf("\nMinimum repeatings=%li", min);
printf("\nMaximum repeatings=%li", max1);
printf("\nMaximum-minimums =%li (less is better)", max1 - min);
if (sum[MAX] != 0) printf("\nFAILS maximum value reached one time every %lg", (double)(1.0*NUM / sum[MAX]));
double power[MAX], data_in[MAX], frecuencia[MAX];
long NFFT;
for (i = 0; i<max;> power_spectrum(data_in, MAX, MAX, power, frecuencia, NFFT);
double mid = 0.0, rms = 0.0, x;
for (i = 0; i mid = 2.0*mid / MAX;
for (i = 0; i rms = sqrt(2.0*rms / MAX);
printf("\nRMS of FFT power =%lg (less is better)", rms);
}
long random0(long max)
{
static long sum = 0;
long x;
x = (rand() + sum) % max;
sum += RAND_MAX;
if (sum>(LONG_MAX - RAND_MAX))
sum = 0;
return x;
}
void srans(long seed = 0)
{
seed = seed % (1024L * 1024L);//limited to 1Meg
long i, k, l = 0; double temp;
for (i = 0; i for (i = 0; i<(8L * 1024L * 1024L); i++)
{
k = i % RANS_MAX; l = (k + rand()) % RANS_MAX;
SWAP(temp, rans_m[k], rans_m[l]);
//temp=rans_m[k];rans_m[k]=rans_m[l];rans_m[l]=temp;
}
for (i = 0; i<(seed % RANS_MAX); i++) //seed is applied, but not longer than 128k times
rans();
}
double rans()
{
static int index = -1;
index = (index + 1) % RANS_MAX; //index must be < RANS_MAX
return rans_m[index];
}
//input: data_in[0..number_of_samples]
//output:power[0..NFFT/2] and frecuencia[0..NFFT/2]. CAUTION!!: fix power[] and frecuencia[] dimension to 2 x number_of_samples to avoid stack overflow!!!!
//NFFT=real size of power[] and frecuencia[] that is greater than number_of_samples in the form 2^N
//sample_rate=samples per second of the acquisition system
void power_spectrum(double data_in[], long number_of_samples, long sample_rate, double power[], double frecuencia[], long &NFFT)
{
long i;
double *data;
double *X;
NFFT = (long)pow(2.0, ceil(log((double)number_of_samples) / log(2.0)));
/* allocate memory for NFFT complex numbers (note the +1) */
X = (double *)malloc((2 * NFFT + 1) * sizeof(double));
/* Storing x(n) in a complex array to make it work with four1.
This is needed even though x(n) is purely real in this case. */
for (i = 0; i<number_of_samples;> {
X[2 * i + 1] = data_in[i];
X[2 * i + 2] = 0.0;
}
/* pad the remainder of the array with zeros (0 + 0 j) */
for (i = number_of_samples; i<nfft;> {
X[2 * i + 1] = 0.0;
X[2 * i + 2] = 0.0;
}
data = X;
long n, mmax, m, j, istep;
double wtemp, wr, wpr, wpi, wi, theta;
double tempr, tempi;
n = NFFT << 1;
j = 1;
for (i = 1; i < n; i += 2)
{
if (j > i)
{
tempr = data[j]; data[j] = data[i]; data[i] = tempr;
tempr = data[j + 1]; data[j + 1] = data[i + 1]; data[i + 1] = tempr;
}
m = n >> 1;
while (m >= 2 && j > m)
{
j -= m;
m >>= 1;
}
j += m;
}
mmax = 2;
while (n > mmax)
{
istep = 2 * mmax;
theta = TWOPI / (1.0*mmax);
wtemp = sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for (m = 1; m < mmax; m += 2)
{
for (i = m; i <= n; i += istep)
{
j = i + mmax;
tempr = wr*data[j] - wi*data[j + 1];
tempi = wr*data[j + 1] + wi*data[j];
data[j] = data[i] - tempr;
data[j + 1] = data[i + 1] - tempi;
data[i] += tempr;
data[i + 1] += tempi;
}
wr = (wtemp = wr)*wpr - wi*wpi + wr;
wi = wi*wpr + wtemp*wpi + wi;
}
mmax = istep;
}
//I have now in data[] the real and imaginary parts of the FFT
//I calculate now Power and frequency:
double kk = 0.5*number_of_samples;//It must be 0.5*NFFT, but we introducted samples only in the 0-number_of_samples segments, others were filled with zeroesd
for (i = 1; i <= NFFT / 2; i++)
{
frecuencia[i - 1] = 1.0*(i - 1)*sample_rate / NFFT;
power[i - 1] = data[2 * i + 1] * data[2 * i + 1] + data[2 * i + 2] * data[2 * i + 2];
power[i - 1] = sqrt(power[i - 1]) / kk;
}
free(X);
}
double drand48()
{
long ki;
double uu;
ki = (long)((C_DRAND* yz_DRAND) / A_DRAND);
yz_DRAND = C_DRAND * yz_DRAND - ki * A_DRAND;
uu = yz_DRAND / (A_DRAND - 1);
return uu;
}
</int></time.h></math.h></limits.h></stdlib.h></random></stdio.h>这篇关于它是否存在具有平坦功率谱的随机发生器？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
