1 引言
在许多数字信号处理系统中，如语音或音频信号处理中，有限冲激响应（FIR）滤波器是最常用的组件之一，它完成信号预调、频带选择和滤波等功能。FIR滤波器虽然在截止频率的边沿陡峭性能上不及无限冲激响应（IIR）滤波器，但是却具有严格的线性相位特性，稳定性好，能设计成多通带（或多阻带）滤波器组，所以能够在数字信号处理领域得到广泛的应用。
Matlab是mathswork公司推出的一套高性能的数值计算和可视化软件，它集数值分析、信号运算、矩阵运算、信号处理和图形显示于一体，具有很强的绘图功能。利用它自带的丰富的函数和信号处理工具箱，可以很方便地设计和仿真FIR滤波器组。
2 用窗函数法设计FIR滤波器
一个离散时间系统H（z）=B（z）/A（z），若分母多项式A（z）的系数a1=…=aN=0,那么该系统即变成FIR系统，即
（1）
显然，系数b0，b1，…，bM即是该系统的单位抽样响应h(0) ，h(1) ，…，h(M) ，且当n>M时，h(n)≡0。所对应的FIR滤波器的转移函数为
（2）
H（z）的频率响应应近似为Hd(ejw)，而且是线性相位的。
FIR滤波器有许多行之有效的设计方法，如窗函数法、频率抽样法及最佳一致逼近法等，其中窗函数法很适合于带通滤波器组的设计。
一个截止频率为ωc的理想低通滤波器Hd(ejw)，其冲激响应hd(n)是无限长和非因果的。对hd(n)加窗（窗函数用w(n)表示，窗宽为有限值M）的结果，便得到一个FIR滤波器h(n)，它是对所要求的理想低通滤波器的逼近。如果窗函数w(n)和FIR滤波器h(n)的频率特性分别用W(ejw)和H(ejw)表示，则H(ejw)将等于W(ejw)和Hd(ejw)的卷积。由于窗函数频谱主瓣不是无限窄，因而FIR滤波器的频率特性在通带和阻带之间形成了一个过渡带。由于窗函数频谱中存在着旁瓣，因而FIR滤波器的频率特性在通带和阻带内产生了一些起伏振荡的波纹。
窗函数法设计的基本思想是对于给定的滤波器技术要求，选择滤波器长度M和具有最窄主瓣宽度和尽可能最小的旁瓣衰减的某个窗函数w(n)。常用的窗函数有矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、哈明（Hamming）窗、布莱克曼（Blackman）窗和凯塞（Kaiser）窗，他们的主要区别在于参数主瓣3dB带宽B、最大边瓣峰值A（dB）和边瓣谱峰渐进衰减速度D（Db/oct）之间的差别。
3 FIR滤波器组的设计与Matlab仿真
3.1 FIR滤波器组的设计
用窗函数法设计FIR滤波器组的基本思路为: (1)根据所要求的分析频率范围和频率分辨率选择一组分析频率ωk及相对应的原型理想低通滤波器的截止频率ωck（k=0，1，…，L-1），L为滤波器组的通道数，得出理想带通滤波器组的复合频率响应。(2)用窗函数法近似实现各个理想带通滤波器。实际的复合频率响应等于理想复合频率响应与设计窗的频率特性的卷积。在用窗函数法设计FIR滤波器时，还需要知道通带内的纹波幅度δ，它可以根据阻带衰减量As的要求来确定，通常的典型值为
As=-20lgδ=40∽60 （3）
如果选择Kaiser窗，参数β的计算公式为
当As >50时，β=0.1102(As-8.7)；
当21<As <50时，β=0.5845(As-21)0.5+0.07886(As-21)。 （4）
滤波器的阶数为
M=(As-7.95)/14.36△F+1 (5)
式中△F是归一化过渡带宽，即△F=△f•T，而
△f=△ω/2π （6）
过渡带宽△ω的选择根据对每个滤波器的频率分辨率要求来决定，必须满足限制条件
△ωk<2ωck (7)
如果各通道采用相同的设计窗，保证中心频率及带宽的选择能覆盖-π到π频率范围（不管滤波器的中心频率和带宽如何分布），那么最终得到的实际复合响应一定是理想的，即具有平坦幅度和线性相位。
在数字语音处理中，语音信号取样频率f s=9.6kHz，频谱分析范围为200Hz到3200Hz，设计一FIR带通滤波器组，使阻带内衰减为60dB。选取滤波器通道数L=15，均匀分布，则理想的原型低通滤波器的截止频率为
f c=(3200-200)/(2x15)=100Hz
中心频率为
f k =(200k+100)Hz, k=1,2,…,15
采用Kaiser窗函数设计，要求δ=10-3，取As=-20lgδ=60，则由式（4）、（5）、（6）和（7）可计算出相应的参数为
β=0.1102(As-8.7)=5.65326
△f<2f c=200Hz
△F<△f•T=200Hz/9600Hz
M=(As-7.95)/14.36△F+1≈175
则带通滤波器组的实际复合冲激响应为
h(n)=w(nd) δ(n-nd)
其中延时nd =（M-1）/2=87。
3.2 Matlab仿真
在Matlab平台上进行FIR滤波器组的仿真主要有两种途径：一是基于数据流的仿真，一是基于时间流的仿真。基于数据流的仿真是用Matlab函数以命令行的形式实现，在仿真过程中不能改变初始数据或仿真参数，若想改变只有终止仿真重新操作。基于时间流的仿真是用动态仿真工具Simulink来实现，它提供了很多模块库，只需将所需模块调出进行连接即可，初始数据或仿真参数在仿真过程中可以改变，以便动态地观察系统随之引起的变化。
两种方法相比较而言，基于数据流的仿真有较高的灵活性，对一些特殊的功能模块，由于有Matlab大量函数的支持，用软件编程实现比较容易。而Simulink提供的模块库虽然也很多，但仍然不能满足通信系统的要求，有一定的限制。因此，我们用函数编程来仿真。对前面提到的数字语音处理中的FIR带通滤波器组用Matlab进行仿真，根据系统所给参数编程，其源程序为如下：
%ideal_lp.m
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);
%freqz_m.m
function [db,mag,pha,grd,w]=freqz_m(b,a)
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);
%firkaiser.m
hh=0;
L=15;
As=60;
beta=0.1102*(As-8.7);
delta_f=(3200-200)/L;
delta_F=delta_f/9600;
M=round((As-7.95)/(14.36*delta_F)+1);
for i=1:L
w1=200*i*2*pi/9600;
w2=200*(i+1)*2*pi/9600;
n=[0:1:(M-1)];
hd=ideal_lp(w2,M)-ideal_lp(w1,M);
w_kaiser=(kaiser(M,beta))';
h=hd.*w_kaiser;
hh=hh+h;
[db,mag,pha,grd,w]=freqz_m(h,[1]);
delta_w=2*pi/1000;
subplot(2,1,1);plot(w/pi,db);title('Magnitude Response in dB');grid
axis([0 1 -70 10]);xlabel('frequency in pi units');ylabel('Decibels')
hold on
end
hold off
[db,mag,pha,grd,w]=freqz_m(hh,[1]);
delta_w=2*pi/1000;
subplot(2,1,2);plot(w/pi,db);title('Magnitude Response in dB');grid
axis([0 1 -100 10]);xlabel('frequency in pi units');ylabel('Decibels')
仿真的结果如图1所示，上半部分为各个通道独立的幅频特性，下半部分为滤波器组复合幅频特性。从图中可以看出，阻带衰减达63dB，满足技术指标的要求。
将源程序中滤波器组的幅频特性纵轴输出参数稍做调整，则可以得到图2所示的幅频特性，从图中可以看出通带纹波幅度δ小于0.02dB（即小于10-3），完全满足技术指标的要求。
4 结束语
MATLAB在数字信号处理中的应用越来越广泛，用MATLAB来设计和仿真FIR滤波器组即简便又易于调整参数，是以往其他设计和仿真方法所无法比拟的，是工程技术人员进行系统设计和仿真的得力助手。
参考文献：
[1][美]维纳·K·恩格尔著。Digital Signal Processing Using MATLAB。西安交通大学出版社，2002年
[2][美]Shoichiro Nakamura著。Numerical Analysis and Graphic Visualization with MATLAB (Second Edition)。电子工业出版社，2002年
[3]姚天任著。数字语音处理。华中科技大学出版社，2002年
扫一扫
111
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
