发布时间：2023-05-14 19:09:10
基于MATLAB的FIR数字高通滤波器分析和设计
雷学堂;徐火希
以FIR数字高通滤波器为例,详细分析时域卷积运算和频域加权算法的物理意义。并利用MATLAB的声音处理函数作为数据接口,利用多媒体播放器作为交互界面,利用MATLAB的FDA-Tool作为滤波器设计工具,设计一组语音高通滤波器,通过对比滤波前后的语音效果,可加深对数字信号处理的认识。
【作者单位】：黄冈师范学院物理科学与技术系 湖北黄冈438000
【关键词】：数字信号处理;数字滤波;高通滤波;MATLAB
【基金】：黄冈师范学院扰燃尘科学研究、青年科研基金资助项目(03CQ61)
【分类号】：TN713.4
【DOI】：cnki:ISSN:1009-3907.0.2006-10-009
【正文快照】：
0引言滤波就是有选择性地提取或去掉(或削弱)某一段或某几段频率范围内的信号,数字滤波器是一种用来过滤时间离散信号的数字系统,它是通过对抽样数据进行数学处理来达到选频目的。数字滤波器可分为IIR(无限冲激响应)和FIR(有限冲激响应)两种结构,FIR滤波器最大的优点是可设计成
数字信号处理（Digital Signal Processing，简称DSP）是一门涉及许多学科而又广泛应用于许多领域的新兴学科。20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。在过去的二十多年时间里，数字信号处理已经在通信等领域得到极为广泛的应用。
数字信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、滤波、估值、增强、压缩、识别等处理，以得到符合人们需要的信号形式。
数字信号处理是围绕着数字信号处理的理论、实现和应用等几个方面发展起来的。数字信号处理在理论上的发展推动了数字信号处理应用的发展。反过来，数字信号处理的应用又促进了数字信号处理理论的提高。而数字信号处理的实现则是理论和应用之间的桥梁。
虽然数字信号处理的理论发展迅速，但在20世纪80年代以前，由于实现方法的限制，数字信号处理的理论还得不到广泛的应用。直到20世纪70年代末80年代初世界上第一片单片可编程DSP芯片的诞生，才将理论研究结果广泛应用到低成本的实际系统中，并且推动了新的理论和应用领域的发展。可以毫不夸张地说，DSP芯片的诞生及发展对近20年来通信、计算机、控制等领域的技术发展起到十分重要的作用。
2 FIR数字滤波器的设计原理
2.1 数字滤波器
数字滤波器的可分为FIR(有限脉冲响应)和IIR(无限脉冲响应)两种。IIR滤波器的系统函数是两个Z的多项式的有理分式，而FIR滤波器的分母为1，即只有一个缓禅分子段睁多项式。
数字滤波器的理想幅频特性如图2.1所示。在0到 的全部字段上，其幅值为1的区域为通带。其余为阻带，其幅值为0。根据Wc1和Wc2取值不同可以分为四种类型：
朋友，找了半个小时总算有了成果，希望对你有帮助，谢谢。
1.FIR数字滤波器原理
假设理想低通滤波器的截止频率为ωc＝2πfc，且具有线性相位，群延时为a，即频率响应：
航空重力勘探理论方法及应用
表示成幅度函数和相位函数形式：
航空重力勘探理论方法及应用
则幅度函数：
航空重力勘探理论方法及应用
在通带范围｜ω｜ ≤ωc（截止频率）内Hd（ejω）的幅度为1，相位为－ωα；对应的时间域（或空间域）滤波函数为：
航空重力勘探理论方法及应用
有限脉冲响应FIR（Finite Impulse Response）数字滤波器要求用有限长的单位冲击响应h（n）来逼近无限长的理想滤波器的单位冲击响应hd（n），最常用和有效的方法就是用一个有限长（长度为N）的“窗函数”序列w（n）来截取hd（n）的主要成分（陈玉东，2005）：
航空重力勘探理论方法及应用
实际上是用有限长的h（n）去逼近hd（n），通过这种方式得到的频率响应H（ejω）近似于理想频率响应Hd（ejω）（在频率域内采用均方差最小准则逼近）。按照线性相位滤波器的约束要求，h（n）必须是偶对称的，其对称中心应为它长度的一半：h（n）=h（N-1-n），而且
；所以同时要求窗函数w（n）也必须是关于中心偶对称：w（n）=w（N-1-n）。
2.几种常见窗函数
（1）矩形窗
长度为N的矩形窗函数为：
航空重力勘探理论方法及应用
（2）三角形窗（Bartlett）
长度为N的三角形窗函数为：
航空重力勘探理论方法及应用
（3）汉宁窗（Hanning）
长度为N的汉宁窗函数为：
航空重力勘探理论方法及应用
（4）海明窗（Hamming）
为使得旁瓣更小，可将汉宁窗改进成海明窗，长度为N的海明窗函数为：
航空重力勘探理论方法及应用
（5）布拉克曼窗（Blackman）
为进一步有效抑制旁瓣，可以再加上余弦的二次谐波分量，得到长度为N的布拉克曼窗函数为：
航空重力勘探理论方法及应用
（6）凯泽窗（Kaiser）
长度为N的凯泽窗函数为：
航空重力勘探理论方法及应用
其中I0（x）为第一类变形零阶贝塞尔函数，α＝（N-1）/2。β是一个可自由选择的参数，它可同时调整窗函数谱主瓣宽度与旁瓣幅值；β越大，则窗函数w（n）变化越快、变得越窄拦蚂，频谱旁瓣就越小，但主瓣宽度相应增加。一般选择4＜β＜9，相当于窗函数频谱旁瓣幅度与主瓣幅度的比值由3.1%变到0.047%。β＝0时相当于矩形窗（陈玉东，2005）。
3.窗函数FIR滤波器
式（7-4-3）至式（7-4-8）窗函数都满足关于中心偶对称的线性简猜埋相位滤波器的约束要求，结合式（7-4-1）至式（7-4-2）可以得到相应窗函数的FIR低通数字滤波器函数（郭志宏，罗锋，等，2007）:
航空重力勘探理论方法及应用
航空重力勘探理论方法及应用
用该滤波器窗口对时间域（或空间域）长度为M的数据序列逐点进行窗口滑动卷积求和计算（实际处理时窗口中点作为输出计算点，则一边损失半个滤波窗口数据），就可获得FIR滤波后的数据（郭志宏，段树岭，等，2009）：
航空重力勘探理论方法及应用
h（n）为滤波器系数，x（n）、y（n）分别为输入、输出数据序列。
4.窗函数FIR滤波试验
（1）GT-1A型航空重力数据
图7-4-1至图7-4-2分别为GT-1A型航空重力系统获得的一条原始未滤波、100 s和60 s滤波自由空间重力异常测线数据，其中飞机的飞行速度约60m/s，剖面图横轴为测线基准点号，基准点间距约30 m。图7-4-1中GT-1A型系统航空原始未滤波自由空间重力测线数据的高频干扰非常之严重，噪声幅度在－5 000×10-5m·s-2至5000×10-5m·s-2的大范围内变化，而幅度通常只有（10-3～10-4）m·s-2的由密度和构造变化等地质因素引起的重力异常信号（图7-4-2）则完全淹没在高频干扰中。图7-4-2中GT-1A型系统航兆仿空100 s、60 s滤波自由空间重力测线数据是采用GT-1A型航空重力系统自带软件模块由图7-4-1的航空原始未滤波自由空间重力测线数据获得的滤波数据，滤波后高频干扰已基本消除，油气和矿产地球物理勘查所需的重力异常则较好的显现出来。
图7-4-1 GT-1A型航空重力系统原始未滤波自由空间重力异常
（2）几种窗函数FIR滤波试验
根据式（7-4-1）至式（7-4-10），我们研制了窗函数法FIR数字滤波计算软件，用各种窗函数FIR滤波器对图7-4-1的GT-1A航空原始未滤波自由空间重力测线数据分别进行了截止波长为100 s、60 s长度（按v=60m/s的航速计算，截止波长A。分别为6km、3.6km，按fc=v/λc计算的截止频率分别为0.01 Hz、0.0167 Hz）的低通滤波试验计算，试验结果见图7-4-3至图7-4-8。为了图形对比方便，各剖面图中仍然保留了测线边部两端的半个滤波窗口数据，这些数据由于存在边部效应，因而是不准确的，实际应用时应该去掉。从试验结果图可以看到，矩形窗和三角窗FIR滤波后异常整体形状虽然也与图7-4-2类似，但其上叠加了高频扰动，尤其是矩形窗FIR滤波结果，这就是通常所说的“吉布斯”振荡效应（陈玉东，2005）。如果在图7-4-3至图7-4-4的基础上，采用空间域非线性曲率滤波方法（郭志宏，刘浩军，等，2003），用中国国土资源航空物探遥感中心的“空中探针”系统（刘浩军，薛典军，等，2003）中的滤波软件进一步处理，则可获得消除扰动后接近图7-4-2效果的异常数据。从汉宁窗、海明窗、布拉克曼窗以及凯泽窗FIR滤波试验结果看到，通过选择合适的窗口长度、截至波长等滤波参数，基本都获得了令人满意的效果。
图7-4-2 GT-1A型航空重力系统100s、60s滤波自由空间重力异常
图7-4-3 矩形窗FIR低通滤波截止波长100s、60s航空自由空间重力异常
表7-4-1为图7-4-3至图7-4-8所示的各种窗函数FIR低通滤波截止波长100 s、60 s长度航空自由空间重力异常与图7-4-2所示的GT-1A型航空重力系统100 s、60 s滤波自由空间重力异常（作为标准）的比较，通过两者之差值的统计结果来衡量吻合程度。从统计表中可以看到，除了矩形窗、三角窗外，其他几种窗函数FIR低通滤波结果的差异值都在±1×10-5m·s-2以内，均方差值则多数为0.3×10-5m·s-2左右，可见吻合程度还是比较好的。
图7-4-4 三角窗FIR低通滤波截止波长100s、60s航空自由空间重力异常
5.结论
1）通过选择合适的窗形、窗口长度、滤波参数，窗函数法FIR低通数字滤波器可以在航空重力数据的滤波处理中发挥应有的作用。
图7-4-5 汉宁窗FIR低通滤波截止波长100s、60s航空自由空间重力异常
图7-4-6 海明窗FIR低通滤波截止波长100s、60s航空自由空间重力异常
图7-4-7 布拉克曼窗FIR低通滤波截止波长100s、60s航空自由空间重力异常
图7-4-8 凯泽窗（β＝6）FIR低通滤波截止波长100s、60s航空自由空间重力异常
2）为了获得与GT-1A型航空重力系统100 s、60 s低通滤波（60m/s航速）对应的自由空间重力测线数据，所选择汉宁、海明、布拉克曼、凯泽窗的长度通常为400点（2 Hz采样率），FIR低通滤波对应的截止频率分别为0.01 Hz、0.0167 Hz。
3）窗函数法不但可以设计FIR低通滤波器，还可设计FIR高通、带通、带阻滤波器等。通常一个高通滤波器相当于一个全通滤波器减去一个低通滤波器；一个带通滤波器相当于两个低通滤波器相减；而一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器。
表7-4-1 窗函数FIR滤波试验结果与GT-1A系统滤波结果的差值统计
4）除了窗函数法FIR低通滤波器，其他诸如等波纹法FIR低通滤波器、无限脉冲响应IIR低通滤波、Kalman滤波等方法（周坚鑫，刘浩军，等，2001；陈玉东，2005）均可用于航空重力数据的低通数字滤波处理中。
我举个例子好了。蠢锋升
矩形窗的窗函数是w=boxcar(n)；
其中n是窗长度。
输入这三行代码就可以看到矩形窗的频率响应了。
n=100;
w=boxcar(n);
fvtool(w)；
如果你想看其他窗函数的频率响应，把boxcar换掉就可以了。
（1）矩形窗基橡（rectangle
window）
调用格式：w=boxcar(n)，根据长度
n
产生一个矩形窗
w。
（2）三角窗（triangular
window）
调用格式：w=triang(n)，根据长度
n
产生一个三角窗
w。
（3）汉宁窗（hanning
window）
调用格式：w=hanning(n)，根据长度
n
产生一个汉宁窗
w。
（4）海明窗（hamming
window）
调用格式：w=hamming(n)，根据长度
n
产生一个海明窗
w。
（5）布拉克曼窗（blackman
window）
调用格式：w=blackman(n)，根据长度
n
产生一个布拉克曼窗
w。
（6）恺撒窗（kaiser
window）
调用格式：w=kaiser(n,beta)，根据长度
n
和影响窗函数旁瓣的β参数产生一个带老恺撒窗w。
参考的函数就是这几个
我举个例子好了。 矩形窗的窗函数是w=boxcar(n)； 其中n是窗长度。
输入这三行代码就可以看到矩形窗的频率响应了。
n=100;
w=boxcar(n);
fvtool(w)；
如果你想看其他窗函数的频率响应槐烂御，把boxcar换掉就可以了。
（1）矩形窗（Rectangle Window） 调用格式：w=boxcar(n)，根据长度 n 产生一个矩形窗 w。
（2）历空三角窗（Triangular Window） 调用格式：w=triang(n)，根据长度 n 产生一个三角窗 w。
（3）汉宁窗（Hanning Window） 调用格式：w=hanning(n)，根据长度 n 产生一个汉宁窗 w。
（4）海明窗（Hamming Window） 调用格式：w=hamming(n)，根据长度 n 产生一个海明窗 w。
（5）布拉克曼窗（Blackman Window） 调用格式：w=blackman(n)，根据长度 n 产生一个布拉克曼窗 w。
（6）恺撒窗（Kaiser Window） 调用格式：w=kaiser(n,beta)，根据长度 n 和铅岩影响窗函数旁瓣的β参数产生一个恺撒窗w。
