www.pudn.com > 滤波图生成系统filter_1.34.zip > calculate.h

#include "math.h" 
//用混合同余法产生（a，b）区间上的随机数yi. 
double uniform(double a,double b,long *seed) 
{ 
	double t; 
	*seed=2045*(*seed)+1; 
	*seed=*seed-(*seed/1048576)*1048576; 
	t=(*seed)/1048576.0; 
	t=a+(b-a)*t; 
	return(t); 
} 
 
double gauss(double mean,double sigma,long *s) 
{ 
	int i;double x,y; 
    for(x=0,i=0;i<12;i++) 
    { 
		x+=uniform(0.0,1.0,s); 
	} 
    x=x-6.0; 
	y=mean+x*sigma; 
	return(y); 
} 
 
void sinwn(double a[],double f[],double ph[],int m,double fs,double snr,long seed,double x[],int n) 
{ 
	int i,k; 
	double z,pi,nsr; 
	pi=4.0*atan(1.0); 
	z=snr/10.0; 
	z=pow(10.0,z); 
	z=1.0/(2*z); 
	nsr=sqrt(z); 
	for(i=0;i0;i--) 
		{ 
			re=br; 
			im=bi; 
			br=(re+b[i])*zr-im*zi; 
			bi=(re+b[i])*zi+im*zr; 
		} 
		ar=0.0; 
		ai=0.0; 
		for(i=n;i>0;i--) 
		{ 
			re=ar; 
			im=ai; 
			ar=(re+a[i])*zr-im*zi; 
			ai=(re+a[i])*zi+im*zr; 
		} 
		br=br+b[0]; 
		ar=ar+1.0; 
		numr=ar*br+ai*bi; 
		numi=ar*bi-ai*br; 
		den=ar*br+ai*br; 
		x[k]=numr/den; 
		y[k]=numi/den; 
		switch(sign) 
		{ 
		case 1: 
			{ 
				temp=sqrt(x[k]*x[k]+y[k]*y[k]); 
				y[k]=atan2(y[k],x[k]); 
				x[k]=temp; 
				break; 
			} 
		case 2: 
			{ 
				temp=x[k]*x[k]+y[k]*y[k]; 
				y[k]=atan2(y[k],x[k]); 
				x[k]=10.0*log10(temp); 
			} 
		} 
	} 
}