gaussiannoisegenerator.cpp

QPSK调制Modem的C实现及文档说明

/* This file is written in C. This is obtained from Jean Au. */

#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 134456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349

double ran2(long  *idum)
{
 static long ix1,ix2,ix3;
 static double r[98];
 double temp;
 static int iff=0;
 int j;

   if (*idum<0 || iff ==0){
	 iff=1;
	 ix1=(IC1-(*idum)) % M1;
	 ix1=(IA1*ix1+IC1) % M1;
	 ix2=ix1 % M2;
	 ix1=(IA1*ix1+IC1) % M1;
	 ix3=ix1 % M3;
	 for (j=1;j<=97;j++) {
	   ix1=(IA1*ix1+IC1) % M1;
	   ix2=(IA2*ix2+IC2) % M2;
	   r[j]=(ix1+ix2*RM2)*RM1; 
	  }
	  *idum=1; 
    }
    ix1=(IA1*ix1+IC1) % M1;
    ix2=(IA2*ix2+IC2) % M2;
    ix3=(IA3*ix3+IC3) % M3;
    j=1+ ((97*ix3)/M3);
    temp=r[j]; 
    r[j]=(ix1+ix2*RM2)*RM1; 
    return temp; 
} 

/*********************************************************************/

double Gaussian(long *idum2)
{ 
 static int iset=0; 
 static double gset;
 double fac,r,v1,v2; 
 double ran1();

   if (iset ==0) { 
	 do { 
		v1=2.0*ran2(idum2)-1.0;
		v2=2.0*ran2(idum2)-1.0;
		r=v1*v1+v2*v2; 
	 }  while (r>= 1.0 || r==0.0); 
	 fac=sqrt(-2.0*log(r)/r); 
	 gset=v1*fac; 
	 iset=1; 
	 return v2*fac; 
	}else { 
	  iset=0; 
	  return gset; 
	} 
}