实现用c++产生高斯噪声

实现用c++实现高斯噪声的产生，可以对图像进行添加，可以设定高斯噪声的均值和方差。

#include 
#include 
#include 
#include 
#include 
#include 
#include 

#define N_gauss		2000	//产生1000个点高斯序列
#define N_perpoint	100		//每个高斯点需要100个均匀分布的随机数
#define N_h			2000	//产生1000个点H（n）
#define N_ans		4000	//输出序列长度
#define PI			3.1415926f
#define GAUSS		0		//在printans时的type，使用代表gauss
#define ANS			1		//在printans时的type，使用代表ans

void init_gauss（double gauss[]）			//清gauss[N_gauss]序列全部为0
{
	int i;
	for （i=0; i	{
		gauss[i] = 0;
	}
}

void gauss_gen（double x[] double gauss[]）//产生长度为N_gauss的高斯序列
{
	int		i n stime;
	long	ltime;
	for（i=0; i	{
		ltime = time（NULL）;
		stime = （unsigned int）ltime;	
		stime = stime+i;				//随机的stime
		srand（stime）;					//获取随机种子next
		for（n=0; n			x[n] = 0;
		for（n=0; n		{
		//初始化随机数序列x[N_perpoint]为随机数范围0~1
			x[n] = （double）rand（）/（RAND_MAX）;	
			gauss[i] = gauss[i] + x[n];		//gauss范围-N_perpoint~N_perpoint
		}
		//使gauss有正有负，对半分，gauss范围-N_perpoint/2~N_perpoint/2
		gauss[i] = gauss[i] - （N_perpoint/2）;
		gauss[i] = gauss[i]/（N_perpoint）;	//gauss范围-1/2~1/2
	}
}

void dft（double x[] int len double rx[] double ix[] double y[]）//离散傅里叶变换
{
	int k n;
	for （k=0; k	{
		rx[k] = 0;
		ix[k] = 0;
		y[k] = 0;
		for （n=0; n		{
			rx[k] = rx[k] + x[k]*cos（2*PI*n*k/len）;
			ix[k] = ix[k] - x[k]*sin（2*PI*n*k/len）;
		}
		y[k] = sqrt（rx[k]*rx[k] + ix[k]*ix[k]）;
	}
}

void bp_gauss_gen（double rh[] double ih[]）//通过已有的高斯噪声序列，产生带通高斯白噪声
{
	int i;
	for （i=0; i<500; i++）
	{
		rh[i] = 0;
		ih[i] = 0;
	}
	for （i=1501; i	{
		rh[i] = 0;
		ih[i] = 0;
	}
}

//频域相乘
void multiply（double rx[] double ix[] double rh[] double ih[] double real[] double imag[]）
{
	int i;
	for（i=0;i<（N_gauss+N_h）;i++）
    {
		real[i]=rx[i]*rh[i]-ix[i]*ih[i];
		imag[i]=rx[i]*rh[i]+ix[i]*rh[i];
    }
}

void idft_（double rx[] double ix[] double y[]）//离散傅里叶逆变换
{
	int		k n length;
	length = N_gauss+N_h;
	double	ry[（N_gauss+N_h）] iy[（N_gauss+N_h）];
	for （k=0; k	{
		ry[k] = 0;
		iy[k] = 0;
		y[k] = 0;
		for （n=0; n		{
			ry[k] = ry[k] + rx[n]*cos（2*PI*n*k/length）-ix[n]*sin（2*PI*n*k/length）;
			iy[k] = iy[k] + rx[n]*sin（2*PI*n*k/length）+ix[n]*cos（2*PI*n*k/length）;
		}
		y[k] = sqrt（ry[k]*ry[k] + iy[k]*iy[k]）;
	}
}

void printans（double x[] bool type）//各种输出
{
	int i len;
	char* str;
	switch（type）
	{
		case 0:						//如果输出gauss
			str = “gauss“;
			len = N_gauss;
			break;
		case 1:						//如果输出ans
			str = “ans“;
			len = N_gauss+N_h;
			break;
	}
	for （i=0; i	{
		if （i%3 == 0）
		{
			cout<		}
		cout<			<			<	}
	if （str == “gauss“）
	cout<}

void output（double ans[]

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件     249900  2009-12-08 21:41  高斯噪声\Debug\gaosi.exe

     文件     305516  2009-12-08 21:41  高斯噪声\Debug\gaosi.ilk

     文件      23665  2009-12-08 21:41  高斯噪声\Debug\gaosi.obj

     文件     341448  2009-12-08 21:03  高斯噪声\Debug\gaosi.pch

     文件     590848  2009-12-08 21:41  高斯噪声\Debug\gaosi.pdb

     文件      50176  2009-12-08 21:41  高斯噪声\Debug\vc60.idb

     文件      61440  2009-12-08 21:41  高斯噪声\Debug\vc60.pdb

     文件       4024  2009-12-08 21:41  高斯噪声\gaosi.cpp

     文件       3389  2009-12-08 21:25  高斯噪声\gaosi.dsp

     文件        535  2009-12-08 22:43  高斯噪声\gaosi.dsw

     文件      33792  2009-12-08 22:43  高斯噪声\gaosi.ncb

     文件      48640  2009-12-08 22:43  高斯噪声\gaosi.opt

     文件        741  2009-12-08 21:41  高斯噪声\gaosi.plg

     目录          0  2009-12-08 21:41  高斯噪声\Debug

     目录          0  2009-12-08 22:43  高斯噪声

----------- ---------  ---------- -----  ----

              1714114                    15