数字信号处理滤波器设计

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

2、R 滤波器设计过程FIR 滤波器的设计：1. 设计步骤：（ 1）写出理想低通滤波器的传输函数和单位脉冲响应。（ 2）根据题目已有的模拟滤波器频率进行归一化求得数字滤波器的频率，并对题目所给的通带波动和阻带波动进行转换，可求得通带截止频率，阻带截止频率，通带波动，阻带衰减，（ 3）选择用布莱克曼窗函数来设计FIR 滤波器，通过查阅资料可知道布莱克曼窗的过渡带宽为，最小衰减为，则可以通过 matlab 中 ceil 函数向上取整的方式以及阻带频率和通带频率之差所求的过渡带频率来确定滤波器的长度,在经过调试可以取得N=99 。当然，在设计中可以之间用blackman() 这个函数，但是在这里我们会利

3、用布莱克曼窗函数的公式构造布莱克曼窗：.（ 4）通过求得的矩形窗函数和单位脉冲响应进行时域上的乘积即可求得FIR 滤波器的时域响应，再利用DTFT 直接求出，并且可以画出幅度响应曲线（ dB）及相位响应曲线。（ 5）根据时域所求的函数做 FFT并由此画出幅度响应曲线（ dB）及相位响应曲线，对这两种方法的结果进行比较并验证设计结果是否满足设计指标。3. 设计过程：（ 1）求出，并做的 DTFT 从而直接求出的幅度曲线响应（ dB）和相位响应曲线。由设计步骤中可以知道以及自己通过公式所构造的布莱克曼窗函数即可求得。然后将经过 DTFT 即可求得，进而求得幅度曲线响应和相位响应曲线。波形如下：.图

4、一的波形图二 理想滤波器冲击响应hdn 的波形图三的波形的前 100 个值Columns 1 through 13-0.0000-0.0000-0.0000-0.00000.00000.00010.0000.-0.0001-0.0002-0.00010.00020.00040.0003Columns 14 through 26-0.0003-0.0008-0.00050.00050.00140.0010-0.0008-0.0023-0.00160.00120.00360.0026-0.0016Columns 27 through 39-0.0053-0.00410.00210.00770.00

5、62-0.0026-0.0109-0.00930.00310.01540.0139-0.0035-0.0220Columns 40 through 52-0.02110.00390.03270.0340-0.0043-0.0549-0.06480.00450.13910.27290.32880.27290.1391Columns 53 through 650.0045-0.0648-0.0549-0.00430.03400.03270.0039-0.0211-0.0220-0.00350.01390.01540.0031Columns 66 through 78-0.0093-0.0109-0

6、.00260.00620.00770.0021-0.0041-0.0053-0.00160.00260.00360.0012-0.0016Columns 79 through 91-0.0023-0.00080.00100.00140.0005-0.0005-0.0008-0.00030.00030.00040.0002-0.0001-0.0002Columns 92 through 99-0.00010.00000.00010.0000-0.0000-0.0000-0.0000-0.0000.图四的幅度响应曲线图四的相位响应曲线（ 2）根据时域所求的函数做 FFT并由此画出幅度响应曲线（dB

7、）及相位响应曲线。.图五的幅度响应曲线图六的相位响应曲线由题目可知滤波器的指标：通带波动为0.174dB，阻带衰减为 -70dB ，通带截止频率为 0.8549Hz，阻带截止频率为1.2110Hz。通过 DTFT直接求所得到的 FIR滤波器的指标是：通带波动为 0.001dB，阻带衰减为 -73.14 dB 。通过 FFT 求所得到的 FIR滤波器的指标是： 通带波动为 0.0014dB，阻带衰减为 -75 dB 。.图七的幅度响应曲线 ()由图六可以知道这个滤波器两个正负尖峰频率差值为：1.2612-0.80381=0.45739Hz设计的的 FIR滤波器的过渡带宽为：1.2110-0.85

8、49=0.3561Hz图八 DTFT和 FFT所求的幅度响应 (dB)曲线由该图可以看出 DTFT 所求幅度响应（ dB）的曲线和 FFT所求的幅度响.应（ dB）的曲线是基本一致的。图八 DTFT和 FFT所求的相位响应曲线综上所述： 1. 所设计的布莱克曼窗FIR滤波器在阻带衰减上接近了-70(dB)的指标，但是在通带波动上0.02(dB)略小于所要求的0.174(dB)的指标。2. FIR滤波器的过渡带宽小于两个正负尖峰频率差值，满足要求。 3. DTFT和 FFT两种方法所设计的滤波器结果及指标基本没有区别，但是DTFT的计算量会比 FFT相对大，而 FFT所求的幅度和相位响应的波形是

9、对称的，需要经过处理才能达到要求。 4. DTFT所算出来的指标会相对来说会更接近于题目中要求的指标。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）计算样本模拟滤波器的阶次可求得 N=17 。（ 3）发现 Rp3(dB)，由于巴特沃斯滤波器

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

11、器系统函数分子分母系数来检验衰减指标。3. 设计过程：.（ 1）将要求的各种指标代入巴特沃斯函数中以求得阶数和3（dB）衰减的截止频率，并由此确定模拟滤波器系统函数的分子分母系数，再利用impinvar 将模拟滤波器系统函数的分子分母系数转换成数字滤波器系统函数的分子分母系数。然后利用 freqz （）函数求出数字滤波器频域上的幅度响应曲线（dB）和相位响应曲线。图一的幅度响应曲线图二的相位响应曲线（ 2）将（ 1）中所求的数字滤波器系统函数的分子分母的系数利用imz 函数即可.求得系统的冲激响应即时域函数。图三的波形的前 100 个值Columns 1 through 13-0.0000-0

12、.00000.00000.00000.00040.00510.03090.10830.24340.36100.33520.1315-0.1050Columns 14 through 26-0.1764-0.05330.09240.0972-0.0112-0.0796-0.03560.04150.0494-0.0062-0.0414-0.01650.0239Columns 27 through 390.0250-0.0060-0.0226-0.00670.01440.0124-0.0050-0.0123-0.00220.00860.0060 -0.0037-0.0066Columns 40 th

13、rough 52-0.00040.00500.0028-0.0026-0.00350.00030.0029.0.0012-0.0017-0.00180.00040.00160.0005Columns 53 through 65-0.0010-0.00090.00040.00090.0002-0.0006-0.00040.00030.00050.0000-0.0004-0.00020.0002Columns 66 through 780.0003-0.0000-0.0002-0.00010.00010.0001-0.0000-0.0001-0.00000.00010.0001-0.0000-0.

14、0001Columns 79 through 91-0.00000.00000.0000-0.0000-0.000000000000Columns 92 through 1040000000000000（ 3）根据时域所求的函数 h(n)做 FFT并由此画出幅度响应曲线（ dB）及相位响应曲线。图四的幅度响应曲线.图五的相位响应曲线由题目可知滤波器的指标：通带波动为0.087dB，阻带衰减为 26dB，通带截止频率为 1.1398Hz，阻带截止频率为1.5387Hz。由于求利用 impz 函数求 h(n)只是数字滤波器的单位冲激响应逼近模拟滤波器的单位冲激响应的过程， 再进行 FFT变换所得到

15、的数字滤波器的指标误差相对较大，所以在这里只检验直接求频域响应所得到的数字滤波器的指标。利用书上检验衰减指标的代码并将参数代入可求得：通带波动为0.0634 dB，阻带衰减为 26dB，非常接近题目中要求的指标。图六 直接求和 FFT所求的幅度响应（ dB）曲线.图七 直接求和 FFT所求的相位响应曲线由以上两幅图可以知道直接求所得到的幅度响应曲线及相位响应曲线在w/pi=00.7 这个范围内的波形基本保持一致，但是在 w/pi 0.7 之后 FFT所得到的幅度响应的曲线的谐波较多，而且呈现的也不是线性相位。综上所述，直接求所设计的低通滤波器基本能够满足题目要求的指标，在比较了直接求和 FFT

16、设计的数字低通滤波器后， 认为在设计低通巴特沃斯滤波器上还是用直接求的方式较为精确。设计过程中遇到的问题和解决方法：1. FIR滤波器中在求窗长度的时候 11pi/(ws-wp) 得不到整数。解决方法：用 ceil 函数向上取整的方式得到整数。2. 在做 FFT变换的时候发现结果是呈对称的形式而且曲线呈现的波形不够完整。解决方法：首先将波形的范围延长一倍， 然后采用想要的 x 的范围去限制 matlab画出的 y 的波形，就可以得到想要的一般的波形。3. 如何克服矢量维度不一致的问题。.解决方法：先看两个矢量的维度是多少， 然后可以在那个维度少的矢量后边用补零函数进行调整使得两个函数的矢量一致

17、。心得体会：通过这次设计 FIR滤波器和 IIR 滤波器的经历，我真心觉得自己的matlab 水平真的有了一个质的飞跃， 我承认网络上有很多类似的资料，但是要真正达到题目的要求还是需要自己动脑去做的。感觉这个过程就像是建立模型和代码实现的过程，我也是按照书本上两个滤波器的模型来设计并且实现的，整个过程发现了不少的问题， 有些是关于课程的问题， 比如加窗类型的选择以及调试，在达到过渡带小于正负尖峰之间那个指标的过程真的是调了很久，最后才确立了布莱克曼窗，参数也是改了好多遍。有些是关于matlab 的基础问题，因为我接触matlab的时间不多，一开始比较吃力，但是我觉得经过了这个项目之后确实对ma

18、tlab有了进一步的认识， 其实挺希望这个项目能稍微增加难度然后放在课程一开始让学生研究，我觉得这样的效果会更好。程序源代码：FIR 滤波器close allf1=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_w=ws-wp;M=ceil(c/delt_w);N=1*M+1%窗长度n=0:M;%利用课本公式得到布莱克曼窗win=(0.42-0.5*

19、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',12);ylabel('ht','fontsize',12);grid on.figure(2)% h

20、,w=freqz(ht,1,1024);%离散时间傅里叶变换得到幅度% W=w/pi;H=20*log10(abs(h);% hold on% plot(W,H);k=1:1024;n=0:M;w=(pi/1024)*k;W=w/pi; h=ht*(exp(-j*pi/1024).(n'*k);%DTFT 变换plot(w/pi,20*log10(abs(h);xlabel('w/pi');ylabel('幅度响应曲线 dB');%直接求所得到的频率幅度响应grid onhold online(0,1,-70,-70,'color',&#

21、39;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',10);% 对应处打上显目的红点plot(W(dots),H(dots),'.r',

22、'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),'','',num2str(round(W(dots)*100)/100),.',

23、9;,num2str(round(H(dots)*100)/100),'dB','fontsize',15);%标示 ws 在频响图%的位置hold offfigure(3)plot(W,angle(h);axis(0,0.5,-4,4);xlabel('w/pi');ylabel('相位响应曲线 ');%直接求所得到的相位响应曲线grid onfigure(4)y=20*log10(abs(fft(ht,1024);%FFT求得的频域响应.x=(0:1023)/(512);y1 = y(x>=0&x<=1);

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

25、'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);y3 = y2(x>=0&x<=1);

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

27、B');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);%标出幅度曲线在正负尖峰的两个点，检测是否满足过渡带的要求hold onaxis(0,2,-0.5,1.5);pl

28、ot(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');yl

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

30、11025%采样频率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);%数字滤波器的频域响应

31、db,mag,pha,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(' 模拟滤波器幅频 ') ;grid onfigure(2);plot(pha);title(' 模拟滤波器相频 ') ;grid onfigure(3

32、);plot(w1/pi,20*log10(abs(H);title(' 数字滤波器幅频 ') ;xlabel('w/pi');ylabel(' 幅度响应曲线dB');grid onfigure(4);plot(w1/pi,angle(H);title(' 数字滤波器相频 ') ;xlabel('w/pi');ylabel(' 相位响应曲线 ');grid onfigure(5);plot(ha);title(' 冲击响应 '); grid onfigure(6);plot(t,hn

33、);title(' 单位抽样响应 ') ;xlabel('n','FontSize',12);ylabel('h(n)','fontsize',12);grid onfigure(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&g

34、t;=0&x<=1),y1);xlabel('w/pi');.ylabel('幅度响应曲线 dB');grid onfigure(9)y2=angle(fft(hn,1024);%时域响应的 FFT变换求相位响应x=(0:1023)/(512);%调整 FFT生成的波形y3 = y2(x>=0&x<=1);plot(x(x>=0&x<=1),y3);xlabel('w/pi');ylabel('相位响应曲线 ');grid onfigure(10)plot(w1/pi,20*l

35、og10(abs(H),'*r',x(x>=0&x<=1),y1,'ob');%将两种方法曲线画在图上比较legend(' 直接求响应的曲线 ','FFT求响应的曲线 ');xlabel('w/pi');ylabel('幅度响应曲线 dB');figure(11)plot(w1/pi,angle(H),'.r',x(x>=0&x<=1),y3,'ob');% 将两种方法曲线画在图上比较.legend(' 直接求响应的曲线

36、 ','FFT求响应的曲线 ');xlabel('w/pi');ylabel('相位响应曲线 ');function C,B,A = dir2par(b,a)M=length(b);N=length(a);r1,p1,C=residuez(b,a);p=cplxpair(p1,10000000*eps);I=cplxcomp(p1,p);r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);if K*2=N;for i=1:2:N-2Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);Brow,Arow=res