#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<len;k++)
	{	px[0]=x[k];
		x[k]=0.0;
		for(i=0;i<=m;i++)
		{	x[k]=x[k]+b[i]*px[i];/*px[i]==x[k-i];i=0...m;*/
		}
		for(i=1;i<=n;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];
	}
}