filterp.c

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

#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; 
		}
	}
}