gainc.c

实现了巴特沃斯滤波。用C语言写的。实现了巴特沃斯滤波。用C语言写的。

#include <stdio.h>
#include "math.h"
void gainc(b,a,n,ns,x,y,len,sign)
int n,ns,len,sign;
double b[],a[],x[],y[];
{
	int i,j,k,n1;
	double ar,ai,br,bi,zr,zi,im,re,den,numr,numi,freq,temp;
	double hr,hi,tr,ti;
	n1=n+1;
	for (k=0;k<len ;k++ )
	{
		freq=k*0.5/(len-1);
		zr=cos(-8.0*atan(1.0)*freq);
		zi=sin(-8.0*atan(1.0)*freq);
		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*n1+i])*zr-im*zi;
				bi=(re+b[j*n1+i])*zi+im*zr;
			}
			ar=0.0;
			ai=0.0;
			for (i=n;i>0 ;i-- )
			{
				re=ar;
				im=ai;
				ar=(re+a[j*n1+i])*zr-im*zi;
				ai=(re+a[j*n1+i])*zi+im*zr;
			}
			br=br+b[j*n1+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);

			}
		}
	}
}