利用窗函数法设计线性相位FIR滤波器

1、西安工程大学课 程 设 计 说 明 书题目：利用窗函数法设计线性相位FIR滤波器学院（系）：计算机科学学院年级专业：电子信息科学与技术01班学 号：*学生姓名：*日期：20131231-20140103论文摘要数字信号处理学科的一项重大进展是关于数字滤波器设计方法的研究。而FIR数字滤波器可以方便地实现线性相位且其群时延不随频率变化的，在任何幅度特性下都能稳定的工作，因此在数字信号处理中占有非常重要的地位。在现代电子系统中，FIR数字滤波器以其良好的线性特性被广泛使用。FIR数字滤波器传统的设计方法有窗函数法、频率抽样法和等波纹逼近法。用窗函数设计FIR数字滤波器就是用有限长的脉冲相应逼近序列

2、，其基本设计思想为：首先选定一个理想的选频滤波器，然后截取它的脉冲响应得到线性相位。本文是以窗函数设计方法为基础的。首先介绍了满足线性相位条件的FIR数字滤波器的基础知识，简述了FIR数字滤波器的设计方法与流程。本课题主要应用MATLAB软件利用窗函数法设计满足线性相位的FIR数字滤波器，并对所设计的滤波器进行仿真，得到滤波器的单位脉冲响应、频率响应特性。关键字：数字滤波器;C语言。窗口函数设计FIR滤波器论文、窗函数法的基本思想窗函数设计的基本思想是要选取某一种合适的理想频率选择性滤波器，然后将它的脉冲响应截断以得到一个线性相位和因果的FIR滤波器。因此这种方法的重点在于选择某种合适的窗函数

3、和一种理想滤波器。对于给定的滤波器技术指标，选择滤波器长度和具有最窄主瓣宽度和尽可能小的旁瓣衰减的某个窗函数。任何数字滤波器的频率响应都是的周期函数，它的傅立叶级数展开式为： (1-1)式中， 其中的为滤波器的归一化的截止频率。傅立叶系数ha(n)实际上就是理想数字滤波器的冲激响应。获得有限冲激响应数字滤波器的一种可能方法就是把无穷级数截取为有限项级数来近似，而吉布斯(Gibbs)现象使得直接截取法不甚令人满意。窗函数法就是用被称为窗函数的有限加权系列(n)来修正式(1-2)的傅立叶级数，以求得要求的有限冲激响应序列h(n)，即有： (1-2)是有限长序列，当nN-1及n0时，。、窗函数的选择

4、工程中比较常用的窗函数有：矩形窗函数、三角形(Bartlett)窗函数、汉宁(Harming)窗函数、汉明(Hamming)窗函数、布莱克曼(Blackman)窗函数和凯塞(Kaiser)窗函数。这几种窗函数的比较见下表所示：窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134/N-12三角窗-258/N-25汉宁窗-318/N-44海明窗 -418/N-53布莱克曼窗-5712/N-74凯塞窗-5710/N-80六种窗函数的基本参数几种常用窗函数对比，窗函数的选择原则是：具有较低的旁瓣幅度，尤其是第一旁瓣幅度；旁瓣幅度下降速度要大，以利增加阻带衰减；主瓣的宽度要窄，以获得较陡的过

5、渡带。通常上述三点很难同时满足。当选用主瓣宽度较窄时，虽然得到较陡的过渡带，但通带和阻带的波动明显增加：当选用最小的旁瓣幅度时，虽能得到匀滑的幅度响应和较小的阻带波动，但过渡带加宽。因此，实际选用的窗函数往往是它们的折衷。在保证主瓣宽度达到一定要求的条件下，适当牺牲主瓣宽度来换取旁瓣波动的减少。就几种窗函数的比较在本课题中我们采用了汉明窗。1、 矩形窗2、 汉宁窗3、 三角窗4、 海明窗5、 布拉克曼窗6、 凯塞窗其中是I。（）第一类零阶修正贝塞耳函数。口是控制窗函数形状的参数，越大，w（n）窗越窄，频谱的旁瓣越小，但主瓣也相应加宽。口的典型值为4p9。P=0时，凯塞窗变为矩形窗；=544时，

6、凯塞窗与海明窗接近；p=8.5时，凯塞窗与布拉克曼窗接近。、FIR滤波器的设计方法1选择窗函数的类型和长度（1）.根据阻带最小衰减选择窗函数的类型原则是：在保证阻带衰减满足要求的情况下，尽量选择主瓣窄的窗函数。（2）.根据过渡带的宽度选择窗函数的长度2按性能指标要求，构造希望频率响应函数 (1-12) 近似为过渡带中心频率，幅度函数衰减一半（-6dB) (1-13)（1）.确定期望滤波器的单位脉冲响应(1-14)（2）.加窗得到设计结果(1-15)、 利用窗函数设计C语言实现：1子函数语句 vold firwin(n，band，fln，fhn，wn，h) 2形参说明 n整型变量。滤波器的阶数。

7、 band整型变量。滤波器的类型。取值为1、2、3和4，分别对应低通、高通、带通和带阻滤波器。 fln双精度实型变量。 fhn双精度实型变量。 对于低通和高通滤渡器，fln：通带边界频率对于带通和带阻滤渡器，fln：通带下边界频率，fhn：通带上边界频率 wn一整型变量。窗函数的类型；取值l到7，分别对应矩形窗，图基窗，三角窗，汉宁窗，海明窗，布拉克曼窗和凯塞窗。 h双精度实型一维数组，长度为(n+1)。存放FIR滤波器的系数。 函数firwin()调用以下函数： 子函数window()，窗函数的计算。 子函数Kalser()，凯塞窗的计算。子函数bessel0()，贝塞耳函数的计算。子函数程

8、序#includemath.h#includestdio.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=

9、7)/只有凯塞窗需要输入系数betaprintf(input beta parameter of Kaiser window(2beta=3) /*先判断用户输入的数据，如果band参数大于3*/wc2=2.0*pi*fhn;switch(band)case 1:for(i=0;i=n2;i+)s=i-delay;/hi=(sin(wc1*s/fs)/(pi*s)*window(wn,n+1,i,beta);/低通,窗口长度=阶数+1，故为n+1hn-i=hi;if(mid=1)hn/2=wc1/pi;/n为偶数时，修正中间值系数break;case 2: for(i=0;i=n2;i+)s=

10、i-delay;hi=(sin(pi*s)-sin(wc1*s/fs)/(pi*s);/高通-/对hi=hi*window(wn,n+1,i,beta);hn-i=hi;if(mid=1)hn/2=1.0-wc1/pi;/对break;case 3:for(i=0;in2;i+)s=i-delay;hi=(sin(wc2*s/fs)-sin(wc1*s/fs)/(pi*s);/带通-/对hi=hi*window(wn,n+1,i,beta);hn-i=hi;if(mid=1)hn/2=(wc2-wc1)/pi;/对break;case 4:for(i=0;i=n2;i+)s=i-delay;h

11、i=(sin(wc1*s/fs)+sin(pi*s)-sin(wc2*s/fs)/(pi*s);/带阻-/对hi=hi*window(wn,n+1,i,beta);hn-i=hi;if(mid=1)hn/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;switc

12、h(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

13、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=2.0*i/(double)(n-1)-1.0;a2=a*a;beta1=beta*sqrt(1.0-a2);b2=bessel0(beta1);w=b2/b1;return(w);static double bessel0(x)double x;in

14、t i;double d,y,d2,sum;y=x/2.0;d=1.0;sum=1.0;for(i=1;i=25;i+)d=d*y/i;d2=d*d;sum=sum+d2;if(d2=3) printf(”input high cutoff frequency fhn”) scanf(”If”，&fh)； printf(”input sample frequence fsn”)：scanf(”If”，&fs)；printf(”select windown”)；printf(“1 -rectangular；2-tapered rectangularn”)printf(“3 -triangular

15、；4Hanningn”)I printf(“5- Hamming)；6Blaekmann”)；printf(“7-Kaisern”)l scanf(”d”，wn)；fl=fl/fs；fh=fh/fs；firwin(n，band，fl，fh，wn，h)；printf(”FIR digital filtern”)；printf(*f impulse response*nn”); n2=n/2;for(i=0;i=n2；i+)j=n-i； printf(”h(2d)一128lf=h(2d)n”，i，hi，j);printf(“ninput file nsme of frequency respons

16、en”);scanf(”s，fname)；if(fpfopen(fname，”w”)=NULL)printf(”cannot open this filen”)； exit(1)； gain(h,c,n,1,x,y,300,2)； for(i=0；i300；i+) freq=0.5*i/299； fprintf(fp,”lf1fn”,freq,xi)； fclose(fp)；设计一个30阶的FIR高通数字滤波器，其通带边界频率为03。选择参数n=30，hand=2，fl=0.3，fs=1，wn=4，采用海明窗函数。运行结果：h(0)= 0.=h(30) h(1)=-0.=h(29) h(2)=

17、 0. =h(28)h(3)= 0.=h(27)h(4)=-0.=h(26)h(5)=-0.=h(25)h(6)= 0.=h(24)h(7)=-0.=h(23)h(8)=-0.=h(22)h(9)= 0.=h(21)h(10)= 0.=h(20)h(11)=-0.=h(19)h(12)= 0.=h(18)h(13)= 0.=h(17)h(14)=-0.=h(16)h(15)= 0.=h(15)滤波器幅度响应如图：总结与体会设计带通滤波器时首先要计算出过渡带，然后查表得到不同窗函数所需要的阶数，不同的窗函数所设计的滤波器的形状各有差异，尤其在主瓣宽度、旁瓣的形状以及主瓣与旁瓣的高度差上有比较明显得差别，实际应用中应根据实际情况，折衷处理，兼顾各项指标。为了这次课程设计，自己自学了数字信号处理领域中窗