用MATLAB设计FIR数字的滤波器某实验八

实验八：用MATLAB 设计FIR 数字滤波器
1、 选择合适的窗函数设计FIR 数字低通滤波器，要求： w p =0.2π，R p =0.05dB ； w s =0.3π，A s =40dB 。

描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

程序清单如下：
wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp; N0=ceil(6.6*pi/deltaw);
N=N0+mod(N0+1,2) %为实现FIR 类型1偶对称滤波器，应确保N 为奇数 windows=(hamming(N))'; wc=(ws+wp)/2;
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
函数 ideal_lp 调用部分如下
function hd=ideal_lp(wc,N) %点0到N-1之间的理想脉冲响应 %wc=截止频率（弧度）
课程名称：数字信号处理 实验成绩： 指导教师： 实 验 报 告
院系： 信息工程学院 班级： 学号： 姓名：
日期： 2011. 11.26
%N=理想滤波器的长度 tao=(N-1)/2; n=[0:(N-1)];
m=n-tao+eps; %加一个小数以避免0作除数
hd=sin(wc*m)./(pi*m);
函数 freqz_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);
程序运行结果如下： N = 67
Rp =
0.0394
As = 52
20
40
60
-0.05
00.050.10.150.20.25实际脉冲响应
n
h (n )
20
40
60
00.20.40.60.8
1窗函数特性
n
w d (n
)
0.20.3
1
-50
-20
-30幅度频率响应
频率（单位：π）
H (e j ω)
0.20.3
1
-3.1416
3.1416
4
相位频率响应
频率（单位：π）
φ(ω)
2、 用凯塞窗设计一个FIR 数字高通滤波器，要求： w p =0.3π，R p =0.1dB ；w s =0.2π，A s =50dB 。

描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

程序清单如下：
N=61;As=50;
wp=0.3*pi;ws=0.2*pi;
beta=0.112*(As-8.7)
windows=kaiser(N,beta);
wc=(ws+wp)/2/pi
b=fir1(N-1,wc,'high',windows);
[db,mag,pha,grd,w]=freqz_m(b,1);
n=0:N-1;dw=2*pi/1000;
As=-round(max(db(1:ws/dw+1))) %检验最小阻带衰减
Rp=-(min(db(wp/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,ws/pi,wp/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,ws/pi,wp/pi,1]);
set(gca,'YTickMode','manual','YTick',[-3.1416,0,3.1416,4]);grid
函数 freqz_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);
程序运行结果如下：
beta =
4.6256
wc =
0.2500
As = 51
Rp =
0.0460
20
40
60
-0.2
00.20.40.6
0.8实际脉冲响应
n
h (n )
0204060
0.5
1
窗函数特性
n
w d (n )
0.2
0.3
1
-50
-20
-3
0幅度频率响应
频率（单位：π）
H (e j ω)
0.20.3
1
-3.1416
3.1416
4相位频率响应
频率（单位：π）
φ(ω)
3、 选择合适的窗函数设计一个FIR 数字带通滤波器，要求：f p1=3.5kHz ，f p2=6.5kHz ，R p =0.05dB ；f s1=2.5kHz ，f s2=7.5kHz ，A s =60dB 。

滤波器采样频率Fs=20kHz 。

描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

程序清单如下：
fp1=3500;fp2=6500; fs1=2500;fs2=7500; Fs=20000;
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=wp1-ws1;
N0=ceil(11*pi/deltaw);
N=N0+mod(N0+1,2) %为实现FIR 类型1偶对称滤波器，应确保N 为奇数 windows=(blackman(N))';
wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2; hd=ideal_lp(wc2,N)-ideal_lp(wc1,N); b=hd.*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,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1]); set(gca,'YTickMode','manual','YTick',[-100,-65,-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,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1]);
set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid
函数 ideal_lp 调用部分如下
function hd=ideal_lp(wc,N) %点0到N-1之间的理想脉冲响应
%wc=截止频率（弧度）
%N=理想滤波器的长度
tao=(N-1)/2;
n=[0:(N-1)];
m=n-tao+eps; %加一个小数以避免0作除数
hd=sin(wc*m)./(pi*m);
函数 freqz_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);
程序运行结果如下：
N =
111
Rp =
0.0034
As =
74
20
40
60
80
100
-0.3
-0.2-0.100.1
0.2
0.30.4实际脉冲响应
n
h (n )
02040
6080100
0.20.40.60.8
1
窗函数特性
n
w d (n )
0.250.350.650.75
1
-100
-65-20
-30幅度频率响应
频率（单位：π）
H (e j ω)
0.250.350.650.75
1
-3.1416
3.1416
相位频率响应
频率（单位：π）
φ(ω)
4、 选择合适的窗函数设计一个FIR 数字带阻滤波器，要求：f p1=1kHz ，f p2=4.5kHz ，R p =0.1dB ；f s1=2kHz ，f s2=3.5kHz ，A s =40dB 。

滤波器采样频率Fs=10kHz 。

描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。

程序清单如下：
fp1=1000;fp2=4500; fs1=2000;fs2=3500; Fs=10000;
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);
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
函数 ideal_lp 调用部分如下
function hd=ideal_lp(wc,N) %点0到N-1之间的理想脉冲响应
%wc=截止频率（弧度）
%N=理想滤波器的长度
tao=(N-1)/2;
n=[0:(N-1)];
m=n-tao+eps; %加一个小数以避免0作除数
hd=sin(wc*m)./(pi*m);
函数 freqz_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);
程序运行结果如下：
N =
31
As =
44
Rp =
0.0938
10
20
30
00.10.20.30.4
0.5实际脉冲响应
n
h (n )
0102030
0.2
0.40.60.81窗函数特性
n
w d (n
)
0.20.40.70.91
-150
-40
-30幅度频率响应
频率（单位：π）
H (e j ω
)
0.20.40.70.91
-3.1416
3.1416
相位频率响应
频率（单位：π）
φ(ω)
5、思考题：列写用窗函数法设计FIR 数字滤波器的基本方法，并写出不同方法设计低通、高通、带通、带阻滤波器的主要程序语句。

答：用窗函数法设计FIR 数字滤波器的基本方法如下：
（1）根据过渡带和阻带衰减指标选择窗函数的类型，估算滤波器的阶数N 。

（2）由数字滤波器的理想频率响应H(e j ω
)求出其单位冲激响应h d (n)。

（3）计算数字滤波器的单位冲激响应h(n)=w(n)h d (n)。

（4）检查设计的滤波器是否满足技术指标。

主要程序语句：
低通滤波器：①hd=ideal_lp(wc,N)
②b=fir1(N-1,wc,windows)
高通滤波器：①hd=ideal_lp(pi,N)-ideal_lp(wc,N) ②b=fir1(N-1,wc,'high',windows)
带通滤波器：①hd=ideal_lp(wc2,N)-ideal_lp(wc1,N) ②b=fir1(N-1,[wc1,wc2],windows)
带阻滤波器：①hd=ideal_lp(wc1,N)+ideal_lp(pi,N)-ideal_lp(wc2,N)
②b=fir1(N-1,[wc1,wc2],'stop',windows)
6、实验总结：
通过本次实验，我知道了采用第一类线性相位滤波器；在由过渡带宽度确定滤波器阶数N时，若N取奇数时，则可以设计所需的任何滤波器，若N取偶数时，则只可以设计低通、带通滤波器。

另外检查设计的滤波器是否满足技术指标也十分关键，有时设计出来的滤波器通带内允许最大衰减Rp、阻带内允许最小衰减As并不满足所给的设计技术指标，
这时则需要重新选择或调整窗函数的类型，估算滤波器的阶数N。