数字信号处理滤波器设计

1、滤波器设计课程设计题目一、FIR（选取合适的窗函数及窗长度）通带06000Hz通带波动士0.01阻带8500HzFs/2阻带波动-70dB抽样频率44100Hz二、IIR（可选巴特沃斯、切比雪夫I型、切比雪夫II型）1、设计数字低通滤波器，满足以下指标：通带边带频率2000Hz通带波动0.01阻带边界频率2700Hz11025Hz阻带波动0.05抽样频率要求:1、求出h(n)和H(ej)4a)用H(ej)直接画出幅度响应曲线(dB)及相位响应曲线；(b)做h(n)的FFT,有此画出幅度响应曲线(dB)及相位响应曲线；(c)并两种方法的结果进行比较。(d)验证你的设计结果满足设计指标。FIR和I

2、IR滤波器设计过程FIR滤波器的设计：1.设计步骤:(1)写出理想低通滤波器的传输函数和单位脉冲响应。?(n=H?/?>?=?2?.-?1*sin?夕?？一、2?-?(2)根据题目已有的模拟滤波器频率进行归一化求得数字滤波器的频率，并对题目所给的通带波动和阻带波动进行转换，可求得通带截止频率？=0.8549?阻带截止频率？?=1.2110?通带波动？?=0.174?阻带衰减？?=70?(3)选择用布莱克曼窗函数来设计FIR滤波器，通过查阅资料可知道布莱克曼11兀，、,，窗的过渡市范为,最小我减为74?则可以通过matlab中ceil函数向上取整的方式以及阻带频率和通带频率之差所求的过渡带

3、频率来确定滤波器的长度，在经过调试可以取得N=99。当然，在设计中可以之间用blackman()这个函数,但是在这里我们会利用布莱克曼窗函数的公式构造布莱克曼窗：2?4?wn=0.42-0.5cos而%+0.08cos与“?(?？(4)通过求得的矩形窗函数和单位脉冲响应？?？n进行时域上的乘积即可求得FIR滤波器的时域响应，再利用DTFT直接求出H(?X并且可以画出幅度响应曲线(dB)及相位响应曲线OOX?严?=DTFTxn=?=-CO(5)根据时域所求的函数做FFT并由此画出幅度响应曲线(dB)及相位响应曲线，对这两种方法的结果进行比较并验证设计结果是否满足设计指标。3.设计过程：(1)求出

4、h(n),并做h(n)的DTFT从而直接求出H(?的幅度曲线响应(dB)和相位响应曲线。由设计步骤中可以知道？?？n=吗箸以及自己通过公式所构造的布莱克曼窗函数即可求得h(n)。然后将h(n)经过DTFT即可求得H(?,进而求得幅度曲线响应和相位响应曲线。波形如下：0.90.BD1D6邕0三D.40.30.2卜广9/1j1VItRPH*11Vj/L-X尸11-_11020304Q6060703D90100Hi图一布莱克曼窗函数的波形图二理想滤波器冲击响应hdn的波形可|JJJ1L1L11国4口暹0-1001020304050n图三h(n)的波形?的前100个值0.00010.0003Colum

5、ns1through13-0.0000-0.0000-0.0000-0.00000.00000.0000-0.0001-0.0002-0.00010.00020.0004Columns14through26-0.0003-0.0008-0.00050.00050.00140.0010-0.0008-0.0023-0.00160.00120.00360.0026-0.0016Columns27through39-0.0053-0.00410.00210.00770.0062-0.0026-0.0109-0.00930.00310.01540.0139-0.0035-0.0220Columns40

6、through52-0.02110.00390.03270.0340-0.0043-0.0549-0.06480.00450.13910.27290.32880.27290.1391Columns53through650.0045-0.0648-0.0549-0.00430.03400.03270.0039-0.0211-0.0220-0.00350.01390.01540.0031Columns66through78-0.0093-0.0109-0.00260.00620.00770.0021-0.0041-0.0053-0.00160.00260.00360.0012-0.0016Colu

7、mns79through91-0.0023-0.00080.00100.00140.0005-0.0005-0.0008-0.00030.00030.00040.0002-0.0001-0.0002Columns92through99-0.00010.00000.00010.0000-0.0000-0.0000-0.0000-0.0000山P郛祖垣百顺蹙31图四H(?的幅度响应曲线图四H(?的相位响应曲线(2)根据时域所求的函数h(n)做FFT并由此画出幅度响应曲线(dB)及相位响应曲线吗器窃冏誓联譬1111rA0701vn旧I111ll1111iii1I1一一!i一J0,3856-乃dBl1

8、ii1j:-:|-i1II!和_11111II14ii!1i1r|Jiffii>1>iiH1i11111111图五H(?的幅度响应曲线图六H(?的相位响应曲线由题目可知滤波器的指标：通带波动为0.174dB,阻带衰减为-70dB,通带截止频率为0.8549Hz,阻带截止频率为1.2110Hz。通过DTFT直接求H(?所得到的FIR滤波器的指标是：通带波动为0.001dB,阻带衰减为-73.14dB。通过FFT求H(?7所得到的FIR滤波器的指标是：通带波动为0.0014dB,阻带衰减为-75dB1.51noci3"1j111II1111mi1h£u1£

9、；111JQJ。名11.2141.6132w?Hz图七H（?了的幅度响应曲线（|H?/?|）由图六可以知道这个滤波器两个正负尖峰频率差值为：1.2612-0.80381=0.45739Hz设计的的FIR滤波器的过渡带宽为:1.2110-0.8549=0.3561Hz.I1!1!，4，4一1I*41-1-I-.II-R4riaJ.4Bu4BJ.lld.n.ril->»»»»sKMIW<1*-20r-g亦回加顿g例皿吟+HFT来第应的曲为卜卜，曲珀应的曲线图八DTFT和FFT所求的幅度响应（dB）曲线由该图可以看出DTFT所求幅度响应（dB）的曲

10、线和FFT所求的幅度响应(dB)的曲线是基本一致的D1TT求啊应的曲线OFFT未晌应的自装一，口0050101502025030350404505,JV/p4图八DTFT和FFT所求的相位响应曲线综上所述：1.所设计的布莱克曼窗FIR滤波器在阻带衰减上接近了-70(dB)的指标，但是在通带波动上0.02(dB)略小于所要求的0.174(dB)的指标。2.FIR滤波器的过渡带宽小于两个正负尖峰频率差值，满足要求。3.DTFT和FFT两种方法所设计的滤波器结果及指标基本没有区别，但是DTFT的计算量会比FFT相对大，而FFT所求的幅度和相位响应的波形是对称的，需要经过处理才能达到要求。4.DTFT

11、所算出来的指标会相对来说会更接近于题目中要求的指标。IIR滤波器的设计：1.设计步骤：(1)选择采用冲激响应不变法设计巴特沃斯低通滤波器，已知通带边带频率为2000Hz,阻带边带频率为2700Hz,通带波动为0.01,阻带波动为0.05,抽样频率为11025Hz,根据对频率的归一化以及对阻带波动通带波动转换成(dB)的形式可得，Wp=1.1398Hz,Ws=1.5387Hz,Rp=0.087(dB),As=26(dB)。(2)计算样本模拟滤波器的阶次100.1?3?-1?(lg100.1?.1)/2lg坛？可求得N=17o(3)发现Rp与(dB),?月?，由于巴特沃斯滤波器归一化低通原型的通带

12、截止频率？?=1,去归一化成任意截止频率时必须是3(dB)衰减处的？为利用通带或者阻带满足衰减指标的公式去求得3(dB)衰减的截止频率。?)？/?100.1?3?_12?”?/100.1?-1在这里选择的是阻带满足衰减指标的公式，可求得？?=3.8003e+03。(4)根据？?利阶数N求模拟低通滤波器的系统函数分子分母系数。(5)冲激响应不变法利用模拟低通滤波器的系统函数和采样频率转换成数字低通滤波器的系统函数。(6)利用数字低通滤波器的系统函数的分子分母系数求得冲激响应即数字低通滤波器的时域响应。(7)根据所求的时域响应作FFT以求得在频率上的响应。(8)检测步骤(5)中的求得的模拟滤波器系

13、统函数分子分母系数以及数字滤波器系统函数分子分母系数来检验衰减指标。3.设计过程：(1)将要求的各种指标代入巴特沃斯函数中以求得阶数和3(dB)衰减的截止频率，并由此确定模拟滤波器系统函数的分子分母系数，再利用impinvar将模拟滤波器系统函数的分子分母系数转换成数字滤波器系统函数的分子分母系数。然后利用freqz()函数求出数字滤波器频域上的幅度响应曲线(dB)和相位响应曲线。i聃字汨转器帽阴20萼髭祖归等蚂强OoQQQQ口4&凛，K4J1-1144-图一H(?了的幅度响应曲线II1!,飞.iB-_.V-一I1i11o我字派技界珅吊图二H(?的相位响应曲线(2)将(1)中所求的数字

14、滤波器系统函数的分子分母的系数利用imz函数即可求得系统的冲激响应即时域函数h(n)。06电也扫u珀应图三h(n)的波形?的前100个值Columns1through13-0.0000-0.00000.00000.00000.00040.00510.03090.10830.24340.36100.33520.1315-0.1050Columns14through26-0.1764-0.05330.09240.0972-0.0112-0.0796-0.03560.04150.0494-0.0062-0.0414-0.01650.0239Columns27through390.0250-0.006

15、0-0.0226-0.00670.01440.0124-0.0050-0.0123-0.00220.00860.0060-0.0037-0.0066Columns40through52-0.00040.00500.0028-0.0026-0.00350.00030.00290.0012-0.0017-0.00180.00040.00160.0005Columns53through65-0.0010-0.00090.00040.00090.0002-0.0006-0.00040.00030.00050.0000-0.0004-0.00020.0002Columns66through780.000

16、3-0.0000-0.0002-0.00010.00010.0001-0.0000-0.0001-0.00000.00010.0001-0.0000-0.0001Columns79through91-0.00000.00000.0000-0.0000-0.000000000000Columns92through1040000000000000(3)根据时域所求的函数h(n)做FFT并由此画出幅度响应曲线(dB)及相位响应曲线。图四H(?的幅度响应曲线图五H(?的相位响应曲线由题目可知滤波器的指标：通带波动为0.087dB,阻带衰减为26dB,通带截止频率为1.1398Hz,阻带截止频率为1.5

17、387Hz。由于求利用impz函数求h(n)只是数字滤波器的单位冲激响应逼近模拟滤波器的单位冲激响应的过程，再进行FFT变换所得到的数字滤波器的指标误差相对较大，所以在这里只检验直接求频域响应所得到的数字滤波器的指标。利用书上检验衰减指标的代码并将参数代入可求得：通带波动为0.0634dB,阻带衰减为26dB,非常接近题目中要求的指标。20-T-1IIIVI直接啊G|-FFT/由一的旧线，图六直接求和FFT所求的幅度响应(dB)曲线京田团湾吊_收图七直接求和FFT所求的相位响应曲线由以上两幅图可以知道直接求所得到的幅度响应曲线及相位响应曲线在w/pi=00.7这个范围内的波形基本保持一致，但是

18、在w/pi>0.7之后FFT所得到的幅度响应的曲线的谐波较多，而且呈现的也不是线性相位。综上所述，直接求H(?3所设计的低通滤波器基本能够满足题目要求的指标，在比较了直接求H(?和FFT设计的数字低通滤波器后，认为在设计低通巴特沃斯滤波器上还是用直接求H(?的方式较为精确。设计过程中遇到的问题和解决方法：1 .FIR滤波器中在求窗长度的时候11pi/(ws-wp)得不到整数。解决方法：用ceil函数向上取整的方式得到整数。2 .在做FFT变换的时候发现结果是呈对称的形式而且曲线呈现的波形不够完整。解决方法：首先将波形的范围延长一倍，然后采用想要的x的范围去限制matlab画出的y的波形，

19、就可以得到想要的一般的波形。3 .如何克服矢量维度不一致的问题解决方法：先看两个矢量的维度是多少，然后可以在那个维度少的矢量后边用补零函数进行调整使得两个函数的矢量一致。心得体会：通过这次设计FIR滤波器和IIR滤波器的经历，我真心觉得自己的matlab水平真的有了一个质的飞跃，我承认网络上有很多类似的资料，但是要真正达到题目的要求还是需要自己动脑去做的。感觉这个过程就像是建立模型和代码实现的过程，我也是按照书本上两个滤波器的模型来设计并且实现的，整个过程发现了不少的问题，有些是关于课程的问题，比如加窗类型的选择以及调试，在达到过渡带小于正负尖峰之间那个指标的过程真的是调了很久，最后才确立了布

20、莱克曼窗，参数也是改了好多遍。有些是关于matlab的基础问题，因为我接触matlab的时间不多，一开始比较吃力，但是我觉得经过了这个项目之后确实对matlab有了进一步的认识，其实挺希望这个项目能稍微增加难度然后放在课程一开始让学生研究，我觉得这样的效果会更好。程序源代码：FIR滤波器closeallf1=6000;%通带截止频率f2=8500;%阻带截止频率f=44100;%抽样频率wp=2*pi*f1/f%归一化ws=2*pi*f2/falphap=0.174;%通带波动dBalphas=70;%阻带衰减dBc=11*pi;%Blackman窗过渡带宽wc=(ws+wp)/2;delt_

21、w=ws-wp;M=ceil(c/delt_w);N=1*M+1%窗长度n=0:M;%利用课本公式得到布莱克曼窗win=(0.42-0.5*cos(2*pi*n)/(N-1)+0.08*cos(4*pi*n)/(N-1);win1=win.'figure(1)n=-0.5*M:0.5*M;%低通滤波器的时域响应hd=sin(wc*n)./(pi*n);hd(find(n=0)=wc*cos(wc*0)/pi;ht=hd.*win1'%加窗过程plot(n,ht,'.-')%画冲激响应系数图xlabel('n','FontSize'

22、,12);ylabel('ht','fontsize',12);gridonfigure(2)%h,w=freqz(ht,1,1024);%离散时间傅里叶变换得到幅度%W=w/pi;H=20*log10(abs(h);%holdon%plot(W,H);k=1:1024;n=0:M;w=(pi/1024)*k;W=w/pi;h=ht*(exp(-j*pi/1024)A(n'*k);%DTFT变换plot(w/pi,20*log10(abs(h);xlabel('w/pi');ylabel('幅度响应曲线dB');%直接求所

23、得到的频率幅度响应gridonholdonline(0,1,-70,-70,'color','r');%标出相应的限制范围line(0,1,-0.174,-0.174,'color','r');dotp=round(mean(find(w>wp-0.1&w<wp+0.1)%找到点wp的坐标dots=round(mean(find(w>ws-0.1&w<ws+0.1)%找到点ws的坐标plot(W(dotp),H(dotp),'.r','MarkerSize',

24、10);%对应处打上显目的红点plot(W(dots),H(dots),'.r','MarkerSize',10);text(W(dotp),H(dotp),'','',num2str(round(W(dotp)*1000)/1000),.',',num2str(round(H(dotp)*1000)/1000),'dB','FontSize',15);%标示wp在频响图%中的位置，从中可以看来相应增益text(W(dots),H(dots),'',''

25、;,num2str(round(W(dots)*100)/100),.',',num2str(round(H(dots)*100)/100),'dB','fontsize',15);%标示ws在频响图%的位置holdofffigure(3)plot(W,angle(h);axis(0,0.5,-4,4);xlabel('w/pi');ylabel('相位响应曲线');%直接求所得到的相位响应曲线gridonfigure(4)y=20*log10(abs(fft(ht,1024);%FFT求得的频域响应x=(0:10

26、23)/(512);y1=y(x>=0&x<=1);plot(x(x>=0&x<=1),y1);%对所求的对称图形进行处理xlabel('w/pi');ylabel('幅度响应曲线dB');%FFT求得的频域响应曲线gridonholdon;plot(0.3856,-75,'.r','MarkerSize',10);%标出plot(0.2722,-0.0014,'.r','MarkerSize',10);text(0.3856,-75,'0.3856&

27、#39;,'-75','dB','FontSize',15);text(0.2722,-0.0014,'0.2722','-0.0014','dB','FontSize',15);line(0,1,-70,-70,'color','r');%标出相应的限制范围line(0,1,-0.174,-0.174,'color','r');figure(5)y2=angle(fft(ht,1024);x=(0:1023)/(512

28、);y3=y2(x>=0&x<=1);plot(x(x>=0&x<=1),y3);axis(0,0.5,-4,4);xlabel('w/pi');ylabel('相位响应曲线');%FFT求得的相位响应曲线gridonfigure(10)plot(W,H,'*r',x(x>=0&x<=1),y1,'ob');legend('DTFT求响应的曲线','FFT求响应的曲线);将两种方法曲线画在图上比较xlabel('w/pi');yla

29、bel('幅度响应曲线dB');figure(11)plot(W,angle(h),'.r',x(x>=0&x<=1),y3,'ob');axis(0,0.5,-4,4);legend('DTFT求响应的曲线','FFT求响应的曲线);将两种方法曲线画在图上比较xlabel('w/pi');ylabel('相位响应曲线');figure(12)plot(w,abs(h);%标出幅度曲线在正负尖峰的两个点，检测是否满足过渡带的要求holdonaxis(0,2,-0.5,1.

30、5);plot(0.8038,1,'.r','MarkerSize',10);plot(1.261,2.88e-05,'.r','MarkerSize',10);text(0.8038,1,'0.8038','1','','FontSize',15);text(1.261,2.88e-05,'1.261','2.88e-05','','FontSize',15);xlabel('w/Hz'

31、);ylabel('H幅度曲线')；figure(13)plot(n,win1);%画出窗函数的图像gridon;xlabel('n');ylabel('w(n)');figure(14)n=0:M;plot(n,hd);%画出理想滤波器冲激响应的图像gridon;xlabel('n');ylabel('hd(n)');IIR滤波器clc;clearallOP=2*pi*2000;%模拟通带截止频率OS=2*pi*2700;%模拟阻带截止频率Rp=0.087;%通带波动dBAs=26;%阻带衰减dBFs=11025

32、%采样频率Wp=OP/Fs;%归一化Ws=OS/Fs;N,OC=buttord(OP,OS,Rp,As,'s');%求出阶数和3dB衰减时的截止频率%模拟滤波器系统函数分子分母%模拟滤波器系统函数分子分母%数字滤波器系统函数分子分母%直接形式转化为并联形式%检验衰减指标%数字滤波器的频域响应b1,a1=u_buttap(N,OC);b,a=butter(N,OC,'s');bz,az=impinvar(b,a,Fs)C,B,A=dir2par(bz,az);w0=Wp,Ws;Hx=freqz(bz,az,w0);H,w=freqz(bz,az);db,mag,p

33、ha,grd,w=freqz_m(bz,az);%求模拟滤波器振幅，相位等ha,x,t=impulse(b1,a1)%模拟滤波器时域响应hn,t=impz(bz,az)%数字滤波器时域响应w1=zeros(1,512-length(w),w;%矢量维度一致dbHx=-20*log10(abs(Hx)/max(abs(H);%检验衰减指标figure(1);plot(mag);title('模拟滤波器幅频');gridonfigure(2);plot(pha);title('模拟滤波器相频');gridonfigure(3);plot(w1/pi,20*log10

34、(abs(H);title('数字滤波器幅频,);xlabel('w/pi');ylabel('幅度响应曲线dB');gridonfigure(4);plot(w1/pi,angle(H);title('数字滤波器相频');xlabel('w/pi');ylabel('相位响应曲线');gridonfigure(5);plot(ha);title('冲击响应");gridonfigure(6);plot(t,hn);title('单位抽样响应,);xlabel('n'

35、;,'FontSize',12);ylabel('h(n)','fontsize',12);gridonfigure(7):plot(w/pi,db);hn=hn;zeros(128-length(hn),1);figure(8);y=20*log10(abs(fft(hn,1024);%时域响应的FFT变换求幅度响应x=(0:1023)/(512);%调整FFT生成的波形y1=y(x>=0&x<=1);plot(x(x>=0&x<=1),y1);xlabel('w/pi');ylabel（

36、'幅度响应曲线dB'）;gridonfigure(9) y3=y2(x>=0&x<=1);y2=angle(fft(hn,1024);%时域响应的FFT变换求相位响应x=(0:1023)/(512);%调整FFT生成的波形plot(x(x>=0&x<=1),y3);xlabel('w/pi');ylabel（'相位响应曲线'）;gridonfigure(10)plot(w1/pi,20*log10(abs(H),'*r',x(x>=0&x<=1),y1,'ob&#

37、39;);%将两种方法曲线画在图上比较legend（'直接求响应的曲线','FFT求响应的曲线'）;xlabel('w/pi');ylabel（'幅度响应曲线dB'）;figure(11)plot(w1/pi,angle(H),'.r',x(x>=0&x<=1),y3,'ob');%将两种方法曲线画在图上比较legend('直接求响应的曲线,'FFT求响应的曲线');xlabel('w/pi');ylabel('相位响应曲线');functionC,B,A=dir2par(b,a)M=length(b);N=length(a);r1,p