实验11 用MATLAB设计FIR数字滤波器

课程名称：数字信号处理 实验成绩： 指导教师：

实 验 报 告

院系： 信息工程学院 班级： 电信二班 学号： 姓名： 日期：

实验11 用MATLAB设计FIR数字滤波器

一、实验目的：

1、加深对窗函数法设计FIR数字滤波器的基本原理的理解。

2、学习用MATLAB语言的窗函数法编写设计FIR数字滤波器的程序。 3、了解MATLAB语言有关窗函数法设计FIR数字滤波器的常用函数用法。

二、实验内容及步骤

2、选择合适的窗函数设计FIR数字低通滤波器，要求： wp=0.2π，Rp=0.05dB； ws=0.3π，As=40dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。 分析：根据设计指标要求，并查表11-1，选择汉宁窗。程序清单如下： function hd=ideal_lp(wc,N)

wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp; tao=(N-1)/2; n=[0:(N-1)];

m=n-tao+eps; hd=sin(wc*m)./(pi*m);

function[db,mag,pha,grd,w]=freqz_m(b,a); [H,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); grd=grpdelay(b,a,w);

wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp; wc=(ws+wp)/2;

第 1 页 共 9 页

N0=ceil(6.6*pi/deltaw); N=N0+mod(N0+1,2) windows=(hanning(N))'; hd=ideal_lp(wc,N); b=hd.*windows;

[db,mag,pha,grd,w]=freqz_m(b,1); n=0:N-1;dw=2*pi/1000;

Rp=-(min(db(1:wp/dw+1))) %检验通带波动

As=-round(max(db(ws/dw+1:501))) %检验最小阻带衰减 subplot(2,2,1);stem(n,b);

axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应'); xlabel('n');ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis([0,N,0,1.1]);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db);

axis([0,1,-80,10]);title('幅度频率响应');

xlabel('频率（单位：\\pi）');ylabel('H(e^{j\\omega})'); set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]); set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid subplot(2,2,4);plot(w/pi,pha); axis([0,1,-4,4]);title('相位频率响应');

xlabel('频率（单位：\\pi）');ylabel('\\phi(\\omega)'); set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]); set(gca,'YTickMode','manual','YTick',[-3.1416,0,3.1416,4]);grid N =67 Rp =0.0706 As = 44

第 2 页 共 9 页

实际脉冲响应10.20.10020406000窗函数特性wd(n)h(n)0.5204060n幅度频率响应0-343.1416n相位频率响应H(ej)-20()0-50-3.141600.20.3频率（单位：）100.20.3频率（单位：）1

3、用凯塞窗设计一个FIR数字高通滤波器，要求： wp=0.3π，Rp=0.1dB；ws=0.2π，As=50dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。 程序清单如下：

function hd=ideal_lp(wc,N) tao=(N-1)/2; n=[0:(N-1)]; m=n-tao+eps; hd=sin(wc*m)./(pi*m);

function[db,mag,pha,grd,w]=freqz_m(b,a); [H,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); grd=grpdelay(b,a,w);

第 3 页 共 9 页

wp=0.3*pi;ws=0.2*pi; deltaw=wp-ws; N0=ceil(6.6*pi/deltaw); N=N0+mod(N0+1,2)

windows=(kaiser(N,7.865))'; wc=(ws+wp)/2;

hd=ideal_lp(pi,N)-ideal_lp(wc,N); b=hd.*windows;

[db,mag,pha,grd,w]=freqz_m(b,1); n=0:N-1;dw=2*pi/1000;

Rp=-(min(db(wp/dw+1:501))) %检验通带波动 As=-round(max(db(1:ws/dw+1))) %检验最小阻带衰减 subplot(2,2,1);stem(n,b);

axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应'); xlabel('n');ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis([0,N,0,1.1]);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db);

axis([0,1,-100,2]);title('幅度频率响应');

xlabel('频率（单位：\\pi）');ylabel('H(e^{j\\omega})'); set(gca,'XTickMode','manual','XTick',[0,ws/pi,wp/pi,1]); set(gca,'YTickMode','manual','YTick',[,-40-20,-3,0]);grid subplot(2,2,4);plot(w/pi,pha); axis([0,1,-4,4]);title('相位频率响应');

xlabel('频率（单位：\\pi）');ylabel('\\phi(\\omega)'); set(gca,'XTickMode','manual','XTick',[0,ws/pi,wp/pi,1]); set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid N =67 Rp = 0.2321

第 4 页 共 9 页

As =32

实际脉冲响应0.80.61窗函数特性wd(n)h(n)0.40.20-0.20204060n幅度频率响应0-30.500204060n相位频率响应3.1416H(ej)()0-60-3.141600.20.3频率（单位：）100.20.3频率（单位：）1

4、选择合适的窗函数设计一个FIR数字带通滤波器，要求：fp1=3.5kHz，fp2=6.5kHz，Rp=0.05dB；fs1=2.5kHz，fs2=7.5kHz，As=60dB。滤波器采样频率Fs=20kHz。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

分析：根据设计指标应选择布莱克曼窗。程序清单如下： function[db,mag,pha,grd,w]=freqz_m(b,a); [H,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); grd=grpdelay(b,a,w); fp1=3.5;fp2=6.5; fs1=2.5;fs2=7.5; Fs=20;

ws1=fs1/(Fs/2)*pi;ws2=fs2/(Fs/2)*pi;

第 5 页 共 9 页

wp1=fp1/(Fs/2)*pi;wp2=fp2/(Fs/2)*pi; deltaw=wp1-ws1; N0=ceil(11*pi/deltaw);

N=N0+mod(N0+1,2) %为实现FIR类型1偶对称滤波器，应确保N为奇数 windows=blackman(N);

wc1=(ws1+wp1)/2/pi;wc2=(ws2+wp2)/2/pi; b=fir1(N-1,[wc1,wc2],windows); [db,mag,pha,grd,w]=freqz_m(b,1); n=0:N-1;dw=2*pi/1000;

Rp=-(min(db(wp1/dw+1:wp2/dw+1))) %检验通带波动 ws0=[1:ws1/dw+1,ws2/dw+1:501];%建立阻带频率样点数组 As=-round(max(db(ws0))) %检验最小阻带衰减 subplot(2,2,1);stem(n,b);

axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应'); xlabel('n');ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis([0,N,0,1.1]);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db);

axis([0,1,-150,10]);title('幅度频率响应');

xlabel('频率（单位：\\pi）');ylabel('H(e^{j\\omega})'); set(gca,'XTickMode','manual','XTick',[0,fs1,fp1,fp2,fs2,500]); set(gca,'YTickMode','manual','YTick',[-150,-40,-3,0]);grid subplot(2,2,4);plot(w/pi,pha); axis([0,1,-4,4]);title('相位频率响应');

xlabel('频率（单位：\\pi）');ylabel('\\phi(\\omega)');

set(gca,'XTickMode','manual','XTick',[0,fs1,fp1,fp2,fs2,500]); set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid N =111 Rp =0.0034 As =74

第 6 页 共 9 页

实际脉冲响应0.40.21窗函数特性wd(n)050100n幅度频率响应h(n)0-0.20.50050100n相位频率响应0-3-403.1416H(ej)()0-3.1416-1500频率（单位：）0频率（单位：）

5、选择合适的窗函数设计一个FIR数字带阻滤波器，要求：fp1=1kHz，fp2=4.5kHz，

Rp=0.1dB；fs1=2kHz，fs2=3.5kHz，As=40dB。滤波器采样频率Fs=10kHz。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

分析：根据设计指标选择汉宁窗。程序清单如下： function hd=ideal_lp(wc,N) tao=(N-1)/2; n=[0:(N-1)]; m=n-tao+eps; hd=sin(wc*m)./(pi*m);

function[db,mag,pha,grd,w]=freqz_m(b,a); [H,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);

第 7 页 共 9 页

grd=grpdelay(b,a,w);

fp1=1;fp2=4.5; fs1=2;fs2=3.5; Fs=10;

ws1=fs1/(Fs/2)*pi;ws2=fs2/(Fs/2)*pi; wp1=fp1/(Fs/2)*pi;wp2=fp2/(Fs/2)*pi; wp=[wp1,wp2];ws=[ws1,ws2]; deltaw=ws1-wp1; N0=ceil(6.2*pi/deltaw);

N=N0+mod(N0+1,2) %为实现FIR类型1偶对称滤波器，应确保N为奇数 windows=(hanning(N))';

wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;

hd=ideal_lp(wc1,N)+ideal_lp(pi,N)-ideal_lp(wc2,N);%建立理想带阻 b=hd.*windows;

[db,mag,pha,grd,w]=freqz_m(b,1); n=0:N-1;dw=2*pi/1000;

wp0=[1:wp1/dw+1,wp2/dw+1:501];%建立通带频率样点数组 As=-round(max(db(ws1/dw+1:ws2/dw+1))) %检验最小阻带衰减 Rp=-(min(db(wp0))) %检验通带波动 subplot(2,2,1);stem(n,b);

axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应'); xlabel('n');ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis([0,N,0,1.1]);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db);

axis([0,1,-150,10]);title('幅度频率响应');

xlabel('频率（单位：\\pi）');ylabel('H(e^{j\\omega})');

set(gca,'XTickMode','manual','XTick',[0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1]); set(gca,'YTickMode','manual','YTick',[-150,-40,-3,0]);grid subplot(2,2,4);plot(w/pi,pha);

第 8 页 共 9 页

axis([0,1,-4,4]);title('相位频率响应');

xlabel('频率（单位：\\pi）');ylabel('\\phi(\\omega)');

set(gca,'XTickMode','manual','XTick',[0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1]); set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid N =31 As =44 Rp =0.0938

实际脉冲响应10.4窗函数特性wd(n)0102030h(n)0.200.500102030n幅度频率响应0-3-403.1416n相位频率响应H(ej)()0-3.1416-15000.20.40.70.91频率（单位：）00.20.40.70.91频率（单位：）

第 9 页 共 9 页