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

				#include"math.h"
				/*b,a,m,n:为滤波器参数； x,len:先为输入序列，后为输出序列参数；px,py:0;*/
				void filter(b,a,m,n,x,len,px,py)
				/*b:双精度实型一维数组，长度为（m+1）,存放滤波器分子多项式系数b(i);
				a:双精度实型一维数组，长度为（n+1）,存放滤波器分母多项式系数a(i);
				m:滤波器分子多项式阶数；
				n:滤波器分母多项式阶数；
				x:长度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[],px[],py[];
				{	int k,i;
					for(k=0;k					{	px[0]=x[k];
						x[k]=0.0;
						for(i=0;i						{	x[k]=x[k]+b[i]*px[i];/*px[i]==x[k-i];i=0...m;*/
						}
						for(i=1;i						{	x[k]=x[k]-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]=x[k];
					}
				}