计算机噪声

序

　　大家多少都遇到过关于“随机数”和“噪声”的问题，随机数和噪声直接关系到仿真结果的好坏。而计算机产生的随机数是伪随机数，不真正的随机，那么我们使用的随机数，包括Matlab产生的随机数的质量到底如何呢？什么时候使用什么样的随机数可以满足系统的需要呢？下面我就我搜集到的资料对计算机噪声做一点介绍，希望能对大家有帮助。

　　如果你是初学者或者想对计算机噪声有一个全面的梳理的话，请从第一部分开始；如果你需要了解计算机噪声的产生方法分类以及各种测试标准那么请从第二部分开始；如果你只需要了解什么时候该用什么样的算法，可以从第三部分开始，这里包含了常用算法的比较；当然如果你对计算机噪声有很好的了解的话，附件里有一些测试过的代码资源，或许对你有帮助。当然阅读这篇文章的时候最好有网络，并且能动手去搜索一些东西的话，你的收获会更多。

 

 

【声明】：十分抱歉这个得放在前面，因为我引用了一些开源的库，^_^。

　　　　1.资料均为网上搜集，代码均为本人验证，测试平台：VS2008 SP1。

　　　　2.部分代码来自开源网站，因此附件中的代码遵守GNU规范中的LPL规范。

 

一、概念

 白噪声（White Noise）——参考维基的解释

　　白噪声的白是对于频域来说的。理想的白噪声具有无限带宽，因而其能量是无限大，在数学处理上比较方便，是系统分析的有力工具。实际应用中，一般只要一个噪声过程所具有的频谱宽度远远大于它所作用系统的带宽，并且在该带宽中其频谱密度基本上可以作为常数来考虑，就可以把它作为白噪声来处理。例如，热噪声和散弹噪声在很宽的频率范围内具有均匀的功率谱密度，通常可以认为它们是白噪声。白噪声的概率密度函数（样本统计直方图）均匀分布，因此又称作均匀白噪声或随机噪声。

 

高斯噪声（Gaussian Noise）——参考维基的解释

　　高斯噪声的所谓高斯是对于幅度的概率密度分布来说的。

 

加性高斯白噪声（Additive White Gaussian Noise，AWGN）——参考维基的解释

　　高斯白噪声中的高斯是指概率分布是正态函数，而白噪声是指它的二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。这是考查一个信号的两个不同方面的问题。

 

随机序列——随机数学中的概念，是指物理世界中随机出现的事件的不确定序列，它是随机产生，不可预测，不会重复的序列。由随机序列构成的信号称为随机噪声，即均匀白噪声。

 

伪随机序列——如果一个序列，一方面它是可以预先确定的，并且是可以重复地生产和复制的；一方面它又具有某种随机序列的随机特性（即统计特性），我们便称这种序列为伪随机序列。伪随机序列是计算机产生各种噪声的基础。

　　伪随机序列具有某种随机特性，但其序列是确定的。它们是通常由数字计算机产生确定序列，然而他们却具有某种随机序列的随机特性。因为同样具有随机特性，无法从一个已经产生的序列的特性中判断是真随机序列还是伪随机序列，只能根据序列的产生办法来判断。伪随机序列系列具有良好的随机性和接近于白噪声的相关函数（δ函数），并且有预先的可确定性和可重复性。这些特性使得伪随机序列得到了广泛的应用，特别是在CDMA系统中作为扩频码已成为CDMA技术中的关键问题。

 

常用的伪随机序列——m序列、Gold 序列、Walsh序列等。 

 

（伪）随机数产生器（RNG）——是指产生（伪）随机序列（通常是整型）的算法或方法。

 

二、噪声产生方法概述

  

关于伪随机序列（白噪声）的产生——目前有以下几种常用的算法：

● 线性同余算法（Matlab 4.0及之前版本rand函数采用）

● Shift-Register 方法

● Lagged-Fibonacci算法（Matlab 5.0及之后版本rand函数采用）

● Monte Carlo 方法

● Mid-square method, 1940s

● Mersenne Twister, 1988

● Combined Generators (Combine the outputs of two or more RNGs, eg, using)

  ◆ The universal generators

     Combine 2 generators: Lagged Fibonacci, Xn+1 = (Xn - k) mod 1677213

  ◆ The KISS generator (Keep It Simple, Stupid)

     Combine three simple generators: congruential generator, 3-shift generator, Multiply-with-carry generator

● Blum-Blum-Shub (BSS) generators, 1986

● The HAVEGE generator

 

伪随机序列衡量标准 

● Fast, especially in simulation 快

● Well distributed (pass all statistical tests known)

● Independent 独立的

● Portable and reproducible 在不同的计算机能重复产生 (for verifying simulation results)

● Long periods (for deterministic RNGs) 长周期

● Unpredictable and irreproducible (for cryptography)

● Security 保密 (for cryptography)

● Large seed spaces (for deterministic RNGs) 种子的选择要够多 

 

伪随机序列的质量检验

 

拟合优度测试（Goodness-of-fit test），用于测量产生的样本分布与期望分布之间的差异。

有以下几种常用的测试方法：

● Pearson's chi-square test

● The Kolmogorov-Smirov (KS) test （最常用）

● The Anderson and Darling (AD) test

 

统计检验—— The collision test （碰撞检验）

● Knuth's collection test

● NIST 建议的16种用于检验密码安全的统计测试

  ◆ Frequency (Monobit) Test

  ◆ Frequency Test within a Block

  ◆ Runs Test

  ◆ Tests for the longest Run of Ones in a Block

  ◆ Binary Matrix Rank Test

  ◆ Discrete Fourier Transform (Spectral) Test

  ◆ Non-overlapping Template Matching Test

  ◆ Overlapping Template Matching Test

  ◆ Maurer’s “Universal Statistical” Test

  ◆ Lampel-Ziv Compression Test

  ◆ Linear Complexity Test

  ◆ Serial Test

  ◆ Approximate Entropy Test

  ◆ Cumulative Sums (Cusum) Test

  ◆ Random Excursions Test

  ◆ Random Excursions Variant Test

● Diehard RNGs examining package （最常用）

  ◆ Birthday Spacings （三项必测之一，通过此三项测试即可通过其它所有测试）

  ◆ GCD （三项必测之一）

  ◆ Gorilla 大猩猩 （三项必测之一）

  ◆ Overlapping Permutations

  ◆ Binary Rank nn

  ◆ Binary Rank 68

  ◆ Monkey Tests OPSO, OQSO, DNA

  ◆ Count the 1’s

  ◆ Count the 1’s specific

  ◆ Parking Lot

  ◆ Minimum Distance

  ◆ Random Spheres

  ◆ The Squeeze

  ◆ Overlapping Sums

  ◆ Runs Up and Down

  ◆ The Craps

 

将伪随机数序列转为其它分布

 

一个RNG的输出是均匀分布（如[0, 2^32-1]）的随机整数。但在实际应用中，我们经常使用其它分布的随机数，比如：

● Uniform in [0, 1) （均匀分布）

● Normal (正态分布) （即高斯噪声）

● Exponential (指数分布)

● Gamma (伽玛分布)

● Poisson (泊松分布)

● Binomial (二项式分布)

 

这就需要转换。常用的转换方法有：

● Density distribution function (分布密度函数法)

● The acceptance-rejection method (接受-拒收法)

● The Monty Python Method  (拼凑法)

● The Ziggurat method （速度快）

● The alias method for generating discrete variates

● The straightforward table look up method (查表法)

 

三、噪声产生方案

　　这里以产生最常用的正态分布的伪随机序列（高斯白噪声）为主线，介绍涉及的各种操作和知识。高斯白噪声的产生通常分为两步：先产生均匀分布的白噪声，然后通过均匀分布的白噪声获得高斯白噪声。

 

　　【注】从这里开始我将开始引用英文版的维基百科的解释的链接，因为我注意到维基的英文版是解释最全面而且更新比较及时的，也因为维基的解释已经比较全面，所以我只针对对我们的有用的部分做必要的强调和针对性的说明。

 

1.均匀分布的白噪声的产生

　　古老的线性同余序列LCG(linear congruential generator)代表了最好最朴素的伪随机数产生器算法。参考维基的LCG解释。公式中参数c,m,a直接影响了产生的伪随机数的质量。一般而言，高LCG的m是2的指数次幂(一般2^32或者2^64)，因为这样取模操作截断最右的32或64位就可以了。多数编译器的库中使用了该理论实现其伪随机数发生器rand()。部分编译器使用的各个参数值在维基“Parameters in common use”这一段有描述，其中初始值一般可以通过代码用当时的时间来产生（例如C/C++用srand()函数）。rand()函数产生的随机数速度很慢，关键是随机性非常差，几乎不能胜任一般的应用。也就是如果你用C/C++然后使用srand()+rand()的模式产生随机数的话，效果是很差的，首先他的周期不够长，其次均匀性也不好。

　　LCG不能用于随机数要求高的场合，例如不能用于Monte Carlo模拟，不能用于加密应用。
　　LCG有一些严重的缺陷，例如如果LCG用做N维空间的点坐标，这些点最多位于m1/n超平面上(Marsaglia定理)，这是由于产生的相继X(n)值的关联所致。另外一个问题就是如果m设置为2的指数，产生的低位序列周期远远小于整体。
　　一般而言，输出序列的基数b中最低n位，bk = m (k是某个整数)，最大周期bn。但是有些场合LCG有很好的应用，例如内存很紧张的嵌入式中，电子游戏控制台用的小整数，使用高位可以胜任。

 

　　如果需要高质量的伪随机数，内存充足(约2kb)，Mersenne twister算法是个不错的选择。参考维基中MT的解释。Mersenne twister产生随机数的质量几乎超过任何LCG。不过一般Mersenne twister的实现使用LCG产生种子。
　　Mersenne twister是Makoto Matsumoto (松本)和Takuji Nishimura (西村)于1997年开发的伪随机数产生器，基于有限二进制字段上的矩阵线性再生。可以快速产生高质量的伪随机数，修正了古老随机数产生算法的很多缺陷。 Mersenne twister这个名字来自周期长度通常取Mersenne质数这样一个事实。常见的有两个变种Mersenne Twister MT19937和Mersenne Twister MT19937-64。

　　Mersenne Twister有很多长处，例如：周期2^19937 - 1对于一般的应用来说，足够大了，序列关联比较小，能通过很多随机性测试。用Mersenne twister算法实现的伪随机数版本非常多。例如boost库中的高质量快速随机数产生器就是用Mersenne twister算法原理编写的。

　　在新的C++标准中建议如果不知道如何选择随机数的时候，使用mt19937作为随机数产生器。它是快速的，并有可接受的质量。

　　当然还有Matlab使用的lagged_fibonacci607算法，算法描述参考维基中LF的解释。值得注意的是，它是直接产生均匀分布的浮点数的，而且比较快速。下面是我搜集到的一些RNG的参数：

 

 发生器                循环长度       近似内存需求    近似相对速度   注释
minstd_rand               2^31-2          sizeof(int32_t)           40             -
rand48                      2^48-1          sizeof(uint64_t)         80               -
lrand48（C库）         2^48-1                   -                     20          全局状态
ecuyer1988                近似2^61       2*sizeof(int32_t)        20              -
kreutzer1986                 ?             1368*sizeof(uint32_t)   60            -
hellekalek1995            2^31-1         sizeof(uint32_t)             3          良好的均匀分布in several dimensions
mt11213b                   2^11213-1     352*sizeof(uint32_t)   100         良好的均匀分布in up to 350 dimensions
mt19937                    2^19937-1     625*sizeof(uint32_t)   100         良好的均匀分布in up to 623 dimensions
lagged_fibonacci607    近似2^32000   607*sizeof(double)      150           -
lagged_fibonacci1279  近似2^67000  1279*sizeof(double)    150             -
lagged_fibonacci2281  近似2^120000 2281*sizeof(double)    150             -
lagged_fibonacci3217  近似2^170000 3217*sizeof(double)     150            -
lagged_fibonacci4423  近似2^230000 4423*sizeof(double)    150             -
lagged_fibonacci9689  近似2^510000  9689*sizeof(double)    150            -

 

【资源】

 

　　均匀分布的代码资源如图，一个工程包含了几个项目：（需要的人可以留下邮箱或者给我发邮件索要。我的邮箱：gloomy2110@163.com）

　　◆关于MT算法我找到了一个略加修改的C/C++版本，文件格式是*.hpp，主要内容类似于Boost库中Random类中的东西，可以供大家参考。

　　◆Feigenbaum算法是在网上遇到的一个算法，他的LCG生成是传统LCG加上混洗技术完成的，然后将生成的随机数作为著名的混沌迭代——费根包姆迭代结合起来，最终生成“真正意义上”的随机数，本人只对程序进行了验证，保证可以运行，初步试验认为随机数质量不错。

　　◆微软发布的Visual Studio 2008的Services Pack 1 包含了它之前发布的feature pack以及完整的TR1支持，关于TR1中的random的支持可以参考MSDN。TR1支持是对C++标准修订的一个响应，所以包含了C++标准中新添加进来的Boost库（一个很好的C++标准库，查看Boost的主页）中的一些随机数生成算法。因此，使用VS2008SP1及以上版本的IDE可以直接调用库函数生成符合要求的随机数。附件里我给出了一个例子，具体的还有很多新的东西大家可以在需要的时候好好挖掘MSDN。

　　

2.高斯白噪声的产生，目前有以下几种常用的方法：

 ● 直接按照分布生成的方法

　　这种方法是按照Gauss分布的正态分布公式直接计算出需要的点数的Gauss序列，直接用公式计算即可，只是这样的速度和精度并不理想，而且如果需要的点序列较长的情况下，会占用大量的内存，见到过莫名的栈溢出现象，一般只推荐短长度序列噪声的产生。

● Box-Muller 方法（Matlab 4.0及之前版本采用 <randn函数>）

　　最常见的算法是Box-Muller算法，参考维基的解释，以及一个变种。这个方法是个更一般的方法，可以一次产生两个正交分布的Gauss序列，速度比较快，但是理论上有当序列长度为无穷大时，分布才是真正Gauss的，工程上一般可以采用此方法。

● Ziggurat 方法（Matlab 5.0及之后版本采用 <randn函数>）

 　　最新的Matlab使用的是Ziggurat算法，算法描述参考维基的解释。由于算法没有速度缓慢的浮点运算，因此算法速度非常快。

●ratio-method

　　这个算法是GSL（又是一个很好的开源库，查看GSL的主页）中自带的一个算法，具体算法性能和随机数的质量有待于考究，不过鉴于GSL的性质，我认为这个算法和Box-Müller算法在一个等级上。

　　下面的表格是我搜集到的有人做的后面三种算法的在一秒的时间内生成的随机数的数量比较（单位：百万），RNG一栏表示的是GSL中不同的随机数生成器。

 

RNG	Ziggurat	Box-Müller	ratio method	improvement
borosh13	8.89	2.31	2.39	3.7
cmrg	2.84	1.34	1.3	2.1
coveyou	8.88	2.39	2.07	3.7
fishman18	1.23	0.882	0.863	1.4
fishman20	5.88	2.21	2.13	2.7
fishman2x	4.35	1.69	1.51	2.6
gfsr4	9.03	2.19	2.29	4.0
knuthran	6.47	2.21	2.27	2.8
knuthran2	0.662	0.527	0.499	1.3
lecuyer21	6.61	2.22	2.17	3.0
minstd	5.61	2.16	2.1	2.6
mrg	4.38	1.73	1.66	2.5
mt19937	8.28	2.01	2.14	3.9
mt19937_1999	8.28	2.02	2.14	3.9
mt19937_1998	8.27	2.01	2.12	3.9
r250	9.1	2.2	2.05	4.1
ran0	5.61	2.17	2.1	2.6
ran1	5.94	2.1	2.04	2.8
ran2	4.4	1.68	1.63	2.6
ran3	6.69	2.05	2.12	3.2
rand	8.47	2.37	2.37	3.6
rand48	8.14	0.887	0.84	9.2
random128-bsd	7.91	2.2	2.3	3.4
random128-glibc2	7.72	2.18	2.31	3.3
random128-libc5	7.94	2.2	2.31	3.4
random256-bsd	8.02	2.18	2.29	3.5
random256-glibc2	8	2.2	2.28	3.5
random256-libc5	7.72	2.18	2.29	3.4
random32-bsd	8.14	2.24	2.34	3.5
random32-glibc2	7.93	2.25	2.34	3.4
random32-libc5	8.19	2.24	2.34	3.5
random64-bsd	7.71	2.2	2.32	3.3
random64-glibc2	7.71	2.21	2.33	3.3
random64-libc5	7.68	2.2	2.32	3.3
random8-bsd	8.45	2.37	2.36	3.6
random8-glibc2	8.45	2.37	2.36	3.6
random8-libc5	8.47	2.37	2.37	3.6

 

*All measurements were done on a Pentium 4 with 2 GHz, the code was compiled with gcc -O2. The last column gives the ratio of the values for the Ziggurat method and for the better one of the other two methods.

 

【资源】

 

 

　　关于Gauss分布的噪声的生成资源是一个工程，包含了几个项目：

　　（需要的人可以留下邮箱或者给我发邮件索要。我的邮箱：gloomy2110@163.com）

　　◆ZigguratTab并不是用来生成噪声的，而是用来生成一些代码（主要是表格），因为Ziggurat算法有查表操作，因此这个程序是专门用来生成表格的，当然为了方便起见，代码里已经为我们完成了变量的申请等等，只是他用的是printf，我们需要把输出重定向到文件里，当然文件的后缀是*.c或者*.cpp，这个是真正意义上的用代码写代码了，呵呵～

◆ZigGauss是从GSL库里面分离出来单独形成的一个工程，工程性质为测试Ziggurat算法生成的随机数的质量，算法实现主要参考ZigGauss.c即可。（注意这里的LCG是使用mt19937产生的，具体实现参考mt.c，它也包含了mt19937以及两个变种的实现，^_^。）

　　◆TimeGauss也是从GSL库里面分离出来单独形成的一个工程，工程性质为测试各种算法的速度，上面的表格可以用这个工程验证，不过要注意只能验证基于我移植的MT系列LCG算法的Gauss随机数生成速度，BoxGauss.c为Box-Muller算法实现，RatGauss.c为ratio-method算法实现，ZigGauss.c为Ziggurat算法实现。（注：这里我对原来的timeGauss.c进行了改动，由于GSL是Linux系统下的东西，移植的时候我尽量保持原代码不变，但是在这个文件中，做了比较大的改动，替换了原来的时间处理，当然我采用的是Windows提供的7种计时方式中计时精度最高的一种，虽然这个精度和你的CPU的时钟频率有关，不过幸运的是最差的误差也不会超过1us）

　　◆tr1Gauss是利用TR1的库生成Gauss分布的噪声的工程，也是其中最简单的一个，TR1还提供了很多其他分布的噪声的生成模板，但是具体是什么算法我并没有找到，具体的用法请大家参考MSDN尽力挖掘吧。