#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.存放濾波器輸入序列；分塊處理時，用于表示當前塊內的濾波器的輸入序列；
y:長度len,存放濾波器輸出序列；分塊處理時，用于表示當前塊內的濾波器的輸出序列；
len:輸入序列和輸出序列長度；分塊處理時，表示塊的長度；
px:體積ns×（n＋1）；分塊處理時，用于保存前一塊濾波時的（n＋1）個輸入序列值px[j][]={x(k),x(k-1),...,x(k-m)}；
py:體積ns×（n＋1）；分塊處理時，用于保存前一塊濾波時的n個輸出序列值py[j][]={y(k-1),y(k-2),...,y(k-n)}；
px、py為了防止x[k]過長而設，要初始化為0.*/
void filterp(b,a,n,ns,x,y,len,px,py)
int n,ns,len;
double a[],b[],x[],y[],px[],py[];
{	int i,j,k,nl;
	double sum;
	nl=n+1;
	for(k=0;k<len;k++)  /*初始化y[];*/
	{	y[k]=0.0;	}
	for( j=0;j<ns;j++)
	{	for(k=0;k<len;k++)
		{	px[j*nl+0]=x[k];   /*存放濾波器輸入序列*/
			sum=b[j*nl+0]*px[j*nl+0];  /*y(k)=b(j,0)*x(k);*/
			for(i=1;i<=n;i++)
			{	sum+=b[j*nl+i]*px[j*nl+i]-a[j*nl+i]*py[j*nl+i]; /*y(k)+=b(j,i)*x(k-i)-a(j,i)*y(k-i); i=1...n;*/
			}
			/*檢驗輸出序列x[]的穩定性;*/
			if(fabs(sum)>1.0e10)
			{	printf("This is an unstable filter!\n");
				exit(0);
			}
			/*px，py向右平移一位；*/
			for(i=n;i>=2;i--)
			{	px[j*nl+i]=px[j*nl+i-1];
				py[j*nl+i]=py[j*nl+i-1];
			}
			px[j*nl+1]=px[j*nl+0]; 
			py[j*nl+1]=sum; /*存放運算后的輸出序列x[k];*/
			y[k]=y[k]+sum; 
		}
	}
}