gainc.c

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

#include"math.h"
/*b:双精度实型二维数组，长度为ns*（n+1）,存放滤波器分子多项式系数b[j][i]表示第j个n阶节的分子多项式的第i个系数；
a:双精度实型二维数组，长度为ns*（n+1）,存放滤波器分母多项式系数a[j][i]表示第j个n阶节的分子多项式的第i个系数;
n:滤波器每节的阶数；
ns:滤波器的n阶节数：L；
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）*/
void gainc(b,a,n,ns,x,y,len,sign)
int n,ns,len,sign;
double b[],a[],x[],y[];
{	int i,j,k,nl;
	double ar,ai,br,bi,zr,zi,im,re,den,numr,numi,freq,temp;
	double hr,hi,tr,ti;
	nl=n+1;
	for(k=0;k<len;k++)
	{	/*freq为归一化频率*/   
		freq= k*0.5/(len-1); /*GAINC EXAMPLE*/
		/*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);*/                  
		x[k]=1.0;
		y[k]=0.0;
		for(j=0;j<ns;j++)
		{	br=0.0;
			bi=0.0;
			for(i=n;i>0;i--)
			{	re=br;
				im=bi;
				br=(re+b[j*nl+i])*zr-im*zi;
				bi=(re+b[j*nl+i])*zi+im*zr;
			}
			ar=0.0;
			ai=0.0;
			for(i=n;i>0;i--)
			{	re=ar;
				im=ai;
				ar=(re+a[j*nl+i])*zr-im*zi;
				ai=(re+a[j*nl+i])*zi+im*zr;
			}
			br=br+b[j*nl+0];
			ar=ar+1.0;
			numr=ar*br+ai*bi;
			numi=ar*bi-ai*br;
			den=ar*ar+ai*ai;
			hr=numr/den;
			hi=numi/den;
			tr=x[k]*hr-y[k]*hi;
			ti=x[k]*hi+y[k]*hr;
			x[k]=tr;
			y[k]=ti;
		}
		switch(sign)
		{	case 1:
			{	temp=sqrt(x[k]*x[k]+y[k]*y[k]);
				if(temp!=0.0)
				{	y[k]=atan2(y[k],x[k]);	}
				else 
				{	y[k]=0.0;	}
				x[k]=temp;
				break;
			}
			case 2:
			{	temp=x[k]*x[k]+y[k]*y[k];
				if(temp!=0.0)
				{	y[k]=atan2(y[k],x[k]);	}
				else
				{	temp=1.0e-40;
					y[k]=0.0;
				}
				x[k]=10.0*log10(temp);
			}
		}
	}
}