MATLAB基于窗函数设计的带通滤波器.doc

. .
课程设计任务书

学生XX： 专业班级： 通信0705
指导教师： 工作单位： 信息工程学院
题 目:基于窗函数法设计的数字带通FIR滤波器

初始条件：
① MATLAB软件
‚数字信号处理与图像处理根底知识

要求完成的主要任务:
利用MATLAB仿真软件系统结合窗函数法设计的数字带通FIR滤波器在数字信号处理平台上〔PC机﹑MATLAB仿真软件系统和TC++编程环境〕进展软件仿真设计，并进展调试和数据分析。

时间安排：
第18周理论设计、实验室安装调试，地点： 鉴主15楼通信实验室一

指导教师签名： 2021 年 1月 8日

系主任〔或责任教师〕签名： 年 月 日

. .word..
. .
目录

摘要I
AbstractII
1 窗函数设计法原理1
2 常见窗函数简介2
2.1 根本窗函数2
下面就几种常用的窗函数展开介绍。3
2.1.1 矩形窗函数3
2.1.2 三角窗函数3
2.2 广义余弦窗4
2.2.1 汉宁窗函数5
2.2.2 海明窗函数5
3 方案设计与论证7
3.1 fdatool设计法7
3.2 程序设计法8
4 窗函数仿真结果分析10
4.1 矩形窗函数仿真结果10
4.2三角形窗函数仿真结果11
4.3 汉宁窗函数仿真结果12
4.4海明窗函数仿真结果13
5 总结与体会14
6参考文献16

. .word..
. .
摘要
现代图像、语音、数据通信对线性相位的要求是普遍的。正是此原因，使得具有线性相位的FIR数字滤波器得到大力开展和广泛应用。
在实际进展数字信号处理时，往往需要把信号的观察时间限制在一定的时间间隔内，只需要选择一段时间信号对其进展分析。这样，取用有限个数据，即将信号数据截断的过程，就等于将信号进展加窗函数操作。而这样操作以后，常常会发生频谱分量从其正常频谱扩展开来的现象，即所谓的“频谱泄漏〞。当进展离散傅立叶变换时，时域中的截断是必需的，因此泄漏效应也是离散傅立叶变换所固有的，必须进展抑制。而要对频谱泄漏进展抑制，可以通过窗函数加权抑制DFT的等效滤波器的振幅特性的副瓣，或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。而在后面的FIR滤波器的设计中，为获得有限长单位取样响应，需要用窗函数截断无限长单位取样响应序列。另外，在功率谱估计中也要遇到窗函数加权问题。由此可见，窗函数加权技术在数字信号处理中的重要地位。

. .word..
. .
Abstract
A modern image, Voice and data munications to the linear phase requirement is universal. It is for this reason, making a linear phase FIR digital filters are to develop and widely used.
In practice, digital signal processing, they often need to observe the time limit for the signal in a certain time interval, only need to select a period of time to analyze the signals. In this way, access to a finite number of data about the process of signal data truncation is equivalent to the signal, the additional window function operation. And this after the operation, often occur from the normal ponent of the spectrum spread spectrum open to the phenomenon of so-called "spectral leakage." When it es to discrete Fourier transform, the time domain truncation is necessary, therefore, leakage effects are also inherent in discrete Fourier transform, and must be inhibited. But would like to suppress the spectral leakage can be inhibited by the weighted window function equivalent DFT filter amplitude characteristics of side-lobe, or the weighted window function so that the limited length of the extension of the input signal cycle, after the border to minimize the degree of discontinuity The method of implementation. In the back of the FIR filter design, in order to obtain finite sampling units to respond to the need to truncate the infinite length window function with unit sampling response sequence. In addition, power spectrum estimation problem also encountered in the weighted window function. Thus, the weighted window function in digital signal processing in an important position.
. .word..
. .
1 窗函数设计法原理
数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应〔iir〕滤波器和有限长冲激响应〔fir〕滤波器。fir数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最正确逼近法等。
用窗函数设计滤波器首先要对滤波器提出性能指标。一般是给定一个理想的频率响应，使所设计的FIR滤波器的频率响应去逼近所要求的理想的滤波器的响应。窗函数法设计的任务在于寻找一个可实现〔有限长单位脉冲响应〕的传递函数。用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列。主要设计步骤为：
〔1〕通过傅里叶逆变换获得理想滤波器的单位脉冲响应hd(n)。
从时域出发，截取有限长的一段冲击响应作为H(z)的系数，冲击响应长度N就是系统函数H(z)的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。要设计一个线性相位的FIR数字滤波器，首先要求理想频率响应。是w的周期函数，周期为，可以展开成傅氏级数：
＝ 〔公式1-1〕
其中是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计FIR DF用的h(n)，因为一般都是无限长、非因果的，物理上无法实现。为了设计出频响类似于理想频响的滤波器，可以考虑用h(n)来近似。
窗函数的根本思想：先选取一个理想滤波器〔它的单位抽样响应是非因果、无限长的〕，再截取〔或加窗〕它的单位抽样响应得到线性相位因果FIR滤波器。这种方法的重点是选择一个适宜的窗函数和理想滤波器。
〔2〕由性能指标确定窗函数W(n)和窗口长度N。
设x(n)是一个长序列，w(n)是长度为N的窗函数，用w(n)截断x(n)，得到N点序列xn(n)，即
xn(n) = x(n) w(n)〔公式1-2〕
在频域上那么有
〔公式1-3〕
(3) 求得实际滤波器的单位脉冲响应h(n)， h(n)即为所设计FIR滤波器系数向量。〔公式1-4〕
由此可见，窗函数w(n)不仅仅会影响原信号x(n)在时域上的波形，而且也会影响到频域内的形状。
2 常见窗函数简介
2.1 根本窗函数
MATLAB信号工具箱主要提供了以下几种窗函数，如表2-1所示：
表2-1 MATLAB窗函数
窗 函 数
定义式
矩形窗〔Boxcar〕

三角窗(Triang)

布莱克曼窗(Blackman)
0nM-1
海明窗(Hamming)
0nM-1
汉宁窗(Hanning)
0nM-1
巴特利特窗〔Bartlett〕

下面就几种常用的窗函数展开介绍。
2.1.1 矩形窗函数
矩形窗(Rectangular Window)函数的时域形式可以表示为：
〔公式2-1〕
它的频域特性为
〔公式2-2〕
Boxcar函数：生成矩形窗
调用方式w = boxcar (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。从功能上讲，该函数又等价于w = ones(n,1)。
2.1.2 三角窗函数
三角窗是最简单的频谱函数为非负的一种窗函数。三角窗函数的时域形式可以表示为：
当n为奇数时
〔公式2-3〕
当n为偶数时
〔公式2-4〕
它的频域特性为：
〔公式2-5〕
三角窗函数的主瓣宽度为，比矩形窗函数的主瓣宽度增加了一倍，但是它的旁瓣宽度却小得多。
Triang函数：生成三角窗
调用方式w = triang(n)；输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。
三角窗也是两个矩形窗的卷积。三角窗函数的首尾两个数值通常是不为零的。当n是偶数时，三角窗的傅立叶变换总是非负数。

2.2 广义余弦窗

汉宁窗、海明窗和布莱克曼窗，都可以用一种通用的形式表示，这就是广义余弦窗。这些窗都是广义余弦窗的特例，汉宁窗又被称为余弦平方窗或升余弦窗，海明窗又被称为改良的升余弦窗，而布莱克曼窗又被称为二阶升余弦窗。采用这些窗可以有效地降低旁瓣的高度，但是同时会增加主瓣的宽度。
这些窗都是频率为0、2π/(N–1)和4π/(N–1)的余弦曲线的合成，其中为窗的长度。通常采用下面的命令来生成这些窗：
〔公式2-6〕
〔公式2-7〕
其中，A、B、C适用于自己定义的常数。根据它们取值的不同，可以形成不同的窗函数，分别是：
汉宁窗 A=0.5，B=0.5，C=0；
海明窗 A=0.54，B=0.54，C=0；
布莱克曼窗 A=0.5，B=0.5，C=0.08；
2.2.1 汉宁窗函数
汉宁窗函数的时域形式可以表示为：
0nM-1 〔公式2-8〕
它的频域特性为：
〔公式2-9〕
其中，为矩形窗函数的幅度频率特性函数。
汉宁窗函数的最大旁瓣值比主瓣值低31dB，但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍，为8π/N。
hanning函数：生成汉宁窗
调用方式
(1) w = hanning(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。
注意：此函数不返回是零点的窗函数的首尾两个元素。
(2) w = hanning(n,'symmetric')：与上面相类似。
(3) w = hanning(n,'periodic')：此函数返回包括为零点的窗函数的首尾两个元素。

2.2.2 海明窗函数
海明窗函数的时域形式可以表示为
0nM-1 〔公式2-10〕
它的频域特性为
〔公式2-11〕
其中，为矩形窗函数的幅度频率特性函数。
海明窗函数的最大旁瓣值比主瓣值低41dB，但它和汉宁窗函数的主瓣宽度是一样大的。
Hamming函数：生成海明窗
调用方式
(1)w = hamming(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。
(2) w = hamming(n,sflag)：参数sflag用来控制窗函数首尾的两个元素值；其取值为symmetric或periodic；默认值为symmetric。

加矩形窗后的频谱和理想频谱可得到以下结论：
加窗使过渡带变宽，过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是。N越大，过渡带越窄、越陡；
过渡带两旁产生肩峰，肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95％，与N无关。工程上习惯用相对衰耗来描述滤波器，相对衰耗定义为
〔公式2-12〕
这样两个肩峰点的相对衰耗分别是0.74dB和-21dB。其中〔-0.0895〕对应的点的值定义为阻带最小衰耗。
以上的分析可见，滤波器的各种重要指标都是由窗函数决定，因此改良滤波器的关键在于改良窗函数。
窗函数谱的两个最重要的指标是：主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为：
旁瓣峰值衰耗＝20lg(第一旁瓣峰值／主瓣峰值) 〔公式2-13〕
为了改善滤波器的性能，需使窗函数谱满足：
1. 较低的旁瓣幅度，尤其是第一旁瓣；
2. 旁瓣幅度要下降得快，以利于增加阻带衰减；
3. 主瓣宽度要窄，这样滤波器过渡带较窄。
但这三点难以同时满足，中选用主瓣宽度较窄时，虽然得到的幅频特性较陡峭，但通带、阻带波动会明显增加；中选用较低的旁瓣幅度时，虽然得到的幅频特性较平缓匀滑，但过渡带变宽。因此，实际的选择往往是取折衷。
. .word..
. .
3 方案设计与论证
用窗函数法设计一个FIR带通滤波器。指标如下：
低端阻带截止频率 wls = 0.2*pi;
低端通带截止频率 wlp = 0.35*pi;
高端通带截止频率 whp = 0.65*pi;
高端阻带截止频率 whs = 0.8*pi
3.1 fdatool设计法
fdatool〔filter design & analysis tool〕是matlab信号处理工具箱里专用的滤波器设计分析工具，matlab6.0以上的版本还专门增加了滤波器设计工具箱〔filter design toolbox〕。fdatool可以设计几乎所有的根本的常规滤波器，包括fir和iir的各种设计方法。它操作简单，方便灵活。

图3-1 fdatool设计界面
fdatool界面总共分两大局部，一局部是design filter，在界面的下半部，用来设置滤波器的设计参数，另一局部那么是特性区，在界面的上半局部，用来显示滤波器的各种特性。design filter局部主要分为：
filter type〔滤波器类型〕选项，包括lowpass〔低通〕、highpass〔高通〕、bandpass〔带通〕、bandstop〔带阻〕和特殊的fir滤波器。
design method〔设计方法〕选项，包括iir滤波器的butterworth〔巴特沃思〕法、elliptic〔椭圆滤波器〕法和fir滤波器的equiripple法、least-squares〔最小乘方〕法、window〔窗函数〕法。
filter order〔滤波器阶数〕选项，定义滤波器的阶数，包括specify order〔指定阶数〕和minimum order〔最小阶数〕。在specify order中填入所要设计的滤波器的阶数〔n阶滤波器，specify order＝n-1〕，如果选择minimum order那么matlab根据所选择的滤波器类型自动使用最小阶数。
frenquency specifications选项，可以详细定义频带的各参数，包括采样频率fs和频带的截止频率。它的具体选项由filter type选项和design method选项决定，例如bandpass〔带通〕滤波器需要定义fstop1〔下阻带截止频率〕、fpass1〔通带下限截止频率〕、fpass2 〔通带上限截止频率〕、fstop2〔上阻带截止频率〕，而lowpass〔低通〕滤波器只需要定义fstop1、fpass1。采用窗函数设计滤波器时，由于过渡带是由窗函数的类型和阶数所决定的，所以只需要定义通带截止频率，而不必定义阻带参数。
window specifications选项，中选取采用窗函数设计时，该选项可定义，它包含了各种窗函数。

3.2 程序设计法
程序法只以Boxcar窗为例如，其他窗函数的程序大体一样，只是在window=Boxcar(N)和hn = fir1(N-1,wc,Boxca(N))处要调自己的窗函数即可。
Boxcar窗设计程序例如：
wls = 0.2*pi; %低端阻带截止频率
wlp = 0.35*pi; %低端通带截止频率
whp = 0.65*pi; %高端阻带截止频率
wc = [wlp/pi,whp/pi]; %通带宽度，求hn时会用
B = wlp-wls; %过渡带
N = ceil(4/0.15); %求N的值，不同的窗函数前面的值不一样
n=0:N-1; %N阶向量

%以下2段是看窗函数以及其频谱的，是窗函数方面的根本情况

window=boxcar(N); %调用窗函数
[h1,w]=freqz(window,1); %求频率响应

figure; %新建一个图形框
stem(window); %离散画图，看窗函数
xlabel('n'); %X轴标签
title('矩形窗函数'); %标题

figure; %新建一个图形框
plot(w/pi,20*log(abs(h1)/abs(h1(1)))); %画频谱图，进展归一化
grid; %显示表格线
xlabel('w/pi'); %X轴标签，归一化
ylabel('幅度(dB)'); %Y轴标签
title('矩形窗函数的频谱'); %标题

%以下2段是关于带通方面的

hn = fir1(N-1,wc,boxcar(N)); %fir1标准频率响应形状，包括低通、带通、高通或带阻特性
[h2,w]=freqz(hn,1,512); %求离散的频响

figure %新建一个图形框
stem(n,hn); %离散画图，带通滤波器的单位脉冲响应
xlabel('n'); %X轴标签
ylabel('h(n)'); %Y轴标签
title('矩形窗函数带通滤波器的单位脉冲响应'); %标题

figure; %新建一个图形框
plot(w/pi,20*log(abs(h2)/abs(h2(1)))); %画带通频谱图，进展归一化
grid; %显示表格线
xlabel('w/pi'); %X轴标签
ylabel('幅度(dB)'); %Y轴标签
title('矩形窗函数带通滤波器的幅度特性'); %标题
根据调试结果，2种设计法所出的结果没有什么大的区别，本次课程设计2种方法都使用，后面所附图片均来自这2种方法所处的仿真结果。

. .word..
. .
4 窗函数仿真结果分析
4.1 矩形窗函数仿真结果

图 4-1 矩形窗函数波形图

图4-2 矩形窗函数的频谱图

图4-3 用矩形窗设计带通滤波器的幅频响应图

图4-4 用矩形窗设计带通滤波器的冲击响应图
4.2三角形窗函数仿真结果

图4-5 三角形窗函数波形图

图4-6 三角形窗函数频谱图

图4-7 用三角形窗设计带通滤波器的幅频响应图

图4-8 用三角形窗设计带通滤波器的冲击响应图
4.3 汉宁窗函数仿真结果

图4-9汉宁窗函数波形图

图4-10 汉宁窗函数频谱图

图4-11 用汉宁窗设计带通滤波器的幅频响应图

图4-12 用汉宁窗设计带通滤波器的冲击响应图

4.4海明窗函数仿真结果

图4-13海明窗函数波形图

图4-14 海明窗函数频谱图

图4-15 用海明窗设计带通滤波器的幅频响应图

图4-16 用海明窗设计带通滤波器的冲击响应图
5 总结与体会
设计带通滤波器时首先要计算出过渡带，然后查表得到不同窗函数所需要的阶数，不同的窗函数所设计的滤波器的形状各有差异，尤其在主瓣宽度、旁瓣的形状以及主瓣与旁瓣的高度差上有比拟明显得差异，实际应用中应根据实际情况，折衷处理，兼顾各项指标。
为了这次课程设计，自己自学了数字信号处理领域中窗函数的有关知识。实际中遇到的离散时间信号总是有限长的，因此不可防止地要遇到数据截断问题。而在信号处理中，对离散序列的数据截断是通过序列与窗函数相乘来实现的。而且，有关滤波器的设计、功率谱估计等根本概念也要用到窗函数。本次课程设计对经常用到的下面6个窗函数：矩形窗函数、三角窗函数、汉宁窗函数、海明窗函函数、布莱克曼窗函数以及凯塞窗函数，先是做了根本概念上的阐释，然后对其MATLAB实现函数做出了说明，最后又结合具体的实例，对这些窗函数的频域特性等进展了介绍。
通过这次学习，我不但掌握了FIR数字滤波器窗函数的根本知识及其实际应用的技巧了，还提高了自己的编程和写报告的能力，收获颇多。

. .word..
. .
6参考文献
[1]陈怀琛·MATLAB在电子信息课程中的应用〔第二版〕·电子工业·2006
[2] 郭仕剑·MATLAB7.X数字信号处理·人民邮电·2006年
[3] 陈贵明·用MATLAB语言处理数字信号与数字图像·科学·2000年
[4] 王家文·MATLAB7.0图形图像处理·国防工业·2006年
[5] 苏金明·MATLAB图形图像·电子工业·2005

本科生数字信号处理成绩评定表

XX

性别

专业、班级

题目：
辩论或质疑记录：

成绩评定依据：

最终评定成绩〔以优、良、中、及格、不及格评定〕

指导教师签字：

年月日
. .word..