filtery.c

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

#include"math.h"
/*b,a,m,n:为滤波器参数； x,len:先为输入序列，后为输出序列参数；px,py:0;*/
void filtery(b,a,m,n,x,y,len,px,py)
/*b:双精度实型一维数组，长度为（m+1）,存放滤波器分子多项式系数b(i);
a:双精度实型一维数组，长度为（n+1）,存放滤波器分母多项式系数a(i);
m:滤波器分子多项式阶数；
n:滤波器分母多项式阶数；
x:长度len.存放滤波器输入序列，分块处理时，用于表示当前块内的滤波器的输入序列；
y:长度len.存放滤波器输出序列，分块处理时，用于表示当前块内的滤波器的输出序列；
len:输入序列和输出序列长度；分块处理时，表示块的长度；
px:长度（m＋1）；分块处理时，用于保存前一块滤波时的（m＋1）个输入序列值px[]={x(k),x(k-1),...,x(k-m)}；
py：长度（n＋1）；分块处理时，用于保存前一块滤波时的n个输出序列值py[]={y(k-1),y(k-2),...,y(k-n)}；
px、py为了防止x[k]过长而设，要初始化为0.*/
int m,n,len;
double a[],b[],x[],y[],px[],py[];
{	int k,i;
	double sum;
	for(k=0;k<len;k++)
	{	px[0]=x[k];
		sum=0.0;
		for(i=0;i<=m;i++)
		{	sum=sum+b[i]*px[i];/*px[i]==x[k-i];i=0...m;*/
		}
		for(i=1;i<=n;i++)
		{	sum=sum-a[i]*py[i];/*py[i]==y[k-i];i=1...n;*/
		}
		/*检查滤波器是否稳定；*/
		if(fabs(x[k])>1.0e10)
		{	printf("This is an unstable filter!\n");
			exit(0);
		}
		/*px,py中数据都向右移一格；*/
		for(i=m;i>=1;i--)
		{	px[i]=px[i-1];
		}
		for(i=n;i>=2;i--)
		{	py[i]=py[i-1];
		}
		py[1]=sum;
		y[k]=y[k]+sum;
	}
}