数字信号处理实验报告6.docx
《数字信号处理实验报告6.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验报告6.docx（15页珍藏版）》请在冰豆网上搜索。
数字信号处理实验报告6
实验报告
学院（系）名称：
计算机与通信工程学院
姓名
学号
专业
班级
实验项目
实验五FIR数字滤波器的设计
课程名称
数字信号处理
课程代码
实验时间
2013年05月31日
实验地点
主校区计算机基础实验室
批改意见
成绩
教师签字：
一，实验目的
加深理解FIR数字滤波器的时域特性和频率特性，掌握FIR数字滤波器的设计原理与设计方法，以及FIR数字滤波器的应用。
二，实验原理
FIR数字滤波器总是稳定的系统，且可以设计成具有线性相位的。
其在数据通信、图像处理、语言信号处理等实际应用领域中得到广泛的应用。
N阶有限冲激响应（FIR）数字滤波器的转移函数为：
，系统的单位脉冲响应h（n）是长度为N的有限长因果序列。
当满足h（n）=h（N-n-1）的对称条件时，该FIR数字滤波器具有线性相位。
FIR数字滤波器的设计方法主要有窗函数法和频率采样法。
1．窗函数法
FIR滤波器的冲激响应就是系统函数各次项的系数，所以设计FIR了不起的方法之一就是：
从时域出发，截取有限长的一段冲激响应作为H（z）系数，冲激响应长度N就是系统函数H（z）阶数。
只要N足够长，并且截取的方法合理，总能够满足频域的要求，这就是FIR滤波器的窗函数设计法。
2、频率采样法
频率采样法是先对理想频响
进行采样，得到样值
再利用插值公式直接求出系统转换函数
，以便实现；或者求出频响
，以便与理想频响进行比较。
在[0,2π]区间上对
进行N点采样，等效于时域以N为周期延拓。
三，实验内容与方法
【例3-5-1】数字滤波器的技术指标如下：
、
、
、
，采用窗函数法来设计一个FIR数字滤波器。
海明窗和布莱克曼窗均可以提供大于50dB的衰减。
如果选用海明窗设计，则它提供较小的过渡带，因此，其具有较小的阶数。
尽管在设计中用不到通带波动值
，但必须检查设计的实际波动，即验证它是否确实在给定容限内。
MATLAB程序如下：
wp=0.4*pi;ws=0.6*pi;deltaw=ws-wp;
N0=ceil（6.6*pi/deltaw）+1;
N=N0+mod（N0+1,2）;
wdham=（hamming（N））';
wc=（ws+wp）/2;
tao=（N-1）/2;
n=[0:
N-1];
m=n-tao+eps;
hd=sin（wc*m）./（pi*m）;
h=hd.*wdham
[H,w]=freqz（h,[1],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（h,[1],w）;
dw=2*pi/1000;
Rp=-（min（db（1:
wp/dw+1）））;
As=-round（max（db（ws/dw+1:
501）））;
n=0:
N-1;
subplot（2,2,1）;stem（n,hd,'.'）;title（'理想脉冲响应'）
axis（[0N-1-0.20.3]）;ylabel（'hd（n）'）;text（N+1,-0.1,'n'）
subplot（2,2,2）;stem（n,wdham,'.'）;title（'海明窗'）
axis（[0N-101.1]）;ylabel（'w（n）'）;text（N+1,0,'n'）
subplot（2,2,3）;stem（n,h,'.'）;title（'实际脉冲响应'）
axis（[0N-1-0.20.3]）;xlabel（'n'）;ylabel（'h（n）'）
subplot（2,2,4）;plot（w/pi,db）;title（'幅度响应（单位：
dB）'）;grid
axis（[01-10010]）;xlabel（'频率（单位：
pi）'）;ylabel（'分贝'）
set（gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]）;
set（gca,'YTickmode','manual','YTick',[-50,0]）
set（gca,'YTickMode','manual','YTickLabels',['50';'0']）
set（gcf,'color','w'）;
计算结果为：
N=35，
=0.04dB,
=52dB,满足要求，并且对通带波动的验证显示为
=0.04dB，它也是满足要求的。
可以看出用海明窗设计的FIR数字滤波器是满足要求的，且时域和频域的曲线如下图所示。
此外，在信号处理工具箱中，MATLAB还直接提供了一个子程序即firl，它利用窗函数法设计FIR滤波器，其标准调用格式为：
B=firl（M,WN,‘type’,window）
其中：
b为待设计的滤波器
系数向量，其长度为N=M+1；M为所选的滤波器阶数；
为滤波器给定的边缘频率，可以是标量，也可以是一个数组；‘type’为滤波器的类型，如高通、带通、带阻等，缺省时为低通；window为选定的窗函数类型，缺省时为Hamming窗。
同上例，MATLAB程序如下：
b=firl（34,1/pi,hamming（35））;
[H,w]=freqz（b,1,512）;
H_db=20*log10（abs（H））;
subplot（2,1,1）;stem（b）;
title（'hamming'）;
subplot（2,1,2）;plot（w,H_db）;
title（'frequency'）;
gridon;
得到的结果如下图所示。
2、频率采样法
【例3-5-2】同例3-5-1数字滤波器的技术指标如下：
、
、
、
，利用频率设计法设计一个FIR数字滤波器。
MATLAB程序如下：
N=35;wp=0.4*pi;ws=0.6*pi;wc=（wp+ws）/2;
N=N+mod（N+1,2）;
N1=fix（wc/（2*pi/N））;N2=N-2*N1-1;
A=[ones（1,N1+1）,zeros（1,N2）,ones（1,N1）];
theta=-pi*[0:
N-1]*（N-1）/N;
H=A.*exp（j*theta）;
h=real（ifft（H））;
wp1=2*pi/N*fix（wc/（2*pi/N）-1）;ws1=wp1+8*pi/N;
[H,w]=freqz（h,[1],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（h,[1],w）;
N=length（h）;L0=（N-1）/2;L=floor（L0）;
n=1:
L+1;ww=[0:
511]*pi/512;
ifall（abs（h（n）-h（N-n+1））<1e-8）
Ar=2*h（n）*cos（（（N+1）/2-n）'*ww）-mod（N,2）*h（L+1）;
type=2-mod（N,2）;
elseifall（abs（h（n）+h（N-n+1））<1e-8）&（h（L+1）*mod（N-2）<=1e-8）
A=2*h（n）*sin（（（N+1）/2-n）'*ww）;
type=4-mod（N,2）;
elseerror（'错误：
这不是线性相位滤波器！
'）
end
dw=2*pi/1000;
Rp=-min（db（1:
fix（wp1/dw）+1））;
As=-round（max（db（fix（ws1/dw）+1:
501）））;
l=0:
N-1;wl=（2*pi/N）*l;
wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];
subplot（2,2,1）;plot（wl（1:
N）/pi,A（1:
N）,'.',wdl,Adl）;
axis（[0,1,-0.1,1.1]）;title（'频率样本'）
xlabel（''）;ylabel（'A（k）'）;
set（gca,'XTickMode','manual','XTick',chop（[0,0.5,1],2））;
set（gca,'YTickMode','manual','YTick',[0,1]）;grid
subplot（2,2,2）;stem（l,h,'.'）;axis（[-1,N,-0.1,0.5]）;grid
title（'脉冲响应'）;ylabel（'h（n）'）;text（N+1,-0.1,'n'）
subplot（2,2,3）;plot（ww/pi,Ar,wl（1:
N）/pi,A（1:
N）,'.'）;
axis（[0,1,-0.2,1.2]）;title（'幅频响应'）
xlabel（'频率（单位：
pi）'）;ylabel（'Ar（w）'）;grid
subplot（2,2,4）;plot（w/pi,db）;axis（[0,1,-50,10]）;
title（'幅度响应'）;xlabel（'频率（单位：
pi）'）;
ylabel（'分贝数'）;grid
set（gca,'XTickMode','manual','XTick',chop（[0,0.5,1],2））;
set（gca,'YTickMode','manual','YTick',[-As,0]）;
set（gca,'YTickLabelMode','manual','YTickLabels',['As';'0']）;
set（gcf,'color','w'）;
经过计算可以得到：
=1.6011dB、
=24dB,所设计的滤波器的时域与频域特性如下图所示。
可以看出，其不满足设计指标。
为此，可以在构造的频率响应间断点处增加一个或若干了过度点来改善带外衰减指标。
在此选择增加两个过度点
、
来改善滤波器的设计。
MATLAB程序如下：
N=35;wp=0.4*pi;ws=0.6*pi;wc=（wp+ws）/2;
N=N+mod（N+1,2）;
N1=fix（wc/（2*pi/N））;N2=N-2*N1-1;
T1=input（'过渡带第一样本T1的值'）;
T2=input（'过渡带第二样本T2的值'）;
A=[ones（1,N1）,T1,T2,zeros（1,N2-2）,T2,T1,ones（1,N1-1）];
A=[ones（1,N1+1）,zeros（1,N2）,ones（1,N1）];
theta=-pi*[0:
N-1]*（N-1）/N;
H=A.*exp（j*theta）;
h=real（ifft（H））;wp1=2*pi/N*fix（wc/（2*pi/N）-1）;ws1=wp1+8*pi/N;
[H,w]=freqz（h,[1],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（h,[1],w）;
N=length（h）;L0=（N-1）/2;L=floor（L0）;
n=1:
L+1;ww=[0:
511]*pi/512;
ifall（abs（h（n）-h（N-n+1））<1e-8）
Ar=2*h（n）*cos（（（N+1）/2-n）'*ww）-mod（N,2）*h（L+1）;
type=2-mod（N,2）;
elseifall（abs（h（n）+h（N-n+1））<1e-8）&（h（L+1）*mod（N-2）<=1e-8）
A=2*h（n）*sin（（（N+1）/2-n）'*ww）;
type=4-mod（N,2）;
elseerror（'错误：
这不是线性相位滤波器！
'）
end
dw=2*pi/1000;
Rp=-min（db（1:
fix（wp1/dw）+1））;
As=-round（max（db（fix（ws1/dw）+1:
501）））;
l=0:
N-1;wl=（2*pi/N）*l;
wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];
subplot（2,2,1）;plot（wl（1:
N）/pi,A（1:
N）,'.',wdl,Adl）;
axis（[0,1,-0.1,1.1]）;title（'频率样本'）
xlabel（''）;ylabel（'A（k）'）;
set（gca,'XTickMode','manual','XTick',chop（[0,0.5,1],2））;
set（gca,'YTickMode','manual','YTick',[0,1]）;grid
subplot（2,2,2）;stem（l,h,'.'）;axis（[-1,N,-0.1,0.5]）;grid
title（'脉冲响应'）;ylabel（'h（n）'）;text（N+1,-0.1,'n'）
subplot（2,2,3）;plot（ww/pi,Ar,wl（1:
N）/pi,A（1:
N）,'.'）;
axis（[0,1,-0.2,1.2]）;title（'幅频响应'）
xlabel（'频率（单位：
pi）'）;ylabel（'Ar（w）'）;grid
subplot（2,2,4）;plot（w/pi,db）;axis（[0,1,-50,10]）;
title（'幅度响应'）;xlabel（'频率（单位：
pi）'）;
ylabel（'分贝数'）;grid
set（gca,'XTickMode','manual','XTick',chop（[0,0.5,1],2））;
set（gca,'YTickMode','manual','YTick',[-As,0]）;
set（gca,'YTickLabelMode','manual','YTickLabels',['As';'0']）;
set（gcf,'color','w'）;
运行程序画出幅频特性曲线如图所示。
同窗函数一样，MATLAB信号处理工具箱中提供了频率采样法的设计函数fir2。
它的典型调用法为：
b=fir2（M,f,m）
其次，b为待设计的滤波器系数向量，其长度为N=M+1；M为所选的滤波器阶数；f制定归一化的各频率边界频率，从0到1递增，1对应采样频率的一半，及奎奈斯特频率；m为制定各边界处的幅度响应，因此f和m的长度相等。
f=[01/pi1/pi1];m=[1100];
b=fir2（34,f,m）;
[H,w]=freqz（b,1,128）;
H_db=20*log10（abs（H））;
subplot（2,1,1）;stem（b,'.'）;
title（'实际脉冲响应'）;
subplot（2,1,2）;plot（w/pi,H_db）;
title（'幅度响应'）;
gridon
set（gcf,'color','w'）;
得到结果如图所示。
四、程序设计实验
数字滤波器的技术指标如下：
、
、
、
，分别采用窗函数法和频率采样法设计一个FIR数字滤波器。
采用窗函数法MATLAB程序如下：
wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp;
N0=ceil（6.6*pi/deltaw）+1;
N=N0+mod（N0+1,2）;
wdham=（hamming（N））';
wc=（ws+wp）/2;
tao=（N-1）/2;
n=[0:
N-1];
m=n-tao+eps;
hd=sin（wc*m）./（pi*m）;
h=hd.*wdham
[H,w]=freqz（h,[1],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（h,[1],w）;
dw=2*pi/1000;
Rp=-（min（db（1:
wp/dw+1）））;
As=-round（max（db（ws/dw+1:
501）））;
n=0:
N-1;
subplot（2,2,1）;stem（n,hd,'.'）;title（'理想脉冲响应'）
axis（[0N-1-0.20.3]）;ylabel（'hd（n）'）;text（N+1,-0.1,'n'）
subplot（2,2,2）;stem（n,wdham,'.'）;title（'海明窗'）
axis（[0N-101.1]）;ylabel（'w（n）'）;text（N+1,0,'n'）
subplot（2,2,3）;stem（n,h,'.'）;title（'实际脉冲响应'）
axis（[0N-1-0.20.3]）;xlabel（'n'）;ylabel（'h（n）'）
subplot（2,2,4）;plot（w/pi,db）;title（'幅度响应（单位：
dB）'）;grid
axis（[01-10010]）;xlabel（'频率（单位：
pi）'）;ylabel（'分贝'）
set（gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]）;
set（gca,'YTickmode','manual','YTick',[-50,0]）
set（gca,'YTickMode','manual','YTickLabels',['50';'0']）
set（gcf,'color','w'）;
程序运行结果如图
采用频率采样法MATLAB程序如下：
N=35;wp=0.2*pi;ws=0.3*pi;wc=（wp+ws）/2;
N=N+mod（N+1,2）;
N1=fix（wc/（2*pi/N））;N2=N-2*N1-1;
A=[ones（1,N1+1）,zeros（1,N2）,ones（1,N1）];
theta=-pi*[0:
N-1]*（N-1）/N;
H=A.*exp（j*theta）;
h=real（ifft（H））;
wp1=2*pi/N*fix（wc/（2*pi/N）-1）;ws1=wp1+8*pi/N;
[H,w]=freqz（h,[1],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（h,[1],w）;
N=length（h）;L0=（N-1）/2;L=floor（L0）;
n=1:
L+1;ww=[0:
511]*pi/512;
ifall（abs（h（n）-h（N-n+1））<1e-8）
Ar=2*h（n）*cos（（（N+1）/2-n）'*ww）-mod（N,2）*h（L+1）;
type=2-mod（N,2）;
elseifall（abs（h（n）+h（N-n+1））<1e-8）&（h（L+1）*mod（N-2）<=1e-8）
A=2*h（n）*sin（（（N+1）/2-n）'*ww）;
type=4-mod（N,2）;
elseerror（''）
end
dw=2*pi/1000;
Rp=-min（db（1:
fix（wp1/dw）+1））;
As=-round（max（db（fix（ws1/dw）+1:
501）））;
l=0:
N-1;wl=（2*pi/N）*l;
wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];
subplot（2,2,1）;plot（wl（1:
N）/pi,A（1:
N）,'.',wdl,Adl）;
axis（[0,1,-0.1,1.1]）;title（'频率样本'）
xlabel（''）;ylabel（'A（k）'）;
set（gca,'XTickMode','manual','XTick',chop（[0,0.5,1],2））;
set（gca,'YTickMode','manual','YTick',[0,1]）;grid
subplot（2,2,2）;stem（l,h,'.'）;axis（[-1,N,-0.1,0.5]）;grid
title（'脉冲响应'）;ylabel（'h（n）'）;text（N+1,-0.1,'n'）
subplot（2,2,3）;plot（ww/pi,Ar,wl（1:
N）/pi,A（1:
N）,'.'）;
axis（[0,1,-0.2,1.2]）;title（'幅频响应'）
xlabel（'频率（单位：
pi）'）;ylabel（'Ar（w）'）;grid
subplot（2,2,4）;plot（w/pi,db）;axis（[0,1,-50,10]）;
title（'幅度响应'）;xlabel（'频率（单位：
pi）'）;
ylabel（'分贝数'）;grid
set（gca,'XTickMode','manual','XTick',chop（[0,0.5,1],2））;
set（gca,'YTickMode','manual','YTick',[-As,0]）;
set（gca,'YTickLabelMode','manual','YTickLabels',['As';'0']）;
set（gcf,'color','w'）;
程序运行结果如图
五，实验预习要求
（1）预习实验原理
（2）熟悉实验程序
（3）思考程序设计实验部分程序的编程
六，实验报告要求
（1）在MATLAB中输入程序，验证实验结果，并将实验结果存入指定存储区域中。
（2）对于程序设计实验，要求通过对验证性实验的练习，自行编制完整的实验程序，实现对信号的模拟，并得到实验结果。
（3）在实验报告中写出完整的自编程序，并给出实验结果。
七，思考题
（1）窗函数法和频率采样法的优缺点分别是什么？
答：
窗函数法优点：
1、频域方均误差最小。
2、无稳定性问题。
3、可以设计各种特殊类型的滤波器。
4、方法特别简单。
缺点：
1、不易控制边缘频率。
2、幅频性能不理想。
3、h（n）较长。
频率采样法优点：
1、函数插值法逼近。
2、可以直接从频域出发设计，比较简单。
缺点：
周期不好设定，需要试凑来找到局部最优组合。
（2）在FIR窗函数设计中，为何采用不同特性的窗函数？
选用窗函数的依据是什么？
答：
不同特性窗函数对滤波器频率特性的影响不同。
选用窗函数的依据是主瓣宽度、旁瓣最大峰值电平、旁瓣谱峰衰减速度。
（3）在频率采样法中，如果阻带衰减不够，应采取什么样的措施？
答：
在边界频率处增加一个过渡点。
为了保证过渡带宽不变，将采样点数增加一倍，并将过渡点的采样值进行优化。
升级会员 