窗函数数字信号课设.docx
《窗函数数字信号课设.docx》由会员分享，可在线阅读，更多相关《窗函数数字信号课设.docx（19页珍藏版）》请在冰豆网上搜索。
窗函数数字信号课设
河北科技大学
课程设计报告
学生姓名：
学号：
专业班级：
电子信息工程
课程名称：
数字信号处理课程设计
学年学期
2013——2014学年第一学期
指导教师：
张秀清
2014年6月
课程设计成绩评定表
学生姓名
学号
成绩
专业班级
起止时间
2014-6-9至2014-6-12
设计题目
1.窗函数设计低通滤波器
2.用哈明窗设计FIR带通数字滤波器
指导教师评语
指导教师：
年月日
1.窗函数设计低通滤波器
1.1设计目的1
1.2设计原理推导与计算1
1.3设计内容与要求2
1.4设计源程序与运行结果3
1.5思考题10
2.用哈明窗设计FIR带通数字滤波器
2.1设计要求
14
2.2设计原理和分析
14
2.3详细设计
15
2.4调试分析及运行结果15
2.5心得体会17
参考文献17
1.窗函数设计低通滤波器
1.1设计目的
1.熟悉设计线性相位数字滤波器的一般步骤。
2.掌握用窗函数法设计FIR数字滤波器的原理和方法。
3.熟悉各种窗函数的作用以及各种窗函数对滤波器特性的影响。
4.学会根据指标要求选择合适的窗函数。
1.2设计原理推导与计算
如果所希望的滤波器的理想的频率响应函数为Hde''，则其对应的单位脉
冲响应为
1
hdnd—He，"（4」）
窗函数设计法的基本原理是设计设计低通FIR数字滤波器时，一般以理想低
N-1
2
通滤波特性为逼近函数He「-，即
Hd（ej*
二Hdej'ejd1ZQ：
d‘=空Ja
2—二n-a
用有限长单位脉冲响应序列
hn逼近hdn。
由于hdn往往是无限长序列,
而且是非因果的，所以用窗函数
「n将hdn截断，并进行加权处理，得到:
（4.2）
hn就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函
数Hej「为
（4.3）
N-4
Hej'八hnejn
n=0
式中，N为所选窗函数」n的长度。
用窗函数法设计的滤波器性能取决于窗函数「n的类型及窗口长度N的取
值。
设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。
各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见
）0
窗函数
旁瓣峰值
幅度/dB
过
渡带宽
阻带最小
衰减/dB
过渡带带
宽沧N）
矩形窗
-13
4n
/N
-12
0.9
三角形窗
-25
8n
/N
-25
2.1
汉宁窗
-31
8n
/N
-44
3.1
哈明窗
-41
8n
/N
-53
3.3
不莱克曼窗
-57
12n/N
-74
5.5
凯塞窗
（B=7.865）
-57
10n/N
-80
5
表
（一）各种窗函数的基本参数
这样选定窗函数类型和长度N之后，求出单位脉冲响应hn=hdn・•n,并按照式（4.3）求出He八oHe八是否满足要求，如果He八不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。
1.3设计内容与要求
（一）设计要求：
1.学会计算滤波器各项性能指标及如何来满足给定的指标要求。
2.用MATLAB©言编程实现给定指标要求的滤波器的设计。
3.熟悉MATLA语言，独立编写程序。
4.设计低通FIR滤波器的指标：
通带最大波动
Rp=0.25dB,©p=0.2兀
阻带最小衰减
（二）、设计内容：
A=50dB,0s=0・3兀
1.熟悉各种窗函数,在MATLAB^令窗下浏览各种窗函数，绘出（或打印）所看到的窗函数图。
2.编写计算理想低通滤波器单位抽样响应hd（n）的m函数文件ideal.m。
3.编写计算N阶差分方程所描述系统频响函数H（er）的m函数文件fr.m。
4.根据指标要求选择窗函数的形状与长度No（至少选择两种符合要求的窗
函数及其对应的长度）。
5.编写.m程序文件，通过调用ideal.m和fr.m文件，计算你设计的实际
低通FIR滤波器的单位抽样响应h（n）和频率响应H（ejJ，打印在频率区间［0,
n］上的幅频响应特性曲线
6.验证所设计的滤波器是
H（e啲）~灼，幅度用分贝表示
否满足指标要求。
1.4设计的源程序及运行结果：
1、利用MATLA窗口观察各种窗函数:
%巴特利特窗
w=bartlett（20）;
subplot（3,2,1）;
plot（w）;
stem（w,'y'）;%'y'表示黄色
%stem表示以离散图输出
title（'巴特利特床窗'）;
xlabel（'n'）;%横坐标为n
ylabel（'w（n）'）;%纵坐标为
w（n）
%布莱克曼窗
subplot（3,2,2）;
plot（w）;
stem（w,'b'）;%'b'表示蓝色title（'布莱克曼窗'）;xlabel（'n'）;
ylabel（'w（n）'）;
%矩形窗
w=boxcar（20）;
subplot（3,2,3）;
plot（w）;
stem（w,'r'）;
title（'矩形窗'）;
w=blackman（20）;
xlabel（'n'）;
ylabel（'w（n）'）;
%S明窗
w=hamming（20）;
plot（w）;
stem（w,'m'）;%'m'表示紫色title（'海明窗'）;
xlabel（'n'）;
ylabel（'w（n）'）;
冰宁窗
w=hanning（20）;
subplot（3,2,5）;
plot（w）;
stem（w,'g'）;%'g'表示绿色
title（'汉宁窗'）;
xlabel（'n'）;
ylabel（'w（n）'）;
%凯泽窗
beta=5.6533;
w=kaiser（20,beta）;
subplot（3,2,6）;
plot（w）;
stem（w,'k'）;%'k'表示黑色
title（'凯泽
窗,beta=5.6533'）;
xlabel（'n'）;
ylabel（'w（n）'）;
巴特利特床窗
布莱克曼窗
矩形窗
n
海明窗
n
汉宁窗
凯泽窗，beta=5.6533
常用窗函数的图形
2、理想低通滤波器单位抽样响应hd（n）的m函数文件ideal.m
functionhd=ideal（wc,M）
%理想低通滤波器计算
%hd为0到M-1之间的理想脉冲响应
%w（为截止频率
%丽理想滤波器的长度
alpha=（M-1）/2;
n=0:
M-1;
m=n-alpha+eps;
hd=sin（wc*m）./（pi*m）;
3、N阶差分方程所描述的系统频响函数的m函数文件fr.m
function[db,mag,pha,gfd,w]=fr（b,a）
%求解系统响应
%db为相位振幅（db）
%ma为绝对振幅
%pha为相位响应
%grd为群延时
%w为频率样本点矢量
%匕为Ha（z）分析多项式系数（对FIR而言，b=h）
%&为Hz（z）分母多项式系数（对FIR而言，a=1）
[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）;
gfd=grpdelay（b,a,w）;
4、实际低通滤波器FIR:
%5海明窗设计低通滤波器
wp=0.2*pi;
ws=0.3*pi;
tr_width=ws-wp;
disp（['海明窗设计低通滤波器参数：
']）;
M=ceil（6.6*pi/tr_width）+1;
disp（['滤波器的长度为’,num2str（M）]）;
n=0:
M-1;
wc=（ws+wp）/2;%理想LPF的截止频率
hd=ideal（wc,M）;
w_ham=（hamming（M））';
h=hd.*w_ham;
[db,mag,pha,gfd,w]=fr（h,[1]）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））;%求出实际通带波动
disp（['实际带通波动为’,num2str（Rp）]）;
As=-round（max（db（ws/delta_w+1:
1:
501）））;%求出最小阻带衰减
disp（['最小阻带衰减为-',num2str（As）],'db'）;
%绘图
subplot（1,1,1）
subplot（2,6,1）
stem（n,hd）;
title（'理想冲击响应'）;
axis（[0M-1-0.10.3]）;
ylabel（'hd（n）'）;
subplot（2,6,2）
stem（n,w_ham）;
title（'海明窗'）;
axis（[0M-101.1]）;
ylabel（'w（n）'）;
subplot（2,6,7）
stem（n,h）;
title（'实际冲激响应'）;
axis（[OM-1-0.10.3]）;
xlabel（'n'）;
ylabel（'h（n）'）;
subplot（2,6,8）
plot（w/pi,db）;
title（'幅度响应（db）'）;
axis（[01-10010]）;
grid;
xlabel（'以pi为单位的频率'）;
ylabel（'分贝数'）;
图
（1）海明窗设计的FIR
海明窗设计低通滤波器参数
滤波器的长度为67
实际带通波动为0.03936
最小阻带衰减为-52db
%用布莱克曼窗设计低通滤波器
wp=0.2*pi;
ws=0.3*pi;
tr_width=ws-wp;
disp（['布莱克曼窗设计低通滤波器的参数：
']）;
M=ceil（11.0*pi/tr_width）+1;
disp（['滤波器的长度为’,num2str（M）]）;
n=0:
M-1;
%!
想LPF的截止频率
wc=（ws+wp）/2;
hd=ideal（wc,M）;
w_bla=（blackman（M））';
h=hd.*w_bla;
[db,mag,pha,gfd,w]=fr（h,[1]）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））;%求出实际通带波动
disp（['实际带通波动为’,num2str（Rp）]）;
As=-round（max（db（ws/delta_w+1:
1:
501）））;%求出最小阻带衰减
disp（['最小阻带衰减-',num2str（As）],'db'）;
%绘图
subplot（2,6,3）
stem（n,hd）;
title（'理想冲击响应'）;
axis（[0M-1-0.10.3]）;
ylabel（'hd（n）'）;
subplot（2,6,4）
stem（n,w_bla）;
title（'布莱克曼窗'）;
axis（[0M-101.1]）;
ylabel（'w（n）'）;
subplot（2,6,9）
stem（n,h）;
title（'实际冲激响应'）;
axis（[OM-1-0.10.3]）;
xlabel（'n'）;
ylabel（'h（n）'）;
subplot（2,6,10）
plot（w/pi,db）;
title（'幅度响应（db）'）;
axis（[01-10010]）;
grid;
xlabel（'以pi为单位的频率'）;ylabel（'分贝数'）;
图
（2）布莱克曼窗设计的FIR
布莱克曼窗设计低通滤波器的参数：
滤波器的长度为111
实际带通波动为0.0033304
最小阻带衰减为-73db
5、技术指标比较：
（1）海明窗设计低通滤波器参数：
滤波器的长度为67
实际带通波动为0.03936
最小阻带衰减为-52db
（2）布莱克曼窗设计低通滤波器的参数：
滤波器的长度为111
实际带通波动为0.0033304
最小阻带衰减为-73db
在相同的技术指标下用布莱克曼窗设计的低通滤波器实际带通波动实际带通波动最小，最小阻带衰减，滤波器的长度最大；海明窗和凯泽窗最小阻带衰减差不多，滤波器的长度页差不多，但是海明窗实际波动小于凯泽窗；所以用布莱
克曼窗用设计的FIR最逼近理想单位冲击响应。
这三个窗设计的低通滤波器都符合要求。
1.5思考题：
1.设计线性相位数字滤波器的一般步骤：
（1）、给定所要求的频率响应函数Hd（ejw）以及技术指标阻带衰减§和皿;
（2）、求hd=IDTFT[Hd（ejw）];
由过渡带带宽及带阻最小衰减的要求，利用六种常见的窗函数基本参数的比较表或凯泽窗设计FIR滤波器的经验公式，选择窗函数的形状及N的大小（一般N要通过几次试探而最后确定）；
（3）求得所设计的FIR滤波器的单位抽样响应：
h（n）=hd（n）w（n）,N=0,1,2,3…,N-1
（4）、求Hd（ejw）=DIFT[h（n）],校验是否满足设计要求，如果不满足，则重新设计
2.窗函数有哪些指标要求？
对给定指标要求的低通滤波器，理论计算所需窗函数的长度N。
答：
窗函数的指标要求：
主瓣宽度，旁瓣峰值。
海明窗设计低通滤波器的长度为67,凯泽窗设计低通滤波器的长度为60，布莱克曼窗设计低通滤波器的长度为111。
3.用窗函数法设计FIR滤波器，滤波器的过渡带宽度和阻带衰减与哪些因素有关？
答：
过渡带宽度与窗函数的形状和窗的宽度有关；阻带衰减只有窗函数的形状决定，不受N的影响。
4.计算理想带通滤波器单位抽样响应hd（n）时取N为奇数和N为偶数有没有区别？
你计算时所用的方法是仅适合于N为奇数或偶数的一种还是两种都可以用？
答:
以海明窗为例说明：
滤波器的长度为67
实际带通波动为0.03936
最小阻带衰减为52
滤波器的长度为68
实际带通波动为0.036424
最小阻带衰减为53
滤波器的长度为64
实际带通波动为0.068677
最小阻带衰减为45
N取奇偶数时，低通滤波器的幅度函数是不同的，如上图所示（海明窗设计低通滤波器N取奇偶数时的图），通过比较，当滤波器的长度大于技术指标要求的长度时，选择偶数也是满足要求的，所以海明窗奇偶都满足。
同理，N取奇偶数，选择其他的窗函数也满足设计低通滤波器的指标要求。
比较所选窗长N相同但窗形状不同对滤波器设计结果的影响以及选同一种
窗函数但窗长N不同时对滤波器设计结果的影响，将结论写在报告中。
海明窗
布莱克曼窗
凯泽窗
以海明窗设计的低通
滤波器的长度为67
滤波器的长度为67
滤波器的长度为67
滤波器的长度N=67,
实际带通波动0.03936
实际带通波动0.29523
实带通波动0.047919
用布莱克曼窗和凯泽
最小阻带衰减为52
最小阻带衰减为30
最小阻带衰减为51
窗设计低通滤波器
线性斜率beta=4.5513
以凯泽窗设计的低通
滤波器的长度为60
滤波器的长度为60
滤波器的长度为60
滤波器的长度N=60,
实际带通波动为0.12166
实际带通波动为0.47079
实际带通波动0.053711
用布莱克曼窗和海明
最小阻带衰减为39
最小阻带衰减为26
最小阻带衰减为51
窗窗设计低通滤波器
线性斜率beta=4.5513
以布莱克曼窗设计的
滤波器的长度为
111
滤波器的长度为111
滤波器的长度为111
低通滤波器的长度
实际带通波动为
0.02894
实际带通波动0.0033304
实际带通波动0.034426
N=111,用凯泽窗和海
最小阻带衰减为
59
最小阻带衰减为73
最小阻带衰减为58
明窗窗设计低通滤波
线性斜率beta=4.5513
器
结论:
1、当以海明窗设计的低通滤波器的长度N,用布莱克曼窗和凯泽窗设计低通滤波器，即N=67时，三个窗函数设计的低通滤波器中布莱克曼窗设计的不满
足要求，凯泽窗设计的带通波动大于海明窗设计的低通滤波器带通波动，最小阻
带衰减小于海明窗设计的低通滤波器。
凯泽窗设计的结果没有用海明窗设计结果好。
当N取67、60、111时，用海明窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N<67时不满足设计技术指标。
2、以凯泽窗设计的低通滤波器的长度N,用布莱克曼窗和海明窗窗设计低通滤波器，即N=60时，海明窗设计的低通滤波器的最小阻带衰减不满足指标要求，布莱克曼窗设计的低通滤波器的最小阻带衰减和带通波动都不满足要。
当N取67、60、111时，用布莱克曼窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N<111时不满足设计技术指标。
3、以布莱克曼窗设计的低通滤波器的长度N,用凯泽窗和海明窗窗设计低通滤波器，即N=111时，三个都满足要求。
但是用布莱克曼窗函数设计的低通滤波器的性能更好。
当N取67、60、111时，用凯泽窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N>60时都满足设计技术指标。
2、用哈明窗设计FIR带通数字滤波器
2.1设计要求
针对一个含有5Hz、15Hz和30HZ勺混合正弦波信号x（t），设计滤波系统，滤除5Hz和30HZ的正弦分量，阻带的最小衰减不小于50dB。
1确定采样频率fs、采样长度N,对x（t）进行采样；
2选择合适的窗函数，设计FIR带通数字滤波器；
3通过计算机仿真对滤波器的性能进行分析。
2.2设计原理分析和设计
此题要求也是滤波，但是要求滤掉5Hz和30HZ的信号，而保留15HZ的信号，因此必须采用带通滤波器。
要求阻带的最小衰减不低于50dB,查表可知可选哈明窗。
分析题目可知：
f1=5Hz，T1=1/5s;
f2=15Hz，T2=1/15s;
f3=30Hz，T3=1/30s;
则采样时长tT至少应为0.2s，取tT=1s;
采样频率fs>2fc;可取fs=150;则采样间隔T=1/fs;所以采样长度N=tT/T=150;
数字技术指标可取
wp1=2*pi*14/fs;
wp2=2*pi*23/fs;ws1=2*pi*7/fs;
ws2=2*pi*27/fs;
滤波器设计好后，利用卷积和可得输出信号y（n）=x（n）*h（n），卷积和长度
N=N1+N2-1;
2.3详细设计
源程序代码:
tT=1;%采样时长
fs=150;%采样频率
wp1=2*pi*14/fs;
wp2=2*pi*20/fs;
ws1=2*pi*7/fs;
ws2=2*pi*27/fs;
trwidth=wp1-ws1;%过渡带带宽
2=ceil（8*pi/trwidth）%计算滤波器阶次
3dB截止频率
理想带通滤波器的单位响应
n=0:
N1-1;
wc1=（ws1+wp1）/2;wc2=（ws2+wp2）/2;%计算
alpha=（N1-1）/2;%单位响应的对称中心
m=n-alpha;
hd=sin（wc2*m）./（pi*m）-sin（wc1*m）./（pi*m）;%wham=（hamming（N1））';
h=hd.*wham;%实际带通滤波器单位响应
w=0:
0.01:
pi;
H=freqz（h,1,w）;%单位滤波器幅频响应
T=1/fs;%对输入信号进行采样
N2=fs*tT;
t=（0:
N2-1）*T;
输出信号'）
2.4调试分析及运行结果
代码完成后，滤波器幅频响应、单位响应和输入信号能够正常输出，但是绘制输出信号时用的是plot（y），并未指明自变量，结果绘制的图（如下）是以n为自变量，看起来非常不明了也不清楚。
输出信号
后将绘图语句直接改为plot（t,y），结果出现错误
提示t和y矩阵长度不匹配，翻阅《信号与系统》上册课本重新复习了离散卷积和部分，得知两信号卷积后长度变化规律，于是将绘制输出信号部分代码改为：
y=conv（x,h）;
N=N1+N2-1;
tt=（0:
N-1）*T;
plot（tt,y）;
然后得到以时间为横坐标的输出信号。
运行结果如下：
"■IFigure1I®I回
FileEditViewInsertToolsDesktopWindowHelp*
□4'耳、＜?
倉归必「旦口目Q骂
数字滤浪器幅频响应
0.1020.30.40.50.G0.70.80.91
数学滤漩器单位响应
210
1020304050B0708090
軸入倍号
0.160.2
辐出信号
026
0.3
0.350.4
420-22
0
0.1
0.2
0.3
0.4
0.5
0.6
07
o.a
2.5心得体会：
通过运用MATLABS计低通滤波器，学习了MATLAB语言的编程以及熟悉MATLAB^件环境，会利用MATLAB®写基本程序，熟悉绘制图形命令的正确使用，以及掌握MATLA处理数字信号的有关命令。
在命令窗口运行编写好的m文件，
学习M文件调用其他的M文件，文件之间的调用必须满足这两个文件在一个文件夹里。
学习生成常用窗函数的MATLA语言函数，并通过调用生成各种窗函数MATLAG内部函数，以及编写窗函数M文件，在figure窗口观察常用窗函数的图形。
掌握FIR带通滤波器设计原理，并自己查阅有关带通滤波器的设计的书籍，根据设计指标要求利用MATLA软件设计出符合要求的带通滤波器。
另外，掌握设计FIR滤波器的一般步骤，以及掌握窗函数的选择。
在MATLA窗口观察用不
同窗函数设计的带通滤波器的图形，比较在相同的指标条件下不同的窗函数设计的带通滤波器自之间的差别，以及N取不同值时，同一个窗函数设计的FIR的带通滤波器的技术指标不同。
参考文献
1、王华、李有军、刘建存，《MATLAB电子仿真与应用教程》，国防科技大学出版社，2007
2、闻新、周露、张鸿《MATLAB科学图形构建基础与应用》，科技出版社，2002
3、程佩青，数字信号处理及应用，清华大学出版社，2009
4、王宏，MATLA販其在信号处理中的应用，清华大学出版社，2004
升级会员 