漫谈正态分布的生成

莱利分布（Rayleigh distribution）：又称为瑞利分布，当一个随机二维向量的两个分量呈独立的、有着相同的方差的正态分布时，这个向量的模呈莱利分布。例如，当随机复数的实部和虚部独立同分布于0均值，同方差的正态分布时，该复数的绝对值服从莱利分布。

莱利分布的概率密度函数是^[1]

https://upload.wikimedia.org/math/6/3/2/6325f93ac4e7389a4bce9ca8e7186ef6.png

高斯白噪声：因为对于图像来说就是2维的x-y,因此在x方向上是一个高斯分布，在Y方向上是另外一个高斯分布。需要2个高斯分布(也就是两个高斯分布)，这就是box-muller方法可以提供的2个高斯分布，并且满足白噪声的要求。

# -*- coding: utf-8 -*- 
#加性零均值高斯噪声
#code:myhaspl@myhaspl.com
import cv2
import numpy as np

fn="test2.jpg"
myimg=cv2.imread(fn)
img=cv2.cvtColor(myimg,cv2.COLOR_BGR2GRAY)

param=80
#灰阶范围
grayscale=256
w=img.shape[1]
h=img.shape[0]
newimg=np.zeros((h,w),np.uint8)



for x in xrange(0,h):
  for y in xrange(0,w,2):

 r1=np.random.random_sample()

r2=np.random.random_sample()
z1=param*np.cos(2*np.pi*r2)*np.sqrt((-2)*np.log(r1))
z2=param*np.sin(2*np.pi*r2)*np.sqrt((-2)*np.log(r1))

    fxy=int(img[x,y]+z1)
    fxy1=int(img[x,y+1]+z2)           
    #f(x,y)
    if fxy<<span class="number" style="color: rgb(0, 153, 153);">0:
      fxy_val=0
    elif fxy>grayscale-1:
      fxy_val=grayscale-1
    else:
      fxy_val=fxy
    #f(x,y+1)
    if fxy1<<span class="number" style="color: rgb(0, 153, 153);">0:
      fxy1_val=0
    elif fxy1>grayscale-1:
      fxy1_val=grayscale-1
    else:
      fxy1_val=fxy1
    newimg[x,y]=fxy_val
    newimg[x,y+1]=fxy1_val
  

cv2.imshow('preview',newimg)
cv2.waitKey()
cv2.destroyAllWindows()

http://www.bkjia.com/ASPjc/832055.html

numpy.random.normal()和numpy.random.randn()产生的随机值的分布满足高斯分布,但不是white。 计算信噪比也就是简单的对信号的乘方求和而已，例如如果你有信号x和噪声n，
那么信噪比就是：



10*log10( sum(x**2) / sum(n**2))



根据以上信息可以编写如下程序。
wgn(x, snr)中x为信号，snr为信噪比，返回满足条件的高斯白噪声，只需要：
x += wgn(x, snr)，即可以得到和matlab的awgn相同的效果。



程序中用hist()检查噪声是否是高斯分布，psd()检查功率谱密度是否为常数。
"""



import numpy as np 

import pylab as pl 

def wgn(x, snr):


 snr = 10**(snr/10.0)


 xpower = np.sum(x**2)/len(x)


 npower = xpower / snr


 return np.random.randn(len(x)) * np.sqrt(npower)

t = np.arange(0, 1000000) * 0.1


x = np.sin(t)


n = wgn(x, 6)


xn = x+n # 增加了6dBz信噪比噪声的信号


pl.subplot(211)


pl.hist(n, bins=100, normed=True)


pl.subplot(212)


pl.psd(n)


pl.show()

randn()生成的是 band-limited white noise

画2D 频谱图：

http://stackoverflow.com/questions/15541590/how-do-i-plot-fft-in-numpy

图像产生加性零均值高斯噪声，在灰度图上加上噪声，加上噪声的方式是每个点的灰度值加上一个噪声值，噪声值的产生方式为 Box-Muller算法 生成高斯白噪声。

在计算机模拟中，经常需要生成正态分布的数值。最基本的一个方法是使用标准的正态累积分布函数的反函数。除此之外还有其他更加高效的方法，Box-Muller变换就是其中之一。另一个更加快捷的方法是ziggurat算法。下面将介绍这两种方法。一个简单可行的并且容易编程的方法是：求12个在（0,1）上均匀分布的和，然后减6(12的一半)。这种方法可以用在很多应用中。这12个数的和是Irwin-Hall分布；选择一个方差12。这个随即推导的结果限制在（-6,6）之间，并且密度为12，是用11次多项式估计正态分布。

Box-Muller方法是以两组独立的随机数U和V，这两组数在(0,1]上均匀分布，用U和V生成两组独立的标准常态分布随机变量X和Y:

http://www.bkjia.com/uploads/allimg/140713/12423B424-0.png

http://www.bkjia.com/uploads/allimg/140713/12423A404-1.png。

这个方程的提出是因为二自由度的 卡方分布 很容易由指数随机变量（方程中的lnU）生成。因而通过随机变量V可以选择一个均匀环绕圆圈的角度，用指数分布选择半径然后变换成（正态分布的）x,y坐标。

Box-Muller 是产生随机数的一种方法。Box-Muller 算法隐含的原理非常深奥，但结果却是相当简单。它一般是要得到服从正态分布的随机数，基本思想是先得到服从均匀分布的随机数再将服从均匀分布的随机数转变为服从正态分布。

#include 
#include 
#include 
double generateGaussianNoise(double mu, double sigma)
{
        const double epsilon = std::numeric_limits<</span>double>::min();
        const double two_pi = 2.0*3.14159265358979323846;

        static double z0, z1;
        static bool generate;
        generate = !generate;

        if (!generate)
           return z1 * sigma + mu;

        double u1, u2;
        do
         {
           u1 = rand() * (1.0 / RAND_MAX);
           u2 = rand() * (1.0 / RAND_MAX);
         }
        while ( u1 <= epsilon );

        z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
        z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
        return z0 * sigma + mu;
}

有一个概念错误需要指出：“高斯白噪声的幅度服从高斯分布”的说法是错误的（而是概率密度函数为高斯分布)，高斯噪声的幅度服从瑞利分布。

另外，还必须区分高斯噪声和白噪声两个不同的概念。高斯噪声是指噪声的概率密度函数服从高斯分布，白噪声是指噪声的任意两个采样样本之间不相关，两者描述的角度不同。白噪声不必服从高斯分布，高斯分布的噪声不一定是白噪声。当然，实际系统中的热噪声是我们一般所说的白噪声的主要来源，它是服从高斯分布的，但一般具有有限的带宽，即常说的窄带白噪声，严格意义上它不是白噪声。

在python中有一个np.random.normal()函数可以来产生符合高斯分布的样本（只是符合高斯分布，不是white噪声,功率谱密度不为常数)：

https://groups.google.com/forum/#!topic/python-cn/kUZpnCCFTog

http://blog.csdn.net/a_31415926/article/details/50535279

http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.random.normal.html

http://sysbio.mrc-bsu.cam.ac.uk/group/images/6/6f/Infpy_gp.pdf

https://github.com/scottyhq/roipy/blob/master/noise.py

给加高斯噪声的意思，就是在原图像矩阵上面加一个符合高斯或者叫正态分布特征的矩阵。

生成随机噪声的三个方法,如果我们的目标矩阵是一个r*c的矩阵，要生成一个均值是mean，标准差sigma的随机噪声矩阵，那么是这样

sigma*np.random.randn(r,c)+mean, #输入是两个参数，一个mean，一个sigma。

sigma*np.random.standard_normal((r,c)) #输入是一个tuple，里面包含两个元素。  和上面那个是完全一样的。

numpy.random.normal(mean,sigma, size=(r,c))  #mean，sigma和大小都作为输入放进函数了。

这三个函数。np.random.randn(r,c)和np.random.standard_normal((r,c))  两个函数一样，就是输入参数一个是tuple，一个是两个参数，本质一样的，他们生成的都是标准正态分布，也就是均值是0，sigmg为1的。  如果我们要的是一般正态分布，需要自己再转。

numpy.random.normal(mean,sigma,size)这个直接生成的就是一般正态分布矩阵了。

上面生成了随机噪声矩阵，和原图叠加就好，但是这里面有个bug，不能直接加。

因为image读取出来的图像，类型是uin8的（也就是取值范围是0-255）

我们计算的随机噪声值，可能是或负或正的小数，

这样加起来的和的图像img_with_noise，其实是有可能是负值，有可能超了255的，而且这个np.ndarray类型是普通的float64了。

这个时候用imshow函数来显示，会出现全白全黑的图像bug。

解决方法是在显示前将图像做归一化normalization转成0-255区间的unint类型。 

http://dspguru.com/dsp/howtos/how-to-generate-white-gaussian-noise

How to Generate White Gaussian Noise

by Matt Donadio

White Gaussian Noise (WGN) is needed for DSP system testing or DSP system identification. Here are two methods for generating White Gaussian Noise. Both rely on having a good uniform random number generator. We will assume that the function "uniform()" returns a random variable in the range [0, 1] and has good statistical properties. (Be forewarned, though, that many standard random number generators have poor statistical properties.) Once we can generate noise signals with normal distribution, we can adjust the signal's mean and variance to use it for any purpose.

Central Limit Thorem Method

The Central Limit Theorm states that the sum of N randoms will approach normal distribution as N approaches infinity. We can outline an algorithm that uses this approach as:

   X=0

   for i = 1 to N

      U = uniform()

      X = X + U

   end

      

   

   

      

   X = X - N/2                

   X = X * sqrt(12 / N)       

When the algorithm finishes, X will be our unit normal random. X can be further modified to have a particular mean and variance, e.g.:

   X' = mean + sqrt(variance) * X

The drawback to this method is that X will be in the range [-N, N], instead of (-Infinity, Infinity) and if the calls to uniform are not truly independent, then the noise will no longer be white. Jeruchim, et. al., recommend N >=20 for good results.

Box-Muller Method

The Box-Muller method uses the technique of inverse transformation to turn two uniformly distributed randoms, U1 and U2, into two unit normal randoms, X and Y.

   do                        

      U1=uniform()   

   while U1==0

   U2=uniform()

  

   X=sqrt(-2 * log(U1)) * cos(2pi * U2)

   Y=sqrt(-2 * log(U1)) * sin(2pi * U2)

We can use some trigonometry to simplify this a bit and get:

   do

      U1=uniform()            

      U2=uniform()            

      V1=2 * U1 -1            

      V2=2 * U2 - 1           

      S=V1 * V1 + V2 * V2

   while S >=1

      

   X=sqrt(-2 * log(S) / S) * V1

   Y=sqrt(-2 * log(S) / S) * V2

The above is called the Polar Method and is fully described in the Ross book [Ros88].

X and Y will be unit normal random variables (mean = 0 and variance = 1), and can be easilly modified for different mean and variance.

   X' = mean + sqrt(variance) * X

   Y' = mean + sqrt(variance) * Y 

1. using namespace cv;  
2. using namespace std;  
3.   
4. #define TWO_PI 6.2831853071795864769252866  
5.    
6. double generateGaussianNoise()  
7. {  
8.     static bool hasSpare = false;  
9.     static double rand1, rand2;  
10.    
11.     if(hasSpare)  
12.     {  
13.         hasSpare = false;  
14.         return sqrt(rand1) * sin(rand2);  
15.     }  
16.    
17.     hasSpare = true;  
18.    
19.     rand1 = rand() / ((double) RAND_MAX);  
20.     if(rand1 < 1e-100) rand1 = 1e-100;  
21.     rand1 = -2 * log(rand1);  
22.     rand2 = (rand() / ((double) RAND_MAX)) * TWO_PI;  
23.    
24.     return sqrt(rand1) * cos(rand2);  
25. }  
26.   
27.   
28. void AddGaussianNoise(Mat& I)  
29. {  
30.     // accept only char type matrices  
31.     CV_Assert(I.depth() != sizeof(uchar));  
32.   
33.     int channels = I.channels();  
34.   
35.     int nRows = I.rows;  
36.     int nCols = I.cols * channels;  
37.   
38.     if(I.isContinuous()){  
39.         nCols *= nRows;  
40.         nRows = 1;  
41.     }  
42.   
43.     int i,j;  
44.     uchar* p;  
45.     for(i = 0; i < nRows; ++i){  
46.         p = I.ptr(i);  
47.         for(j = 0; j < nCols; ++j){  
48.             double val = p[j] + generateGaussianNoise() * 128;  
49.             if(val < 0)  
50.                 val = 0;  
51.             if(val > 255)  
52.                 val = 255;  
53.   
54.             p[j] = (uchar)val;  
55.   
56.         }  
57.     }  
58.   
59. }