实际可运行的FIR滤波器系数设计的C语言程序

#include"math.h"

#include"stdio.h"

/*n:滤波器的阶数,n是奇数时可用来设计各种滤波器

band:滤波器类型1、2、3、4分别对应低通，高通，带通，带阻

fln:通带下边界频率

fhn:通带上边界频率

wn:窗函数的类型1-7分别对应

h:长度为n+1,存放滤波器的系数

fs:采样频率*/

void firwin(n,band,fln,fhn,wn,h)

int n,band,wn;

double fln,fhn,h[];

{

int i,n2,mid;

double s,pi,wc1,wc2,beta,delay,fs;

double window();//窗函数的计算

beta=0.0;

if(wn==7)//只有凯塞窗需要输入系数beta

{printf("input beta parameter of Kaiser window(2

scanf("%1f",&beta);

}

pi=4.0*atan(1.0);//pi=PI;

if((n%2)==0)/*如果阶数n是偶数*/

{n2=n/2+1;/**/

mid=1;//

}

else

{n2=n/2;//n是奇数,则窗口长度为偶数

mid=0;//

}

delay=n/2.0;

wc1=2.0*pi*fln;//

if(band>=3) /*先判断用户输入的数据，如果band参数大于3*/

{wc2=2.0*pi*fhn;}

switch(band)

{

case 1:

{for(i=0;i<=n2;i++)

{s=i-delay;//

h[i]=(sin(wc1*s/fs)/(pi*s))*window(wn,n+1,i,beta);//低通,窗口长度=阶数+1，故为n+1

h[n-i]=h[i];

}

if(mid==1)h[n/2]=wc1/pi;//n为偶数时，修正中间值系数

break;

}

case 2:

{for(i=0;i<=n2;i++)

{s=i-delay;

h[i]=(sin(pi*s)-sin(wc1*s/fs))/(pi*s);//高通-//对

h[i]=h[i]*window(wn,n+1,i,beta);

h[n-i]=h[i];

}

if(mid==1)h[n/2]=1.0-wc1/pi;//对

break;

}

case 3:

{for(i=0;i

{s=i-delay;

h[i]=(sin(wc2*s/fs)-sin(wc1*s/fs))/(pi*s);//带通-//对

h[i]=h[i]*window(wn,n+1,i,beta);

h[n-i]=h[i];

}

if(mid==1)h[n/2]=(wc2-wc1)/pi;//对

break;

}

case 4:

{for(i=0;i<=n2;i++)

{s=i-delay;

h[i]=(sin(wc1*s/fs)+sin(pi*s)-sin(wc2*s/fs))/(pi*s);//带阻-//对

h[i]=h[i]*window(wn,n+1,i,beta);

h[n-i]=h[i];

}

if(mid==1)h[n/2]=(wc1+pi-wc2)/pi;

break;

}

}

}

/*返回窗函数的表达式

n:窗口长度

type:选择窗函数的类型

beta:生成凯塞窗的系数*/

static double window(type,n,i,beta)

int i,n,type;

double beta;

{

int k;

double pi,w;

double kaiser();

pi=4.0*atan(1.0);//pi=PI;

w=1.0;

switch(type)

{

case 1:

{w=1.0;//矩形窗

break;

}

case 2:

{k=(n-2)/10;

if(i<=k)

w=0.5*(1.0-cos(i*pi/(k+1)));//图基窗

break;

}

case 3:

{w=1.0-fabs(1.0-2*i/(n-1.0));//三角窗

break;

}

case 4:

{w=0.5*(1.0-cos(2*i*pi/(n-1)));//汉宁窗

break;

}

case 5:

{w=0.54-0.46*cos(2*i*pi/(n-1));//海明窗

break;

}

case 6:

{w=0.42-0.5*cos(2*i*pi/(n-1))+0.08*cos(4*i*pi/(n-1));//布莱克曼窗

break;

}

case 7:

{w=kaiser(i,n,beta);//凯塞窗

break;

}

}

return(w);

}

static double kaiser(i,n,beta)//因凯塞窗的表达式复杂，调用贝塞尔窗表达式

int i,n;

double beta;

{

double a,w,a2,b1,b2,beta1;

double bessel0();

b1=bessel0(beta);

a