random.cpp

给语音加噪声的一个程序

#include "random.h"



double ran0(long *ideum)
    //Park & Miller divelope "the lowest" randnum genitor 
	//  it return 0.0-1.0 
{
	long k;
	double ans;

	*ideum^=MASK;
	k=(*ideum)/IQ;
	*ideum=IA*(*ideum-k*IQ)-IR*k;
	if(*ideum<0) *ideum+=IM;
	ans=(double)AM*(*ideum);
	*ideum^=MASK;
	return ans;
}

double ran1(long *ideum)
   //Park & Miller 's "the lowest " randnum genitor 
   // +Bays-Durham shaffle & protect 
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;

	if(*ideum<=0||!iy){
		if(-(*ideum)<1) *ideum=1;
		else *ideum=-(*ideum);
		for(j=NTAB+7;j>=0;j--){
			k=(*ideum)/IQ;
			*ideum=IA*(*ideum-k*IQ)-IR*k;
			if(*ideum<0) *ideum+=IM;
			if(j<NTAB) iv[j]=*ideum;
		}
		iy=iv[0];
	}
			k=(*ideum)/IQ;
			*ideum=IA*(*ideum-k*IQ)-IR*k;
			if(*ideum<0) *ideum+=IM;
			j=iy/NDIV;
			iy=iv[j];
			iv[j]=*ideum;
			if((temp=AM*iy)>RNMX) return RNMX;
			else return temp;
}


double ran2(long *ideum)
//L'Ecuyer 's long randnum (>2e18)
{
	int j;
	long k;
	static long ideum2=123456789;
	static long iy=0;
	static long iv[NTAB];
	double temp;
	if(*ideum<0) {
		if(-(*ideum)<1) *ideum=1;
		else *ideum=-(*ideum);
		ideum2=(*ideum);
		for(j=NTAB+7;j>=0;j--){
			k=(*ideum)/IQ1;
			*ideum=IA1*(*ideum-k*IQ1)-k*IR1;
			if(*ideum<0) *ideum+=IM1;
			if(j<NTAB) iv[j]=*ideum;
		}
		iy=iv[0];
	}
		k=(*ideum)/IQ1;
		*ideum=IA1*(*ideum-k*IQ1)-k*IR1;
		if(*ideum<0) *ideum+=IM1;
		k=ideum2/IQ2;
		ideum2=IA2*(ideum2-k*IQ2)-k*IR2;
		if(ideum2<0) ideum2+=IM2;
		j=iy/NDIV1;
		iy=iv[j]-ideum2;
		iv[j]=*ideum;
		if(iy<1) iy+=IMM1;
		if((temp=AM1*iy)>RNMX) return RNMX;
		 else return temp;
}


double gasdev(long *ideum)
{
	static int iset=0;
	static double gset;
	double fac,rsq,v1,v2;
	if(iset==0){
		do{
			v1=2.0*ran1(ideum)-1.0;
			v2=2.0*ran1(ideum)-1.0;
			rsq=v1*v1+v2*v2;
		}while(rsq>=1.0||rsq==0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	}else{
		iset=0;
		return gset;
	}
}