根据窗函数法设计的数字带通FIR滤波器.doc

,目录目录1中文摘要21 窗函数设计法原理42 常见窗函数简介72.1 基本窗函数72.1.1 矩形窗函数72.1.2 三角窗函数72.1.3 巴特利特窗函数82.2 广义余弦窗92.2.1 汉宁窗函数92.2.2 海明窗函数102.2.3 布莱克曼窗函数112.2.4 凯 塞 窗112.2.5 切比雪夫窗123 方案设计与论证133.1 fdatool设计法133.2 程序设计法154 窗函数仿真结果分析174.1 矩形窗函数仿真结果174.2三角形窗函数仿真结果184.3 巴特利特窗函数仿真结果194.4 汉宁窗函数仿真结果204.5 海明窗函数仿真结果214.6布莱克曼窗函数仿真结果224.7 凯塞窗函数仿真结果234.8 切比雪夫窗函数仿真结果244.9 所有带通滤波器的比较255 总结与体会266参考文献27中文摘要现代图像、语声、数据通信对线性相位的要求是普遍的。正是此原因，使得具有线性相位的FIR数字滤波器得到大力发展和广泛应用。在实际进行数字信号处理时，往往需要把信号的观察时间限制在一定的时间间隔内，只需要选择一段时间信号对其进行分析。这样，取用有限个数据，即将信号数据截断的过程，就等于将信号进行加窗函数操作。而这样操作以后，常常会发生频谱分量从其正常频谱扩展开来的现象，即所谓的“频谱泄漏”。当进行离散傅立叶变换时，时域中的截断是必需的，因此泄漏效应也是离散傅立叶变换所固有的，必须进行抑制。而要对频谱泄漏进行抑制，可以通过窗函数加权抑制DFT的等效滤波器的振幅特性的副瓣，或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。而在后面的FIR滤波器的设计中，为获得有限长单位取样响应，需要用窗函数截断无限长单位取样响应序列。另外，在功率谱估计中也要遇到窗函数加权问题。由此可见，窗函数加权技术在数字信号处理中的重要地位。1 窗函数设计法原理数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（iir）滤波器和有限长冲激响应（fir）滤波器。iir数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。所以iir滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。fir数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。因此设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲击响应作为H(z)的系数，冲击响应长度N就是系统函数H(z)的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。要设计一个线性相位的FIR数字滤波器，首先要求理想频率响应。是w的周期函数，周期为，可以展开成傅氏级数： （公式1-1）其中是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计FIR DF用的h(n)，因为一般都是无限长、非因果的，物理上无法实现。为了设计出频响类似于理想频响的滤波器，可以考虑用h(n)来近似。窗函数的基本思想：先选取一个理想滤波器（它的单位抽样响应是非因果、无限长的），再截取（或加窗）它的单位抽样响应得到线性相位因果FIR滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器。设x(n)是一个长序列，w(n)是长度为N的窗函数，用w(n)截断x(n)，得到N点序列xn(n)，即xn(n) = x(n) w(n) （公式1-2）在频域上则有 （公式1-3）由此可见，窗函数w(n)不仅仅会影响原信号x(n)在时域上的波形，而且也会影响到频域内的形状。MATLAB信号工具箱主要提供了以下几种窗函数，如表1-1所示：表1-1 MATLAB窗函数窗窗 函 数矩形窗Boxcar巴特利特窗Barlett三角窗Triang布莱克曼窗Blackman海明窗Hamming汉宁窗Hanning凯塞窗Kaiser切比雪夫窗Chebwin加矩形窗后的频谱和理想频谱可得到以下结论：加窗使过渡带变宽，过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是。N越大，过渡带越窄、越陡；过渡带两旁产生肩峰，肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95，与N无关。工程上习惯用相对衰耗来描述滤波器，相对衰耗定义为 （公式1-4）这样两个肩峰点的相对衰耗分别是0.74dB和-21dB。其中（-0.0895）对应的点的值定义为阻带最小衰耗。以上的分析可见，滤波器的各种重要指标都是由窗函数决定，因此改进滤波器的关键在于改进窗函数。窗函数谱的两个最重要的指标是：主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为： 旁瓣峰值衰耗20lg(第一旁瓣峰值主瓣峰值) （公式1-5）为了改善滤波器的性能，需使窗函数谱满足：主瓣尽可能窄，以使设计出来的滤波器有较陡的过渡带；第一副瓣面积相对主瓣面积尽可能小，即能量尽可能集中在主瓣，外泄少，使设计出来的滤波器的肩峰和余振小。但上面两个条件是相互矛盾的，实际应用中，折衷处理，兼顾各项指标。2 常见窗函数简介2.1 基本窗函数数字信号处理领域中所用到的基本窗函数主要有：矩形窗函数、三角窗函数和巴特利特窗函数。下面就对这些窗函数展开介绍。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.1.3 巴特利特窗函数巴特利特窗函数的时域形式可以表示为当n为奇数时 （公式2-6）当n为偶数时 （公式2-7）Bartlett函数：生成巴特利特窗调用方式w = bartlett(n)：(1) 输入参数n是窗函数的长度；(2) 输出参数w是由窗函数的值组成的n阶向量。(3) 巴特利特窗是两个矩形窗的卷积。(4) 此函数生成的巴特利特窗的两端总是0。(5) 该函数与三角窗函数非常类似。当n是奇数时，则有bartlett(n+2)的中间部分等于triang(n)；当n是偶数时，巴特利特窗的傅立叶变换是负数。2.2 广义余弦窗汉宁窗、海明窗和布莱克曼窗，都可以用一种通用的形式表示，这就是广义余弦窗。这些窗都是广义余弦窗的特例，汉宁窗又被称为余弦平方窗或升余弦窗，海明窗又被称为改进的升余弦窗，而布莱克曼窗又被称为二阶升余弦窗。采用这些窗可以有效地降低旁瓣的高度，但是同时会增加主瓣的宽度。这些窗都是频率为0、2/(N1)和4/(N1)的余弦曲线的合成，其中为窗的长度。通常采用下面的命令来生成这些窗： （公式2-8） （公式2-9）其中，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 汉宁窗函数汉宁窗函数的时域形式可以表示为： （公式2-10）它的频域特性为： （公式2-11）其中，为矩形窗函数的幅度频率特性函数。汉宁窗函数的最大旁瓣值比主瓣值低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 海明窗函数海明窗函数的时域形式可以表示为 （公式2-12）它的频域特性为 （公式2-13）其中，为矩形窗函数的幅度频率特性函数。海明窗函数的最大旁瓣值比主瓣值低41dB，但它和汉宁窗函数的主瓣宽度是一样大的。Hamming函数：生成海明窗调用方式(1) w = hamming(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。(2) w = hamming(n,sflag)：参数sflag用来控制窗函数首尾的两个元素值；其取值为symmetric或periodic；默认值为symmetric。2.2.3 布莱克曼窗函数布莱克曼窗函数的时域形式可以表示为 (公式2-14)它的频域特性为 （公式2-15）其中，为矩形窗函数的幅度频率特性函数。布莱克曼窗函数的最大旁瓣值比主瓣值低57dB，但是主瓣宽度是矩形窗函数的主瓣宽度的3倍，为12/N。Blackman函数：生成海明窗调用方式(1) w = blackman (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。(2) w = blackman (n,sflag)：参数sflag用来控制窗函数首尾的两个元素值；其取值为symmetric或periodic；默认值为symmetric。2.2.4 凯 塞 窗上面所讨论的几种窗函数，在获得旁瓣抑制的同时却增加了主瓣的宽度。而凯塞窗定义了一组可调的窗函数，它是由零阶贝塞尔函数构成的，其主瓣能量和旁瓣能量的比例是近乎最大的。而且，这种窗函数可以在主瓣宽度和旁瓣高度之间自由选择它们的比重，使用户的设计变得非常灵活。凯塞窗函数的时域形式可表示为 (公式2-16)其中，是第1类变形零阶贝塞尔函数，是窗函数的形状参数，由下式确定： （公式2-17）其中，为凯塞窗函数的主瓣值和旁瓣值之间的差值(dB)。改变的取值，可以对主瓣宽度和旁瓣衰减进行自由选择。的值越大，窗函数频谱的旁瓣值就越小，而其主瓣宽度就越宽。Kaiser函数：生成凯塞窗调用方式w = kaiser(n,beta)：输入参数n是窗函数的长度；输入参数beta用于控制旁瓣的高度；输出参数w是由窗函数的值组成的n阶向量。n一定时，beta越大，其频谱的旁瓣就越小，但主瓣宽度相应的增加；当beta一定时，n发生变化，其旁瓣高度不会发生变化。2.2.5 切比雪夫窗对于给定的旁瓣高度，切比雪夫窗的主瓣宽度最小。这是因为它的旁瓣具有相同的高度，也就是具有等波纹性。切比雪夫窗在边沿的采样点有尖峰。Chebwin函数：生成切比雪夫窗调用方式w = chebwin(n,r)：输入参数n是窗函数的长度；输入参数r用于控制旁瓣的峰值低于主瓣的分贝数。3 方案设计与论证用窗函数法设计一个FIR带通滤波器。指示如下：低端阻带截止频率 wls = 0.2*pi; 低端通带截止频率 wlp = 0.35*pi; 高端通带截止频率 whp = 0.65*pi;高端阻带截止频率 whs = 0.8*pi3.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（巴特沃思）法、chebyshev type i（切比雪夫i型）法、 chebyshev type ii（切比雪夫ii型）法、elliptic（椭圆滤波器）法和fir滤波器的equiripple法、least-squares（最小乘方）法、window（窗函数）法。filter order（滤波器阶数）选项，定义滤波器的阶数，包括specify order（指定阶数）和minimum order（最小阶数）。在specify order中填入所要设计的滤波器的阶数（n阶滤波器，specify ordern-1），如果选择minimum order则matlab根据所选择的滤波器类型自动使用最小阶数。frenquency specifications选项，可以详细定义频带的各参数，包括采样频率fs和频带的截止频率。它的具体选项由filter type选项和design method选项决定，例如bandpass（带通）滤波器需要定义fstop1（下阻带截止频率）、fpass1（通带下限截止频率）、fpass2 （通带上限截止频率）、fstop2（上阻带截止频率），而lowpass（低通）滤波器只需要定义fstop1、fpass1。采用窗函数设计滤波器时，由于过渡带是由窗函数的类型和阶数所决定的，所以只需要定义通带截止频率，而不必定义阻带参数。magnitude specifications选项，可以定义幅值衰减的情况。例如设计带通滤波器时，可以定义wstop1（频率fstop1处的幅值衰减）、wpass （通带范围内的幅值衰减）、wstop2（频率fstop2处的幅值衰减）。当采用窗函数设计时，通带截止频率处的幅值衰减固定为6db，所以不必定义。window specifications选项，当选取采用窗函数设计时，该选项可定义，它包含了各种窗函数。3.2 程序设计法程序法只以Bartlett窗为示例，其他窗函数的程序大体相同，只是在window=bartlett(N)和hn = fir1(N-1,wc,bartlett(N)处要调自己的窗函数即可。Bartlett窗设计程序示例：wls = 0.2*pi; wlp = 0.35*pi; whp = 0.65*pi; wc = wlp/pi,whp/pi; B = wlp-wls; N = ceil(8/0.15); n=0:N-1; window=bartlett(N); h1,w=freqz(window,1); figure; stem(window); xlabel(n); title(bartlett窗函数); figure; plot(w/pi,20*log(abs(h1)/abs(h1(1); grid; xlabel(w/pi); ylabel(幅度(dB); title(bartlett窗函数的频谱); hn = fir1(N-1,wc,bartlett(N); h2,w=freqz(hn,1,512); figure stem(n,hn); xlabel(n); ylabel(h(n); title(bartlett窗函数的单位脉冲响应); figure; plot(w/pi,20*log(abs(h2)/abs(h2(1); grid; xlabel(w/pi); ylabel(幅度(dB); title(bartlett带通滤波器的幅度特性); 根据调试结果，2种设计法所出的结果没有什么大的区别，本次课程设计2种方法都使用，后面所附图片均来自这2种方法所处的仿真结果。4 窗函数仿真结果分析4.1 矩形窗函数仿真结果图 4-1 矩形窗函数波形图图4-2 矩形窗函数的频谱图图4-3 用矩形窗设计带通滤波器的幅频响应图图4-4 用矩形窗设计带通滤波器的相频响应图图4-5 用矩形窗设计带通滤波器的冲击响应图图4-6 用矩形窗设计带通滤波器的零极点图4.2三角形窗函数仿真结果图4-7 三角形窗函数波形图图4-8 三角形窗函数频谱图图4-9 用三角形窗设计带通滤波器的幅频响应图图4-10 用三角形设计带通滤波器的相频响应图图4-11 用三角形窗设计带通滤波器的冲击响应图图4-12 用三角形窗设计带通滤波器的零极点图4.3 巴特利特窗函数仿真结果图4-13 Bartlett窗函数波形图图4-14 Bartlett窗函数频谱图图4-15 Bartlett窗设计带通滤波器的幅频响应图图4-16 Bartlett窗设计带通滤波器的相频响应图图4-17 Bartlett窗设计带通滤波器的冲击响应图图4-18 Bartlett窗设计带通滤波器的冲击响应图4.4 汉宁窗函数仿真结果图4-19 汉宁窗函数波形图图4-20 汉宁窗函数频谱图图4-21 用汉宁窗设计带通滤波器的幅频响应图图4-22 用汉宁窗设计带通滤波器的相频响应图图4-23 用汉宁窗设计带通滤波器的冲击响应图图4-24 用汉宁窗设计带通滤波器的零极点图4.5 海明窗函数仿真结果图4-25 海明窗函数波形图图4-26 海明窗函数频谱图图4-27 用海明窗设计带通滤波器的幅频响应图图4-28 用海明窗设计带通滤波器的相频响应图图4-29 用海明窗设计带通滤波器的冲击响应图图4-30 用海明窗设计带通滤波器的幅频响应图4.6布莱克曼窗函数仿真结果图4-31 布莱克曼窗函数波形图图4-32 布莱克曼窗函频谱形图图4-33布莱克曼窗设计带通滤波器的幅频响应图图4-34 布莱克曼窗设计带通滤波器的相频响应图图4-35 布莱克曼窗设计带通滤波器的冲击响应图图4-36 用布莱克曼窗设计带通滤波器的零极点图4.7 凯塞窗函数仿真结果图4-37凯塞窗函数波形图图4-38 凯塞窗函数频谱图图4-39 用凯塞窗设计带通滤波器的幅频响应图图4-40 用凯塞窗设计带通滤波器的相频响应图图4-41 用凯塞窗设计带通滤波器的冲击响应图图4-42 用凯塞窗设计带通滤波器的零极点图4.8 切比雪夫窗函数仿真结果图4-43切比雪夫窗函数波形图图4-44切比雪夫窗函数频谱图图4-45 切比雪夫窗设计带通滤波器的幅频响应图图4-46用切比雪夫窗设计带通滤波器的相频响应图图4-47 切比雪夫窗设计带通滤波器的冲击响应图图4-48 用切比雪夫窗设计带通滤波器的零极点图4.9 所有带通滤波器的比较图4-49 用矩形窗设计带通滤波器的幅频响应图图4-50 用矩形窗设计带通滤波器的相频响应图图4-51 用三角形窗设计带通滤波器的幅频响应图图4-52 用三角形窗设计带通滤波器的相频响应图图4-53巴特利特窗设计带通滤波器的幅频响应图图4-54巴特利特窗设计带通滤波器的相频响应图图4-55汉宁窗设计带通滤波器的幅频响应图图4-56汉宁窗设计带通滤波器的相频响应图图4-57海明窗设计带通滤波器的幅频响应图图4-58海明窗设计带通滤波器的幅频响应图图4-59布莱克曼窗设计带通滤波器的幅频响应图图4-60布莱克曼窗设计带通滤波器的相频响应图图4-61凯塞窗设计带通滤波器的幅频响应图图4-62凯塞窗设计带通滤波器的相频响应图图4-63切比雪夫窗设计带通滤波器的幅频响应图图4-64切比雪夫窗设计带通滤波器的相频响应图5 总结与体会设计带通滤波器时首先要计算出过渡带，然后查表得到不同窗函数所需要的阶数，不同的窗函数所设计的滤波器的形状各有差异，尤其在主瓣宽度、旁瓣的形状以及主瓣与旁瓣的高度差上有比较明显得差别，实际应用中应根据实际情况，折衷处理，兼顾各项指标。为了这次课程设计，自己自学了数字信号处理领域中窗函数的有关知识。实际中遇到的离散时间信号总是有限长的，因此不可避免地要遇到数据截断问题。而在信号处理中，对离散序列的数据截断是通过序列与窗函数相乘来实现的。而且，有关滤波器的设计、功率谱估计等基本概念也要用到窗函数。本次课程设计对经常用到的下面8个窗函数：矩形窗函数、三角窗函数、巴特利特窗函数、汉宁窗函数、海明窗函函数、布莱克曼窗函数、凯塞窗函数以及切比雪夫窗函数，先是做了基本概念上的阐释，然后对其MATLAB实现函数做出了说明，最后又结合具体的实例，对这些窗函数的频域特性等进行了介绍。通过这次学习，我不但掌握了FIR数字滤波器窗函数的基本知识及其实际应用的技巧了，还提高了自己的编程和写报告的能力，收获颇多。6参考文献1. 陈怀琛MATLAB及在电子信息课程中的应用(第3版) 电子工业出版社20062. 刘泉 数字信号处理原理与实现 电子工业出版社 20053. 刘泉 信号与系统 高等教育出版社 2006