利用MATLAB仿真软件系统结合窗函数法设计一个数字带通FIR滤波器.docx

目录
目录 1
摘要 I
Abstract II
1 原理说明 1
1.1 数字滤波技术 1
1.2 FIR滤波器 2
1.3 窗函数 3
2 滤波器设计 4
2.1 滤波器设计要求 4
2.2 设计函数的选取 4
2.3 窗函数构造 5
2.4 设计步骤 8
2.5 设计方法 8
3 滤波器测试 15
3.1 滤波器滤波性能测试 15
3.2 滤波器时延测量 16
3.3 滤波器稳定性测量 17
4 心得体会 19
5 参考文献 20
附件一： 21
附件二： 22

摘要

Abstract

1 原理说明
随着信息时代的到来，数字信号处理已经成为一门极其重要的学科和技术，并且在通信、语音、图像、自动控制等众多领域得到了广泛的应用。在数字信号处理中，数字滤波器占有极其重要的地位，它具有精度高、可靠性好、灵活性大等特点。现代数字滤波器可以用软件或硬件两种方式来实现。软件方式实现的优点是可以通过滤波器参数的改变去调整滤波器的性能。
MATLAB是一种面向科学和工程计算的语言，它集数值分析、矩阵运算、信号处理和图形显示于一体，具有编程效率高、调试手段丰富、扩充能力强等特点。MATLAB的信号处理工具箱具有强大的函数功能，它不仅可以用来设计数字滤波器，还可以使设计达到最优化，是数字滤波器设计的强有力工具。

1.1 数字滤波技术
数字滤波，就是通过一定的计算或判断程序减少干扰在有用信号中的比重，故实质上是一种程序滤波。与此对应的就是模拟滤波，由于模拟滤波牵扯到的其他知识太多在此不详细介绍了，模拟滤波主要无源绿波（直接用电阻、电容、电感等不外接电源的元件组成的）与有源滤波（如运算放大器等需要外接电源组成的），其目的是将信号中的噪音和干扰滤去或者将希望得到的频率信号滤出为我所用。数字滤波的出现克服了模拟滤波的很多不足，具有以下优点：
A．是用程序实现的，不需要增加硬设备，所以可靠性高，稳定性好。
B．可以对频率很低的信号实现滤波，克服了模拟滤波的缺陷。
C．可以根据信号的不同，采用不同的滤波方法或参数，具有灵活、方便、功能强的特点。
几种常用的滤波方法：
1． 算术平均值法
2． 中值滤波法
3． 滑动平均值法
4． 限幅滤波法
5． 惯性滤波法
数字滤波技术通过数字滤波器实现，从实现方法上可以分为FIR数字滤波器和IIR数字滤波器，按功能可分为低通滤波器（LPF）、高通滤波器（HPF）、带通滤波器（BPF）和带阻滤波器（BSF）。本文主要对FIR滤波器加以介绍。
1.2 FIR滤波器
，FIR滤波器具有以下主要优点：
1.FIR滤波器具有准确的线性相位；
2.FIR滤波器永远稳定；
3.FIR滤波器设计方法一般是线性的；
4.FIR滤波器在硬件上具有更高的运行效率；
5.FIR滤波器启动传输时间只需要有限时间。
FIR滤波器的主要缺点有：
1.FIR滤波器为达到同样的性能要求需要比IIR滤波器高得多的阶数；
2.相应的FIR滤波器的时延比同等性能的IIR滤波器高很多。
FIR滤波器的硬件实现主要有数字集成芯片，DSP芯片FIR滤波器，可编程FIR滤波器，后两者的实际方法主要通过MATLAB软件进行设计，其设计方法多样，形式灵活，能够满足各种要求，并且不受数字集成芯片规格的限制。
FIR滤波器的设计方法主要有窗函数法、多带和过渡带、约束最小二乘法、任意相应法、升余弦法，其中最常用的是窗函数法。
1.3 窗函数
窗函数法是设计FIR滤波器的最主要方法之一，实际中遇到的离散时间信号总是有限长的，因此不可避免的要遇到数据截短的问题，在信号处理中，对离散序列的截短是通过序列与窗函数相乘来实现的。
在信号处理中，窗函数是一种除在给定区间之外取值均为0的实函数。譬如：在给定区间内为常数而在区间外为0的窗函数被形象地称为矩形窗。任何函数与窗函数之积仍为窗函数，所以相乘的结果就像透过窗口“看”其他函数一样。窗函数在光谱分析、滤波器设计以及音频数据压缩等方面有广泛的应用。
常用的窗函数有矩形窗、巴特利特（Bartlett）窗、三角窗、海明（Hamming）窗、汉宁（Hanning）窗、布莱克曼（Blackman）窗、切比雪夫Chebyshev）窗、凯泽（Kaiser）窗。

2 滤波器设计
2.1 滤波器设计要求
利用MATLAB仿真软件系统结合窗函数法设计一个数字带通FIR滤波器。
要求：分别使用矩形窗、三角形窗、汉明窗、布莱克曼窗、凯泽窗进行设计，并输出滤波器的频率特性。
技术指标：
1． 采样频率为20kHz；
2． 通带边缘频率：fp1=4.5kHz,fp2=6.5kHz;
3． 通带峰值起伏：αp<1dB;
4． 阻带边缘频率：fs1=3k,fs2=7.5k;
5． 最小阻带衰减：As>40dB
2.2 设计函数的选取
MATLAB信号处理工具箱提供了基于加窗的线性相位FIR滤波器设计函数fir1和fir2，fir1函数的调用格式为：
b=fir1(n,Wn)
b=fir1(n,Wn,'ftype')
b=fir1(n,Wn,window)
b=fir1(n,Wn,'ftype',window)
b=fir1(…..,'normalization')
函数参数说明如下：
1．n表示滤波器的阶数
2．'ftype'表示所设计滤波器的类型：
3．'high'表示高通滤波器
4．'stop'表示带阻滤波器
5．'DC-1'表示多通带滤波器，第一频带为通带
6．'DC-0'表示多通带滤波器，第一频带为阻带；默认时为低通或带通滤波器；
7．'window'为窗函数，是长度为n+1的列向量，默认时函数自动取Hamming窗。
该函数实现加窗的线性相位FIR滤波器设计，可以设计标准低通、带通、高通和带阻滤波器（具有任意频率响应的加窗滤波器可以采用fir2进行设计）
2.3 窗函数构造
MATLAB工具箱已经提供了各种窗函数的构造函数，因而窗函数的构造十分方便，下面给出几种常用窗函数的构造方法：
1．矩形窗：利用w=boxcar(n)的形式得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，它的元素由窗函数的值组成。‘w=boxcar(n)’等价于‘w=ones(1,n)’.
2．三角窗：利用w=triang(n)的形式得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，它的元素由窗函数的值组成。w=triang(N-2)等价于bartlett(N)。
3．汉宁窗：利用w=hanning(n)得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，包含了窗函数的n个系数。
4．海明窗：利用w=hamming(n)得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，包含了窗函数的n个系数。它和汉宁窗的主瓣宽度相同，但是它的旁瓣进一步被压低。
5．布拉克曼窗：利用w=blackman(n)得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，包含了窗函数的n个系数。它的主瓣宽度是矩形窗主瓣宽度的3倍，为12*pi/N，但是它的最大旁瓣值比主瓣值低57dB。
6．切比雪夫窗：它是等波纹的，利用函数w=chebwin(N,R)方式设计出N阶的切比雪夫2窗函数，函数的主瓣值比旁瓣值高RdB,且旁瓣是等波纹的。
7．巴特里特窗：利用w=bartlett(n)的形式得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，包含了窗函数的n个系数。
8．凯泽窗：利用w=kaiser(n,beta)的形式得到窗函数。
2.3.1 窗函数设计条件
在使用窗函数设计滤波器时要满足以下两个条件：
1． 窗谱主瓣尽可能地窄，从而可以获得较陡峭的过渡带；
2．尽量减少窗谱的最大旁瓣的相对幅度，及尽可能是能量集中于主瓣，减少峰肩和波纹，进而增加阻带的衰减。
根据工程经验，给定的滤波器指标一般为通带截止频率ωp，阻带截止频率ωs，实际通带Rp，和最小阻带衰减As，窗函数的设计经验公式为：
归一化过渡带： 公式（2-1）
滤波器阶数： 公式（2-2）

2.3.2 窗函数设计条件
实际工程常用的窗函数有五种，即矩形窗、三角窗、汉宁窗、海明窗和凯泽窗。这些窗函数之间的性能比较如表2-1所示。
表2-1 5种窗函数性能比较
窗类型

旁瓣峰值
主瓣峰值
最小阻带衰减
矩形窗
13dB
4π/M
21dB
三角窗
25dB
8π/M
25dB
汉宁窗
31dB
8π/M
44dB
海明窗
41dB
8π/M
53dB
凯泽窗
57dB
12π/M
74dB
常用窗函数绘图比较：
在MATLAB中运行以下代码：

代码2-1：
n=50;
x=1:50;
juxing=boxcar(n); %构造矩形窗
sanjiao=triang(n); %构造三角窗
hanming=hamming(n); %构造汉明窗
bulaikeman=blackman(n); %构造布莱克曼窗
kaize=kaiser(n); %构造凯泽窗
plot(x,juxing,'b.',x,sanjiao,'gx',x,hanming,'r+',x,bulaikeman,'cd',x,kaize,'k*');
legend('矩形窗','三角窗','汉明窗','布莱克曼窗','凯泽窗');
运行结果如图2-1所示：

图2-1 5种窗函数绘图比较
2.4 设计步骤
实际利用窗函数法进行FIR滤波器设计时，依据所给的技术指标一般需要经过以下几个步骤进行设计：
1． 给定理想的频率响应函数Hd(ejw)及技术指标；
2． 求出理想的单位抽样响应hd(n);
3． 根据阻带衰减选择窗函数w(n)
4． 根据过渡带宽度确定N值;
5． 求出所设计的FIR滤波器的单位抽样响应；
6． 计算频率响应，验算指标是否满足要求。
2.5 设计方法
MATLAB作为一款优秀的数值计算软件，本身就内置了丰富的函数，其中便有用于通信仿真的一系列函数，并且MATLAB中还集成了通信设计的工具箱，不管是内置的函数，还是通信工具箱，均有专用于滤波器设计的工具，常用的主要有用函数法设计和用通信工具箱设计，下面分别予以介绍。
2.5.1 依据设计步骤编写M文件设计
此种方法不依赖MATLAB中的滤波器设计函数，而是依据FIR滤波器的设计步骤自己求解理想滤波器的冲击响应，然后用窗函数对冲击响应进行截短，从而得到FIR滤波器，由于低通滤波器设计较为简单，因而可以通过两个低通的理想冲击响应函数相减得到理想带通的冲击响应，再通过窗函数对其截短，从而得到实际可行的FIR滤波器，这里以用汉明窗进行截短实现带通FIR滤波器，关键代码如下（完整代码见附页）：
代码2-2：
wp1=0.3*pi; ws1=0.45*pi;
wp2=0.65*pi; ws2=0.75*pi;
tr_width=ws1-wp1; %求过渡带宽度
M=ceil(6.6*pi/tr_width)+1; %求得所需窗函数的长度
n=[0:1:M-1];
wc1=(ws1+wp1)/2; wc2=(ws2+wp2)/2; %求截止频率
hd=ideal_lp(wc2,M)-ideal_lp(wc1,M); %求得理想带通的冲击响应
w_ham=(hamming(M))'; %得到长度为M的汉明窗
h=hd .* w_ham; %利用窗函数截短
程序运行结果：

图2-2 汉明窗带通滤波器
2.5.2 利用MATLAB自带函数设计
利用原理说明介绍的fir1函数进行设计，这种设计方法只需要给出滤波器的阶数，截止频率，窗函数等参数，MATLAB即可自行完成设计，并可通过freqz函数查看滤波器的幅频响应和相频响应，已验证滤波器是否满足设计要求，下面给出利用fir1函数设计的不同窗函数的数字带通滤波器。
1. 利用矩形窗进行设计
代码2-2：
fs=20000; %设定采样频率
fp1=4500;fp2=6500; %第一截止频率
fs1=3000;fs2=7500; %第二截止频率
As=40; %最小阻带衰减
Ws1=(fp1+fs1)/fs;Ws2=(fp2+fs2)/fs; %截止频率归一化处理
w=(fp1-fs1)/fs; %求归一化过渡带
M=ceil((As-7.95)/(14.36*w)) %计算所需滤波器的阶数
juxing=boxcar(M+1); %生成长度为M+1的矩形窗
boxb=fir1(M,[Ws1,Ws2],juxing); %生成矩形窗设计的fir滤波器
freqz(boxb,1,fs,fs); %绘制幅频和相频响应曲线
运行结果：

图2-3矩形窗fir滤波器幅频和相频响应曲线
从幅频响应上看，通带基本无波纹，阻带中波纹较大，因而阻带较不理想，相频响应曲线在通带内为直线，效果较好，信号失真小。
2. 利用三角窗进行设计
利用三角窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用triang（）函数生成三角窗，程序运行结果如下：

图2-4三角窗设计的fir滤波器幅频和相频响应曲线
3. 利用汉明窗进行设计
利用汉明窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用hamming（）函数生成三角窗，程序运行结果如下：

图2-5汉明窗设计的fir滤波器幅频和相频响应曲线
4. 利用布莱克曼窗进行设计
利用布莱克曼窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用blackman（）函数生成三角窗，程序运行结果如下：

图2-6布莱克曼窗设计的fir滤波器幅频和相频响应曲线
5. 利用凯泽窗进行设计
利用凯泽窗进行设计时，滤波器的参数可以用以下函数求得：调用函数[n,wn,bta,ftype]=kaiserord（f，a，dev，fs）
参数： 1）f为对应的归一化频率
2）a为由f指定的各个频带上的幅值向量，一般只有0和1，和f长度关系为（2*a的长度）—2=（f的长度）
3）devs用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差，长度与a相等。
4）fs为信号的采样频率。
利用该函数修改代码得到凯泽窗设计fir滤波器的代码如下：
代码2-3
fs=20000; %设定采样频率
fp1=4500;fp2=6500; %第一截止频率
fs1=3000;fs2=7500; %第二截止频率
[n,wn,bta,ftype]=kaiserord([fs1,fp1,fp2,fs2],...
[0,1,0],[0.01 0.1087 0.01],fs) %求滤波器参数
b=fir1(n,wn,ftype,kaiser(n+1,bta)); %生成fir滤波器
freqz(b,1,fs,fs); %绘制幅频和相频响应曲线
程序运行结果：

图2-7凯泽窗设计的fir滤波器幅频和相频响应曲线
2.5.3利用MATLAB工具箱设计
MATLAB信号处理工具箱中已集成了用于滤波器设计和分析的工具：FDATool，利用它可以实现数字滤波器的可视化设计与分析，操作简单方便，在MATLAB命令行中输入fdatool命令即可打开滤波器设计工具，运行界面如图2-2所示：

图2-8 FDATool运行界面
通过选择滤波器的类型，设计方法，截止频率和其他一些相应参数后，点击Design即可得到滤波器的幅频响应曲线，在analysis菜单中可以选择查看相频响应，群时延，相时延等曲线，并且利用该工具箱设计的滤波器可以直接在simulink仿真中应用进行信号仿真，十分方便。
利用FDATool设计的凯瑟窗fir滤波器幅频响应和相频响应如下：

图2-9 FDATool设计的凯瑟窗fir滤波器幅频响应

图2-10 FDATool设计的凯瑟窗fir滤波器相频响应

3 滤波器测试
利用MATLAB中提供的filter函数可以选择不同的滤波器对数字信号进行滤波，这里主要对滤波器的滤波性能进行简单测试，这里主要针对滤波器的滤波性能和时延特性的测量。
3.1 滤波器滤波性能测试
由MATLAB模拟生成含有不同频率的数字信号，然后利用设计的滤波器对数字信号进行滤波，为方便观察，模拟生成的信号只含有包含在阻带的两个频率（2000Hz，8000Hz）一个包含于通带的频率（5000Hz），测试代码如下：
代码2-4：
fs=20000;
t=0:1/fs:2;
x=sin(2*pi*2000*t)+sin(2*pi*5000*t)...
+sin(2*pi*8000*t); %生成混合信号
xo=filter(b,2,xn); %用滤波器对信号进行滤波
figure;
nn=5000:5100; %取一段信号
subplot(211);
tt=nn/fs;
plot(tt,x(nn)); %绘制原始信号
axis([0.25,0.255,-4,4]);
ylabel('原始信号');xlabel('时间');
subplot(212);
plot(tt,xo(nn)); %绘制滤波后的信号
axis([0.25,0.255,-0.5,0.5]);
ylabel('滤波后的信号');xlabel('时间');
运行结果：

图3-1 滤波器滤波性能测试波形
从原始型号和滤波后的信号对比可以看出，在用设计的滤波器进行滤波后信号基本成单一频率的正弦波，滤波结果令人满意。
3.2 滤波器时延测量
滤波器的时延特性也是滤波器的性能指标之一，为了观测所设计滤波器的时延，可以利用一个还有一单一冲击的数字信号序列通过滤波器观测滤波器的输出波形，相对于滤波性能的测量只是在原始信号产生方法上有所不同，这里仅给出运行结果：

图3-2 滤波器时延性能测试波形
从图上看，滤波器仍然有一定的时延，这也是fir滤波器的缺点之一，并且时延与滤波器滤波性能相矛盾，滤波器的时延随滤波器阶数的升高而增加，实际设计时要综合考虑两方面的因素。
3.3 滤波器稳定性测量
FIR滤波器的一个突出优点便是它的稳定性，从信号与系统的理论可知，当一个系统的Z域的传递函数的极点都在单位圆内时系统是稳定的，FIR滤波器的传递函数的分母为1，即FIR滤波器一定是稳定的，MATLAB中也提供了专门用于绘制零极点图的函数：zplane()，通过该函数可以轻松绘制系统的零极点图，以凯泽窗设计的FIR滤波器为例，绘制其零极点图如下：

图3-3 凯泽窗FIR滤波器零极点图
从该零极点图可以看出，FIR设计的滤波器是属于无极点的系统，因而系统一定是稳定的。
从上面对滤波器三方面的测试可以看出FIR滤波器的特点，第一个便是其相位曲线，在通带内，相频曲线一定为直线，二是FIR滤波器由于阶数比较高，因而有一定的时延，三是FIR滤波器一定是稳定的。

4 心得体会
本次MATLAB课程设计主要任务是完成FIR滤波器的设计，对我来说这个题目还是很有挑战性的，因为自己对MATLAB中有关滤波器的设计知识了解较少，但正是这种有挑战性的题目才能提高自己的能力，才有研究价值，入手这个题目后我查阅了相关的资料，也从网上获得了不少有关MATLAB设计滤波器的资料，加上自己之间对MATLAB有一定的了解，因而设计思路渐渐明朗，经过自己的不断尝试和探索，终于弄明白了FIR滤波器的工作原理，通过查阅相关资料和研究MATLAB中提供的帮助信息，我也明白了FIR滤波器设计相关的一些函数的使用方法，并用它们来设计FIR滤波器，最终完成了题目。
在设计过程中，我也遇到了很多不懂得地方，程序经常出现错误，尤其是在利用所设计的滤波器对模拟出来的数字序列滤波时，出现很多错误，但经过自己的不断努力和尝试，最终还是解决了问题。
总之，通过本次课程设计，我收获很大，不只是学会了用MATLAB设计FIR滤波器，而是学会了自己学习新知识的一种方法。

5 参考文献
[1] 葛哲学，精通MATLAB.电子工业出版社,2008
[2] 陈亚勇，MATLAB信号处理详解.人民邮电出版社,2008
[3] 维基百科，http://zh.wikipedia.org
[4] 周开利，邓春辉，MATLAB基础及其应用教程.北京大学出版社,2007
[5] 赵静，张瑾，基于MATLAB的通信系统仿真.北京航空航天大学出版社,2006
[6] 宋寿鹏，数字滤波器设计及工程应用.江苏大学出版社,2007
[7] 普埃克著，方艳梅译数字信号处理(第四版).电子工业出版社,2007

附件一：
汉明窗带通滤波器设计源代码：
wp1=0.3*pi; ws1=0.45*pi;
wp2=0.65*pi; ws2=0.75*pi;
tr_width=ws1-wp1; %求过渡带宽度
M=ceil(6.6*pi/tr_width)+1; %求得所需窗函数的长度
n=[0:1:M-1];
wc1=(ws1+wp1)/2; wc2=(ws2+wp2)/2; %求截止频率
hd=ideal_lp(wc2,M)-ideal_lp(wc1,M); %求得理想带通的冲击响应
w_ham=(hamming(M))'; %得到长度为M的汉明窗
h=hd .* w_ham; %利用窗函数截短
%绘图部分
subplot(1,1,1)
subplot(2,21); stem(n,hd,'.'); title('理想冲击响应')
axis([0 M-1 -0.3 0.4]); xlabel('n'); ylabel('hd(n)')
subplot(2,2,2); stem(n,w_ham,'.');title('汉明窗')
axis([0 M-1 0 1.1]); xlabel('n'); ylabel('w(n)')
subplot(2,2,3); stem(n,h,'.');title('实际冲击响应')
axis([0 M-1 -0.3 0.4]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4); plot(w/pi,db);title('幅频响应');
axis([0 1 -100 10]); xlabel('f'); ylabel('dB')
理想冲击响应求解函数：
function hd = ideal_lp(wc,M);
alpha = (M-1)/2;
n = [0:1:(M-1)];
m = n - alpha + eps; %加入一个无穷小量，避免除零
hd = sin(wc*m) ./ (pi*m); %求得理想冲击响应

附件二：
本科生课程设计成绩评定表

姓 名

性 别

专业、班级

题 目：利用MATLAB仿真软件系统结合窗函数法设计一个数字带通FIR滤波器
答辩或质疑记录：

成绩评定依据：

最终评定成绩（以优、良、中、及格、不及格评定）

指导教师签字：

2010年 1月 08 日