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

摘要

随着信息时代的到来，数字信号处理已经成为当今一门及其重要的学科和技术，并且在通信、语音、图像、自动控制等众多领域得到了广泛的应用，在数字信号处理中，数字滤波器占有及其重要的地位，它具有精度高，可靠性好、灵活性大等特点，现代数字滤波器可以用软件和硬件两种方式来实现，软件方式实现的优点是可以通过滤波器参数的改变去调整滤波器的性能。

　　MATLAB是一种面向科学和工程计算的语言，它集数值分析，矩阵运算、信号处理和图形显示于一体，具有编程效率高、调试手段丰富、扩充能力强等特点，MATLAB的信号处理工具箱具有强大的函数功能，它不仅可以用来设计数字滤波器，还可以使设计达到最优化，是数字滤波器设计的强有力工具。本文介绍了如何利用MATLAB仿真软件系统结合窗函数法设计一个数字窗函数法

关键字：数字信号处理 数字滤波器 MATLAB 重要的地位 窗函数法窗函数法

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

1.1 数字滤波技术
数字滤波，就是通过一定的计算或判断程序减少干扰在有用信号中的比重，故实质上是一种程序滤波。与此对应的就是模拟滤波，由于模拟滤波牵扯到的其他知识太多在此不详细介绍了，模拟滤波主要无源绿波（直接用电阻、电容、电感等不外接电源的元件组成的）与有源滤波（如运算放大器等需要外接电源组成的），其目的是将信号中的噪音和干扰滤去或者将希望得到的频率信号滤出为我所用。数字滤波的出现克服了模拟滤波的很多不足，具有以下优点：
A．是用程序实现的，不需要增加硬设备，所以可靠性高，稳定性好。
B．可以对频率很低的信号实现滤波，克服了模拟滤波的缺陷。
C．可以根据信号的不同，采用不同的滤波方法或参数，具有灵活、方便、功能强的特点。
几种常用的滤波方法：
1． 算术平均值法
2． 中值滤波法
3.滑动平均值法
4.限幅滤波法
5.惯性滤波法
数字滤波技术通过数字滤波器实现，从实现方法上可以分为FIR数字滤波器和IIR数字滤波器，按功能可分为低通滤波器（LPF）、高通滤波器（HPF）、带通滤波器（BPF）和带阻滤波器（BSF）。本文主要对FIR滤波器加以介绍。
1.2 FIR滤波器
FIR（Finite Impulse Response）滤波器，即有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位冲激响应是有限的，没有输入到输出的反馈，是稳定的系统。因此，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）窗、汉明（Hanming）窗、布莱克曼（Blackman）窗、切比雪夫（Chebyshev）窗、凯泽（Kaiser）窗。
1.4 MATLAB结合窗函数设计法原理
数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。FIR数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。FIR滤波器具有严格的相位特性，这对于语音信号处理和数据传输是很重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率取样法和切比雪夫等波纹逼近的最优化设计方法。常用的是窗函数法和切比雪夫等波纹逼近的最优化设计方法。
因此设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲击响应作为H(z)的系数，冲击响应长度N就是系统函数H(z)的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。
要设计一个线性相位的FIR数字滤波器，首先要求理想频率响应。是w的周期函数，周期为，可以展开成傅氏级数：
（1-1）
使用上述的传递函数去逼近，一个理想的频率响应的傅立叶反变换：
（1-2）
其中是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计FIR DF用的h(n)，因为一般都是无限长、非因果的，物理上无法实现。为了设计出频响类似于理想频响的滤波器，可以考虑用来近似。
窗函数的基本思想：先选取一个理想滤波器（它的单位抽样响应是非因果、无限长的），再截取（或加窗）它的单位抽样响应得到线性相位因果FIR滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器。
设是一个长序列，是长度为N的窗函数，用截断，得到N点序列，即
（1-3）
在频域上则有
（1-4）
由此可见，窗函数不仅仅会影响原信号在时域上的波形，而且也会影响到频域内的形状。
MATLAB信号工具箱主要提供了以下几种窗函数，如表1-1所示：
表1-1 MATLAB窗函数
窗
窗 函 数
矩形窗
Boxcar
三角窗
Triang
海明窗
Hamming
汉宁窗
Hanning
布莱克曼
Blackman
切比雪夫窗
Chebyshev
凯塞窗
Kaiser
加矩形窗后的频谱和理想频谱可得到以下结论：
加窗使过渡带变宽，过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是。N越大，过渡带越窄、越陡；
过渡带两旁产生肩峰，肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95％，与N无关。工程上习惯用相对衰耗来描述滤波器，相对衰耗定义为：
（1-5）
这样两个肩峰点的相对衰耗分别是0.74dB和-21dB。其中（-0.0895）对应的点的值定义为阻带最小衰耗。
以上的分析可见，滤波器的各种重要指标都是由窗函数决定，因此改进滤波器的关键在于改进窗函数。
窗函数谱的两个最重要的指标是：主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为：
旁瓣峰值衰耗＝20lg(第一旁瓣峰值／主瓣峰值) （1-6）
为了改善滤波器的性能，需使窗函数谱满足：
(1)主瓣尽可能窄，以使设计出来的滤波器有较陡的过渡带。
(2)尽量减少最大旁瓣的相对幅度，也就是能量集中于主瓣，以减小带内、带外波动的最大幅度，增大阻带衰减。
一般来说，以上两点很难同时满足。当选取主瓣宽度很窄时，旁瓣的分量势必增加，从而带内、带外的波动也增加了；当选取最小的旁瓣幅度时，降低了带内、带外的波动，但是过渡带的陡度减小了。所以实际采用的窗函数其特性往往是它们的折中，在保证主瓣宽度达到一定要求的前提下，适当牺牲主瓣宽度来换取旁瓣波动的减小。

2.基本窗函数
基本窗函数包括以下三类，主要有：矩形窗函数、汗明窗窗函数和布莱克曼窗函数。下面就对这些窗函数展开介绍。
2.1 矩形窗函数
矩形窗(Rectangular Window)函数的时域形式可以表示为：
(2-1) （公式2(-(1）
它的频域特性为
(2-2) （公式2-2）
Boxcar函数：生成矩形窗
调用方式w = boxcar (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。从功能上讲，该函数又等价于w = ones(n,1)。

2.2汉宁窗（升余弦窗）
汉宁窗（Hanning）函数时域形式可表示为：
（2-3）
利用傅利叶变换的调制特性，由上式可得汉明窗的平谱函数为：

（2-4） （式2-8）
式中， (公式2-9)
当N远大于1时，上式可近似表示为：
（2-5） （公式2-10）
这三部分之和使旁瓣互相抵消，能量更集中在主瓣，汉明窗函数的最大旁瓣值比主瓣值低31dB，但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍，为8π/N。
hanning函数：生成汉明窗
调用方式：
(1) w = hanning(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。
注意：此函数不返回是零点的窗函数的首尾两个元素。
(2) w = hanming(n,'symmetric')：与上面相类似。
(3) w = hanming(n,'periodic')：此函数返回包括为零点的窗函数的首尾两个元素。
2.3布莱克曼窗函数
布莱克曼窗（Blackman Window)函数的时域形式可以表示为

（2-6）
其频谱函数为

（2-7）

其幅度函数

（2-8）
3方案设计与论证
用窗函数法设计一个FIR带通滤波器。指示如下：
低端阻带截止频率 wls = 0.2*pi；
低端通带截止频率 wlp = 0.35*pi；
高端通带截止频率 whp = 0.65*pi；
高端阻带截止频率 whs = 0.8*pi；
六种窗函数的基本参数如表1-2：
表1-2 窗函数基本参数

窗函数
旁瓣峰值幅度/dB
过渡带宽
阻带最小衰减/dB
矩形窗
-13
4π/N
-21
三角形窗
-25
8π/N
-25
汉宁窗
-31
8π/N
-44
哈明窗
-41
8π/N
-53
布莱克曼窗
-57
12π/N
-74
凯塞窗()
-57
10π/N
-80
以上表格里的参数设置是最佳窗函数设计，根据设计方案的要求，选择一个合适的窗函数进行滤波器的设计，从上表可以看出：最小带阻衰减仅有窗函数决定，不受N的影响，而过渡带的宽度则随窗函数的增加而减小。

3.1 程序设计法
程序法只以hamming窗为示例，其他窗函数的程序大体相同，只是在window= hanning(N)和hn = fir1(N-1,wc, hanning (N))处要调自己的窗函数即可。
hanning 窗设计程序示例：
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= hanning(N);
[h1,w]=freqz(window,1);
figure(1);
stem(window);
axis([0 60 0 1.2]);
grid;
xlabel('n');
title('Hanning窗函数');
figure(2);
plot(w/pi,20*log(abs(h1)/abs(h1(1))));
axis([0 1 -350 0]);
grid;
xlabel('w/pi');
ylabel('幅度(dB)');
title('Hanning窗函数的频谱');
hn = fir1(N-1,wc, hanning (N));
[h2,w]=freqz(hn,1,512);
figure(3);
stem(n,hn);
axis([0 60 -0.25 0.25]);
grid;
xlabel('n');
ylabel('h(n)');
title(‘Hanning窗函数的单位脉冲响应’);
figure(4);
plot(w/pi,20*log(abs(h2)/abs(h2(1))));
grid;
xlabel('w/pi');
ylabel('幅度(dB)');

4 窗函数仿真结果分析
4.1 矩形窗
仿真结果如图4-1所示：

图4-1
4.2 汉宁窗
仿真结果如图4-2所示：

图4-2
4.3布莱克曼窗
仿真结果如图4-3所示：
)

图4-3

5语音信号分析
利用Windows下的录音机，录制一段自己的话音，时间在1s内。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，我们很快理解了采样频率、采样位数等概念。这里我直接采用了一段现成的.wav格式的语音信号。
5.1语音信号的频谱分析
程序如下：
clear;clc;close all;
[S,Fs,Bit] = wavread('F:/SHE.wav'); %读取音频信号
sound(S,Fs);
fnoise1 = 10000; % 加入频率为fnoise的正弦噪声信号
fnoise2 = 50; % 加入频率为fnoise的正弦噪声信号
N = length(S);
T= N/Fs;
t=T/N:T/N:T;
for j=1:N
Noise1(j)=sin(2*pi*fnoise1*t(j));
Noise2(j)=sin(2*pi*fnoise2*t(j));
end
% 对于带通滤波器：
S1 = S+0.1*Noise1'+Noise2';
fs1=1000;
fp1=1200;
fs2=3200;
fp2=3000;
sound(S1,Fs);
dw=2*pi*(fp1-fs1)/Fs;
wls = 2*fs1*pi/Fs;
wlp = 2*fp1*pi/Fs;
whp = 2*fs2*pi/Fs;
wc = [wlp/pi,whp/pi];
B = wlp-wls;
N = ceil(6.2/dw);
n=0:N-1;
h1 = fir1(N-1,wc, hanning (N));
X1 = abs(fft(S));
Z1 = (0:150000)/150000;
figure(1);
plot(Z1,X1(1:150001));
X2 = abs(fft(S1));
Z2 = (0:150000)/150000;
figure(2);
plot(Z2,X2(1:150001));
S2=conv(S1,h1);
sound(S2,Fs);
X3=abs(fft(S2));
Z3 = (0:150000)/150000;
figur(3);
plot(Z3,X3(1:150001));
波形图如图5-1所
图
图51图5-1
5.2对语音信号进行滤波
图形如图5-2所示：

图5-2
5.3语音信号的回放
在Matlab中，函数sound可以对声音进行回放。其调用格式：sound(x，fs，bits)，x为要播放的音频信号，fs采样频率，bits采样位。在运行Matlab程序时，可以听到回放的声音，经过比较，明显感觉滤波前后的声音有变化。原来的语音信号与滤波后的语音信号进行对比，有如下的特点：语音高频成分音质非常尖锐，齿音中，声音有些暗淡；语音低频成分音质沉稳，空间感觉强，语音浑厚；语音中频成分音质有力度，有通透感波形如图5-3所示：
图图1图1图5图图图

图5-3

6.结论与展望
本次MATLAB课程设计主要任务是完成FIR滤波器的设计，对我来说这个题目是很有挑战性的，因为自己对MATLAB中有关滤波器的设计知识了解甚少，为了完成老师布置的任务就不得不要强迫自己动手，整合思路，查找资料，为己所用。平时所学的理论知识只是基础，真正应用软件做设计的时候才能知道自己的局限性。一味停留在老师的教学中自己能做的实在是少之又少。老师只是在较高的层次上为自己的学习指明道路，为数字信号处理的整体概念指出思路。至于具体的某个程序要怎么编写，某个新后要怎么处理，不可能手把手的交给自己。所以就应该学会利用资料，首先就是互联网，然后是图书馆。由于本次课设的时间限制，最合理的资料应该是互联网，快速，方便。搜集到资料以后不能照抄，应该仔细阅读，读懂，然后根据自己的要求改变参数。总之，只有知道怎么自己学习，才能知道怎么自己动手。
在设计过程中，我也遇到了很多不懂得地方，程序经常出现错误，尤其是在利用所设计的滤波器对模拟出来的数字序列滤波时，出现很多错误，但经过自己的不断努力和尝试，最终还是解决了问题。
此次的课程设计收获颇多的，经过此次课程设计不断加深巩固了数字信号处理知识还学会了如何使用MATLAB，同时也加深了同学之间的友谊。更学会了自主学习新知识的一种方法。

参考文献
[1] 葛哲学，精通MATLAB.电子工业出版社,2008
[2] 陈亚勇，MATLAB信号处理详解.人民邮电出版社,2008
[3] 维基百科，http://zh.wikipedia.org
[4] 周开利，邓春辉，MATLAB基础及其应用教程.北京大学出版社,2007
[5] 赵静，张瑾，基于MATLAB的通信系统仿真.北京航空航天大学出版社,2006
[6] 宋寿鹏，数字滤波器设计及工程应用.江苏大学出版社,2007
[7] 高西全，丁玉美．数字信号处理．第3版．北京：西安电子科技大学出版社，2008
[8] 刘泉，阙大顺．数字信号处理原理与实现．北京：电子工业出版社，2005
[9] 张磊，毕靖，郭莲英．MATLAB实用教程．北京：人民邮电出版社，2008

附录
矩形窗
程序如下
clear all;
fs1=1000;
fp1=1200;
fs2=3200;
fp2=3000;
Fs=44100
dw=2*pi*(fp1-fs1)/Fs;
wls = 2*fs1*pi/Fs;
wlp = 2*fp1*pi/Fs;
whp = 2*fs2*pi/Fs;
wc = [wlp/pi,whp/pi];
B = wlp-wls;
N = ceil(1.8/dw);
n=0:N-1;
window= rectwin(N);
[h1,w]=freqz(window,1);
hn = fir1(N-1,wc, rectwin (N));
figure(1)
freqz(window)
figure(2)
freqz(hn)
汉宁窗
程序如下：
clear all;
fs1=1000;
fp1=1200;
fs2=3200;
fp2=3000;
Fs=44100

dw=2*pi*(fp1-fs1)/Fs;
wls = 2*fs1*pi/Fs;
wlp = 2*fp1*pi/Fs;
whp = 2*fs2*pi/Fs;
wc = [wlp/pi,whp/pi];
B = wlp-wls;
N = ceil(6.2/dw);
n=0:N-1;
window= hanning(N);
[h1,w]=freqz(window,1);
hn = fir1(N-1,wc, hanning (N));
figure(1)
freqz(window)
figure(2)
freqz(hn)
布莱克曼窗
程序如下：
clear all;
fs1=1000;
fp1=1200;
fs2=3200;
fp2=3000;
Fs=44100
dw=2*pi*(fp1-fs1)/Fs;
wls = 2*fs1*pi/Fs;
wlp = 2*fp1*pi/Fs;
whp = 2*fs2*pi/Fs;
wc = [wlp/pi,whp/pi];
B = wlp-wls;
N = ceil(11/dw);
n=0:N-1;
window= blackman(N);
[h1,w]=freqz(window,1);
hn = fir1(N-1,wc, blackman(N));
figure(1)
freqz(window)
figure(2)
freqz(hn