filterp example.c

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

#include"stdio.h"
#include"filterp.c"
/*求H(z)的单位冲激响应；
H(z)=(0.08327+0.0239z^-1)/(1-1.5658z^-1+0.6549z^-2)
	  -(0.08327+0.0246z^-1)/(1-1.4934z^-1+0.8392z^-2) */
void main()
{	int i,j,n,ns,len,nblk;
	double b[2][3]={{0.08327,0.0239,0.0},{-0.08327,-0.0246,0.0}};
	double a[2][3]={{1.0,-1.5658,0.6549},{1.0,-1.4934,0.8392}};
	double px[2][3],py[2][3];
	double x[25],y[25],data[100];
	FILE *fp;
	n=2;	/*每节的阶数*/
	ns=2;   /*节数L*/
	len=25;	/*输入和输出序列长度*/
	nblk=4;	/*块*/
	/*输入长度为len的单位冲激响应*/
	x[0]=1.0;
	for(i=1;i<len;i++)
	{	x[i]=0.0;	}
	/*初始化px,py;*/
	for(i=0;i<ns;i++)
	for(j=0;j<=n;j++)
	{	px[i][j]=0.0;
		py[i][j]=0.0;
	}
	/*计算并联型IIR数字滤波；*/
	for(j=0;j<nblk;j++)
	{	filterp(b,a,n,ns,x,y,len,px,py);
		for(i=0;i<len;i++)
		{	data[j*len+i]=y[i];
			x[i]=0.0;
		}
	}
	/*输出结果*/
	printf("Unit Impulse Response.\n");
	if(( fp=fopen("filterp.txt","w"))==NULL)
	{	printf("cannot open file 'filterp.txt'!\n");
		exit(0);
	}
	for(i=0;i<100;i++)
	{	fprintf(fp,"%lf ",data[i]); }
	fclose(fp);
}