数字信号处理实验五FIR滤波器设计

实验五FIR滤波器设计一．实验内容（1）认真复习FIR 数字滤波器的基本概念，线性相位FIR 滤波器的条件和特点、幅度函数特点、零点位置的基本特点与性质；窗函数设计法的基本概念与方法，各种窗函数的性能和设计步骤，线性相位FIR 低通、高通、带通和带阻滤波器的设计方法，频率采样设计法的基本概念和线性相位的实现方法。

（2）掌握几种线性相位的特点，熟悉和掌握矩形窗、三角形窗、汉宁窗、海明窗、布莱克曼窗、凯塞窗设计IIR 数字滤波器的方法，熟悉和掌握频率抽样设计法的线性相位的设计方法，并对各种线性相位的频率抽样法的设计给出调整和改进。

（3）熟悉利用MATLAB 进行各类FIR 数字滤波器的设计方法。

二．实验内容a. 设线性相位FIR 滤波器单位抽样响应分别为h(n )={ -4,1, -1,- 2,5,6,5, -2, -1,1, -4}h(n)={ -4,1, -1,- 2,5,6,6,5,- 2, -1,1,- 4}h(n)={ -4,1,- 1, -2,5,0,- 5,2,1,- 1,4}h(n)={ -4,1, -1, -2,5,6,- 6, -5,2,1, -1,4}分别求出滤波器的幅度频率响应H(ω)，系统函数H(z)以及零极点分布，并绘制相应的波形和分布图。

在matlab中新建函数amplres，代码如下:function[A,w,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))<1e-10)&(h(L+1)*mod(N,2)==0)A=2*h(n)*sin(((N+1)/2-n)'*w);type=4-mod(N,2);else error('错误，这不是线性相位滤波器!')end对第一个单位抽样响应，在matlab中新建函数a1，代码如下：h1=[-4,1,-1,-2,5,6,5,-2,-1,1,-4];M=length(h1);n=0:M-1;[A,w,type,tao]=amplres(h1);typesubplot(2,1,1),stem(n,h1);title('冲激响应h1');ylabel('h(n)');xlabel('n');subplot(2,1,2),plot(w/pi,A);ylabel('A');xlabel('\pi');title('·幅频响应');figurerz=roots(h1)for i=1:8r(i)=1/rz(i);endr'zplane(h1,1);title('h1零极点图');生成结果如下：>>a1type =1rz =-0.9807 + 0.1956i-0.9807 - 0.1956i-0.5578 + 0.8300i-0.5578 - 0.8300i 0.4052 + 1.2374i0.4052 - 1.2374i1.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.7299i0.8218 + 0.0000i1.2169 + 0.0000i对第二个单位抽样响应，在matlab中新建函数a2，代码如下：h2=[-4,1,-1,-2,5,6,6,5,-2,-1,1,-4];M=length(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);endr'zplane(h2,1);title('·h2零极点图');生成结果如下：>>a2type =2rz =1.3120 + 0.0000i -1.0000 + 0.0000i -0.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对第三个单位抽样响应，在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);endr'zplane(h3,1);title('h3零极点图');生成结果如下：>> a3type =3rz =-1.0000 + 0.0000i-0.8732 + 0.4874i-0.8732 - 0.4874i0.1010 + 1.2041i0.1010 - 1.2041i1.0000 + 0.0000i0.8280 + 0.5607i0.8280 - 0.5607i0.0692 + 0.8247i0.0692 - 0.8247ians =-1.0000 + 0.0000i -0.8732 + 0.4874i -0.8732 - 0.4874i 0.0692 + 0.8247i0.0692 - 0.8247i1.0000 + 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)for i=1:8r(i)=1/rz(i);endr'zplane(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.4315i0.9021 - 0.4315i1.0000 + 0.0000i 0.1412 + 0.7189i 0.1412 - 0.7189ians =0.1412 + 0.7189i 0.1412 - 0.7189i -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);% Ideal 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 zero hd = sin(wc*m)./(pi*m);end在matlab中新建函数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-1]wc=(ws+wp)/2;%理想低通滤波器hd=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_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=ws-wpN=ceil((As-7.95)/(14.36*tr_width/(2*pi))+1)+1n=[0:1:N-1]beta=0.1102*(As-8.7)wc=(wp+ws)/2%理想低通的截止频率hd=ideal_lp(wc,N)w_kai=(kaiser(N,beta))'h=hd.*w_kai[db,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,db);title('凯森窗幅度响应(dB)');grid axis([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) 通过技术指标，选择一种窗函数进行设计；(2) 求滤波器的单位抽样响应、频率响应，并绘制波形。