数字信号处理实验五FIR滤波器设计汇总.docx

1、数字信号处理实验五FIR滤波器设计汇总实验五 FIR滤波器设计一实验内容（1）认真复习 FIR 数字滤波器的基本概念，线性相位 FIR 滤波器的条件和特点、幅度函数特点、零点位置的基本特点与性质；窗函数设计法的基本概念与方法，各种窗函数的性能和设计步骤，线性相位 FIR 低通、高通、带通和带阻滤波器的设计方法，频率采样设计法的基本概念和线性相位的实现方法。 （2）掌握几种线性相位的特点，熟悉和掌握矩形窗、三角形窗、汉宁窗、海明窗、布莱克曼窗、 凯塞窗设计 IIR 数字滤波器的方法， 熟悉和掌握频率抽样设计法的线性相位的设计方法，并对各种线性相位的频率抽样法的设计给出调整和改进。（3）熟悉利用

2、MATLAB 进行各类 FIR 数字滤波器的设计方法。 二实验内容a. 设线性相位 FIR 滤波器单位抽样响应分别为 h(n )= -4,1, -1,- 2,5,6,5, -2, -1,1, -4h(n)= -4,1, -1,- 2,5,6,6,5,- 2, -1,1,- 4h(n)= -4,1,- 1, -2,5,0,- 5,2,1,- 1,4h(n)= -4,1, -1, -2,5,6,- 6, -5,2,1, -1,4分别求出滤波器的幅度频率响应 H()，系统函数 H(z)以及零极点分布，并绘制相应的波形和分布图。在matlab中新建函数amplres，代码如下:functionA,w,

3、type,tao=amplres(h)N=length(h);tao=(N-1)/2;L=floor(N-1)/2);n=1:L+1;w=0:500*2*pi/500;if all(abs(h(n)-h(N-n+1)1e-10)A=2*h(n)*cos(N+1)/2-n)*w)-mod(N,2)*h(L+1);type=2-mod(N,2);elseif all(abs(h(n)+h(N-n+1)a1type = 1rz = -0.9807 + 0.1956i -0.9807 - 0.1956i -0.5578 + 0.8300i -0.5578 - 0.8300i 0.4052 + 1.23

4、74i 0.4052 - 1.2374i 1.2169 + 0.0000i 0.8218 + 0.0000i 0.2390 + 0.7299i 0.2390 - 0.7299ians = -0.9807 + 0.1956i -0.9807 - 0.1956i -0.5578 + 0.8300i -0.5578 - 0.8300i 0.2390 + 0.7299i 0.2390 - 0.7299i 0.8218 + 0.0000i 1.2169 + 0.0000i对第二个单位抽样响应，在matlab中新建函数a2，代码如下：h2=-4,1,-1,-2,5,6,6,5,-2,-1,1,-4;M=l

5、ength(h2);n=0:M-1;A,w,type,tao=amplres(h2);typesubplot(2,1,1),stem(n,h2);title(冲激响应h2);ylabel(h(n);xlabel(n);subplot(2,1,2),plot(w/pi,A);ylabel(A);xlabel(pi);title(幅频响应);figurerz=roots(h2)for i=1:8r(i)=1/rz(i);endrzplane(h2,1);title(h2零极点图);生成结果如下：a2type = 2rz = 1.3120 + 0.0000i -1.0000 + 0.0000i -0

6、.8868 + 0.4622i -0.8868 - 0.4622i 0.4778 + 1.1851i 0.4778 - 1.1851i -0.2957 + 0.9553i -0.2957 - 0.9553i 0.7622 + 0.0000i 0.2926 + 0.7258i 0.2926 - 0.7258ians = 0.7622 + 0.0000i -1.0000 + 0.0000i -0.8868 + 0.4622i -0.8868 - 0.4622i 0.2926 + 0.7258i 0.2926 - 0.7258i -0.2957 + 0.9553i -0.2957 - 0.9553i

7、对第三个单位抽样响应，在matlab中新建函数a3，代码如下：h3=-4,1,-1,-2,5,0,-5,2,1,-1,4;M=length(h3);n=0:M-1;A,w,type,tao=amplres(h3);typesubplot(2,1,1),stem(n,h3);title(冲激响应h3);ylabel(h(n);xlabel(n);subplot(2,1,2),plot(w/pi,A);ylabel(A);xlabel(pi);title(幅频响应);figurerz=roots(h3)for i=1:8r(i)=1/rz(i);endrzplane(h3,1);title(h3零

8、极点图);生成结果如下： a3type = 3rz = -1.0000 + 0.0000i -0.8732 + 0.4874i -0.8732 - 0.4874i 0.1010 + 1.2041i 0.1010 - 1.2041i 1.0000 + 0.0000i 0.8280 + 0.5607i 0.8280 - 0.5607i 0.0692 + 0.8247i 0.0692 - 0.8247ians = -1.0000 + 0.0000i -0.8732 + 0.4874i -0.8732 - 0.4874i 0.0692 + 0.8247i 0.0692 - 0.8247i 1.0000

9、 + 0.0000i 0.8280 + 0.5607i 0.8280 - 0.5607i对第四个单位抽样响应，在matlab中新建函数a4，代码如下：h4=-4,1,-1,-2,5,6,-6,-5,2,1,-1,4;M=length(h4);n=0:M-1;A,w,type,tao=amplres(h4);typesubplot(2,1,1),stem(n,h4);title(冲激响应h4);ylabel(h(n);xlabel(n);subplot(2,1,2),plot(w/pi,A);ylabel(A);xlabel(pi);title(幅频响应);figurerz=roots(h4)f

10、or i=1:8r(i)=1/rz(i);endrzplane(h4,1);title(h4零极点图);生成结果如下： a4type = 4rz = 0.2631 + 1.3394i 0.2631 - 1.3394i -0.9505 + 0.3108i -0.9505 - 0.3108i -0.7309 + 0.6825i -0.7309 - 0.6825i 0.9021 + 0.4315i 0.9021 - 0.4315i 1.0000 + 0.0000i 0.1412 + 0.7189i 0.1412 - 0.7189ians = 0.1412 + 0.7189i 0.1412 - 0.7

11、189i -0.9505 + 0.3108i -0.9505 - 0.3108i -0.7309 + 0.6825i -0.7309 - 0.6825i 0.9021 + 0.4315i 0.9021 - 0.4315ib. 设计 FIR 数字低通滤波器，技术指标为：p=0.2，st=0.3，1=0.25dB，2=50dB。 (1) 通过技术指标，选择一种窗函数进行设计； (2) 求滤波器的单位抽样响应、频率响应，并绘制波形。 (3) 选择凯塞窗函数设计该滤波器，并绘制相应的波形图。在matlab中新建函数ideal_lp，代码如下：function hd=ideal_lp(wc,M);% I

12、deal LowPass filter computation% -% hd = ideal_lp(wc,M);% hd = ideal impulse response between 0 to M-1% wc = cutoff frequency in radians% M = length of the ideal filter%alpha = (M-1)/2;n = 0:1:(M-1);m = n - alpha +eps; % add smallest number to avoi divided by zerohd = sin(wc*m)./(pi*m);end在matlab中新建

13、函数freqz_m，代码如下：function db,mag,pha,grd,w = freqz_m(b,a);H,w = freqz(b,a,1000,whole); H = (H(1:1:501); w = (w(1:1:501); mag = abs(H); db = 20*log10(mag+eps)/max(mag); pha = angle(H); grd = grpdelay(b,a,w);在matlab中新建函数b1，代码如下：%用Hamming窗函数设计FIR数字滤波器wp=0.2*pi;ws=0.3*piN=61n=0:1:N-1wc=(ws+wp)/2;%理想低通滤波器h

14、d=ideal_lp(wc,N);%理想低通的冲激响应w_ham=(hamming(N);h=hd.*w_ham;%FIR滤波器冲激响应db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1)%实际的通带衰减As=-round(max(db(ws/delta_w+1:1:501)%实际的最小阻带衰减subplot(221);stem(n,hd);title(理想冲激响应)axis(0 N-1 -0.1 0.3);xlabel(n);ylabel(hd(n)subplot(222);stem(n,w

15、_ham);title(hamming窗)axis(0 N-1 0 1.1);xlabel(n);ylabel(w(n);subplot(223);stem(n,h);title(实际冲激响应)axis(0 N-1 -0.1 0.3);xlabel(n);ylabel(h(n)subplot(224);plot(w/pi,db);axis(0 0.8 -100 0);xlabel(以PI为单位的频率);ylabel(对数幅度/db);实验结果如下：在matlab中新建函数b2，代码如下：%用Kaiser窗函数设计FIR数字滤波器wp=0.2*pi;ws=0.3*pi;As=50tr_width

16、=ws-wpN=ceil(As-7.95)/(14.36*tr_width/(2*pi)+1)+1n=0:1:N-1beta=0.1102*(As-8.7)wc=(wp+ws)/2%理想低通的截止频率hd=ideal_lp(wc,N)w_kai=(kaiser(N,beta)h=hd.*w_kaidb,mag,pha,grd,w=freqz_m(h,1)delta_w=2*pi/1000Rp=-(min(db(1:1:wp/delta_w+1)%实际的通带衰减As=-round(max(db(ws/delta_w+1:1:501)%实际的最小阻带衰减subplot(211);plot(w/pi

17、,db);title(凯森窗幅度响应(dB);gridaxis(0 0.5 -100 0)ylabel(对数幅度/db);xlabel(以/pi为单位的频率)subplot(212);plot(w/pi,pha);title(相位响应);gridaxis(0 0.5 -4 4)ylabel(相位);xlabel(以/pi为单位的频率)实验结果如下：c. 设计 FIR 数字带通滤波器，技术指标为： 下阻带边缘：st1=0.2，s1=60dB，下通带边缘：p1=0.35，p1=1dB； 上通带边缘：p2=0.65，p1=1dB，上阻带边缘：st2=0.8，s2=60dB； (1) 通过技术指标，选

18、择一种窗函数进行设计； (2) 求滤波器的单位抽样响应、频率响应，并绘制波形。在matlab中新建函数c，代码如下：ws1=0.2*pi; wp1=0.35*pi;ws2=0.8*pi;wp2=0.65*pi;Ap=60;Rp=1;tr_width=min(wp1-ws1),(ws2-wp2);M=ceil(11*pi/tr_width);n=0:1:M-1;wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);w_bla=(blackman(M);h=hd.*w_bla;H,W=freqz(h,1);subplo

19、t(2,2,1);stem(n,hd);title(理想脉冲抽样);subplot(2,2,2);stem(n,w_bla);title(布莱克曼窗);subplot(2,2,3);stem(n,h);title(实际脉冲抽样);subplot(2,2,4);plot(W/pi,20*log10(abs(H);title(幅度响应（db）);生成结果如下：d. 设计 FIR 数字带通滤波器，技术指标为： 下阻带边缘：st1=0.2，s1=60dB，下通带边缘：p1=0.4，p1=1dB； 上通带边缘：p2=0.6，p1=1dB，上阻带边缘：st2=0.8，s2=60dB； (1) 通过技术指标

20、，选择一种窗函数进行设计； (2) 求滤波器的单位抽样响应、频率响应，并绘制波形。在matlab中新建函数d，代码如下：ws1=0.2*pi; wp1=0.4*pi;ws2=0.8*pi;wp2=0.6*pi;Ap=60;Rp=1;tr_width=min(wp1-ws1),(ws2-wp2);M=ceil(11*pi/tr_width);n=0:1:M-1;wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);w_bla=( hamming(M);h=hd.*w_bla;H,W=freqz(h,1);subplo

21、t(2,2,1);stem(n,hd);title(理想脉冲抽样);subplot(2,2,2);stem(n,w_bla);title(海明窗);subplot(2,2,3);stem(n,h);title(实际脉冲抽样);subplot(2,2,4);plot(W/pi,20*log10(abs(H);title(幅度响应（db）);生成结果如下：e. 设计 FIR 数字带通滤波器，技术指标为： 下阻带边缘：st1=0.2，s1=20dB，下通带边缘：p1=0.4，p1=1dB； 上通带边缘：p2=0.6，p1=1dB，上阻带边缘：st2=0.8，s2=20dB； (1) 通过技术指标，选

22、择一种窗函数进行设计； (2) 求滤波器的单位抽样响应、频率响应，并绘制波形。在matlab中新建函数e，代码如下： ws1=0.2*pi;wp1=0.4*pi;ws2=0.8*pi;wp2=0.6*pi;Ap=20;Rp=1;tr_width=min(wp1-ws1),(ws2-wp2);M=ceil(11*pi/tr_width);n=0:1:M-1;wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);w_bla=( boxcar(M);h=hd.*w_bla;H,W=freqz(h,1);subplot(2

23、,2,1);stem(n,hd);title(理想脉冲抽样);subplot(2,2,2);stem(n,w_bla);title(矩形窗);subplot(2,2,3);stem(n,h);title(实际脉冲抽样);subplot(2,2,4);plot(W/pi,20*log10(abs(H);title(幅度响应（db）);生成结果如下：f. 设计 FIR 数字高通滤波器，技术指标为：通带截止频率为p=15/27，阻带截止频率为st=11/27，通带最大衰减为1=2.5dB，阻带最小衰减为2=55dB。 (1) 通过技术指标，选择一种窗函数进行设计； (2) 求滤波器的单位抽样响应、频

24、率响应，并绘制波形。系统一：在matlab中新建函数f，代码如下： As=55;ws=11*pi/27;wp=15*pi/27; tr_width=wp-ws;%计算过渡带 M=ceil(As-7.95)*2*pi/(14.36*tr_width)+1)+1; %按凯泽窗计算滤波器的长度disp(滤波器的长度,num2str(M);beta=0.1102*(As-8.7); %计算凯泽窗的beta值n=0:1:M-1;disp(线性相位滤波器,num2str(beta);w_kai=(kaiser(M,beta);%求凯泽窗函数wc=(ws+wp)/2;hd=ideal_lp(pi,M)-id

25、eal_lp(wc,M); %求理想脉冲响应h=hd.*w_kai; db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000;Rp=-(min(db(wp/delta_w+1:1:501);disp(实际通带波动为,num2str(Rp);As=-round(max(db(1:1:ws/delta_w+1);disp(最小阻带衰减为,num2str(As);subplot(2,2,1);stem(n,hd);title(理想脉冲响应);axis(0 M-1 -0.4 0.8);ylabel(hd(n);subplot(2,2,2);stem(n,w_ka

26、i);title(凯泽窗); axis(0 M-1 0 1.1);ylabel(wd(n);subplot(2,2,3);stem(n,h);title(实际脉冲响应);axis(0 M-1 -0.4 0.8);xlabel(n);ylabel(h(n);subplot(2,2,4);plot(w/pi,db);title(幅度响应/dB);axis(0 1 -100 10);grid;xlabel(以pi为单位的频率);ylabel(分贝数/dB); 生成结果如下：滤波器的长度47线性相位滤波器5.1023实际通带波动为0.025017最小阻带衰减为?55滤波器的长度47线性相位滤波器5.1

27、023实际通带波动为0.025017最小阻带衰减为55g. 设计 FIR 数字高通滤波器，技术指标为：通带截止频率为p=0.6，阻带截止频率为st=0.4，通带最大衰减为1=0.25dB，阻带最小衰减为2=40dB。 (1) 通过技术指标，选择一种窗函数进行设计； (2) 求滤波器的单位抽样响应、频率响应，并绘制波形。 在matlab中新建函数c，代码如下：function db, mag, pha, w=freqz_m2(b,a);%Modified version of freqz subroutineH,w=freqz(b,a,1000,whole);H=(H(1:501); w=(w(1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);在matlab中新建函数g，代码如下：Wp=0.6*pi;Ws=0.4*pi;tr_width=Wp-Ws;M=ceil(6.2*pi/tr_width);n=0:1:M-1;Wc=(Ws+Wp)/