gain.c

自己编的数值计算方法的c语言实现源程序。

#include"math.h"

/*计算数字滤波器的频率响应、幅频响应和相频响应*/
void gain(b,a,m,n,x,y,len,sign)
/*b:双精度实型一维数组，长度为（m+1）,存放滤波器分子多项式系数b(i);
a:双精度实型一维数组，长度为（n+1）,存放滤波器分母多项式系数a(i);
m:分子多项式阶数；
n:分母多项式阶数；
x:长度为len，sign＝0时，存放Re[H(w)];sign=1时，存放|H(w)|;sign=2时，存放用分贝表示的|H(w)|；
y:长度为len，sign＝0时，存放Im[H(w)];sign=1和2时，存放psi（w）；
len: 频率响应长度；
sign:sign＝0时，计算Re[H(w)]和Im[H(w)]；sign＝1时，计算|H(w)|和psi（w）；sign＝2时，计算|H(w)|（dB表示）和psi（w）*/
int m,n,len,sign;
double b[],a[],x[],y[];
{	int i,k;
	double ar,ai,br,bi,zr,zi,im,re,den,numr,numi,freq,temp;
	/*归一化频率为freq的响应放在x[k]、y[k]中*/
	for(k=0;k<len;k++)
	{	/*freq为归一化频率*/  
		/*freq= k*0.5/(len-1); /*ATTENTION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FOR GAIN EXAMPLE*/        
		freq= k*1.0/(len-1); /*FOR GAIN EXAMPLE2*/
		/*w=2*pi*freq 为圆周频率；  exp(-jw)=cos(-w)+j*sin(-w)*/
		zr=cos(-8.0*atan(1.0)*freq); /*zr=cos(-w); */
		zi=sin(-8.0*atan(1.0)*freq); /*zi=sin(-w);*/                       
		
		/*求H（w）的分子实部br,虚部bi；分母实部ar,虚部ai；(由三角函数分解的循环乘法求出)*/
		br=0.0;
		bi=0.0;
		for(i=m;i>0;i--)
		{	re=br;
			im=bi;
			br=(re+b[i])*zr-im*zi;/*求sum_{i=1}^{m} b(i)*cos(i*w);*/
			bi=(re+b[i])*zi+im*zr;/*求sum_{i=1}^{m} b(i)*sin(i*w);*/
		}
		ar=0.0;
		ai=0.0;
		for(i=n;i>0;i--)
		{	re=ar;
			im=ai;
			ar=(re+a[i])*zr-im*zi;/*求sum_{i=1}^{m} a(i)*cos(i*w);*/
			ai=(re+a[i])*zi+im*zr;/*求sum_{i=1}^{m} a(i)*sin(i*w);*/
		}
		br=br+b[0];
		ar=ar+1.0;  /*a[0]=1.0;*/
		
		/*(br+i*bi)/(ar+i*ai)=(numr+i*numi)/den*/
		numr=ar*br+ai*bi;
		numi=ar*bi-ai*br;
		den=ar*ar+ai*ai;
		x[k]=numr/den;
		y[k]=numi/den;
		/*sign=1时，y[k]=psi（w）;sign=2时，y[k]=psi（w）（dB表示）*/
		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);
			}
		}
	}
}