淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。
用窗函数wd(n)将hd(n)截断，并进行加权处理，得到
根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。
验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止[9]。
h(n)?hd(n)?(n)如果要求线性相位特性， 则h(n)还必须满足：
h(n)??h(N?1?n)
4.1.2典型的窗函数
1、矩形窗(Rectangle Window)
w(n)?RN(n)
其频率响应和幅度响应分别为：
N?1sin(N?/2)?j?2sin(N?/2)j?W(e)?e，WR(?)?
sin(?/2)sin(?/2)2、三角形窗(Bartlett Window)
N?1?2n,0?n??2 w(n)??N?12nN?1?2?,?n?N?1N?12?其频率响应为：
N?12sin(N?/4)2?j?2j?W(e)?[]e
Nsin(?/2)3、汉宁(Hanning)窗，又称升余弦窗
12n?w(n)?[1?cos()]RN(n)
2N?1其频率响应和幅度响应分别为：
22
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
?j(2?2?W(e)?{0.5WR(?)?0.25[WR(??)?WR(??)]}eN?1N?1?W(?)e?j?aj?N?1)?2
W(?)?0.5WR(?)?0.25[WR(??2?2?)?WR(??)]N?1N?14、汉明(Hamming)窗，又称改进的升余弦窗
2n?w(n)?[0.54?0.46cos()]RN(n)
N?12?2?其幅度响应为：W(?)?0.54WR(?)?0.23[WR(??)?WR(??)]
N?1N?15、布莱克曼(Blankman)窗，又称二阶升余弦窗
2n?4n?w(n)?[0.42?0.5cos()?0.08cos()]RN(n)
N?1N?1其幅度响应为：
W(?)?0.42WR(?)?0.25[WR(??2?2?)?WR(??)]N?1N?1 4?4??0.04[WR(??)?WR(??)]N?1N?16、凯泽(Kaiser)窗
I0(?1?[1?2n/(N?1)]2)w(n)?,0?n?N?1
I0(?)其中：β是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，β越大，过渡带越宽，阻带越小衰减也越大。I0(·)是第一类修正零阶贝塞尔函数。
若阻带最小衰减表示为As??20log10?s，β的确定可采用下述经验公式：
?0????0.5842(As?21)0.4?0.07886(As?21)?0.1102(A?8.7)s?As?2121?As?50 As?50若滤波器通带和阻带波纹相等即δp=δs时，滤波器节数可通过下式确定：
N?As?7.95?1
14.36?F其中：
???s??p?F??
2?2?在MATLAB中，实现矩形窗的函数为boxcar和rectwin，其调用格式如下： w=boxcar(N)
23
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
w=rectwin(N)
其中N是窗函数的长度，返回值w是一个N阶的向量，它的元素由窗函数的值组成。实际上，w=boxcar(N)等价于w=ones(N,1)。
在MATLAB中，实现三角窗的函数为triang，调用格式为：
w=triang(N)
在MATLAB中，实现汉宁窗的函数为hann，调用格式如下：
w=hann(N) w=hann(N, 'sflag')
Hann函数中的参数sflag为采样方式，其值可取symmetric(默认值)或periodic。当sflag＝symmetric时，为对称采样；当sflag＝periodic时，为周期采样，此时hann函数计算N+1个点的窗，但是仅返回前N个点。
在MATLAB中，实现海明窗的函数为hamming，调用格式分别如下：
w=hamming (N) w=hamming (N,'sflag')
其中sflag的用法同上。
在MATLAB中，实现布拉克曼窗的函数为blackman，调用格式如下：
w=blackman (N) w=blackman (N,'sflag')
在MATLAB中，实现切比雪夫窗的函数为chebwin，调用格式为：
w=chebwin (N,r)
其中r 表示切比雪夫窗函数的傅里叶变换旁瓣幅度比主瓣低rdB(其默认值为100dB)，且旁瓣是等纹波的。
在MATLAB中，实现巴特里特窗的函数为bartlett，调用格式为：
w=bartlett (N)
在MATLAB中，实现凯塞窗的函数为kaiser，调用格式为：
w=kaiser (N,beta)
其中beta为窗函数的参数β。
在MATLAB中，提供了基于窗函数法的两类设计函数，即函数fir1和函数fir2。 1、 函数fir1
该函数实现加窗的线性相位FIR数字滤波器，可设计标准低通、带通、高通和带阻滤波器[5]。
其调用格式如下： a. b=fir1(n,Wn) b. b=fir1(n,Wn, 'ftype') c. b=fir1(n,Wn,window) d. b=fir1(n,Wn, 'ftype', window) e. b=fir1(?,'normalization')
n表示滤波器的阶数。'ftype'代表所设计滤波器的类型：'high'表示高通滤波器；'stop'表示带阻滤波器；?DC?-1表示多通带滤波器，第一频带为通带；?DC?-0表示多通带滤波器，第一频带为阻带；默认时代表低通或带通滤波器。?window?为窗函
24
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
数，是长度为N的列向量，默认时函数自动取Hamming窗。b=fir1(n,Wn)可得到n阶低通FIR滤波器，调用后返回维数为n+1的行向量b，它是滤波器的系数。b与FIR滤波器的系统函数有如下关系：
H(z)?b(1)?b(2)z?1???b(n?1)z?n
对于高通和带阻滤波器,n取偶数，?n为滤波器的截止频率，范围是(0,1)；对于带通和带阻滤波器，?n=[?1,?2]，且?1??2；对于多通带滤波器，
?n=[?1,?2,?3,?4]，频段为0????1,?1????2,?2????3,?。
2、函数fir2
该函数用于设计基于窗函数的任意响应FIR滤波器，其频率响应由向量f和向量m共同决定，取值在[0，1]之间；n为滤波器阶数；b向量为返回滤波器的系数；window为窗函数，长度为n+1，默认时为hamming窗；npt为对频率响应进行内插点数，默认时为512；lap参数用于指定fir2在重复频率点附近插入的区域大小。 4.2 FIR滤波器滤波实例 1、窗函数设计低通 程序设计如下： clear;close all
[z1,fs,bits]=wavread('E:\\耿博.wav') y1=z1(1:8192); Y1=fft(y1);
fp=1000;fc=1200;As=100;Ap=1;Fs=8000; wc=2*pi*fc/Fs; wp=2*pi*fp/Fs; wdel=wc-wp; beta=0.112*(As-8.7); N=ceil((As-8)/2.285/wdel); wn= kaiser(N+1,beta); ws=(wp+wc)/2/pi; b=fir1(N,ws,wn); figure(1); freqz(b,1); x=fftfilt(b,z1); X=fft(x,8192);
25
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
figure(2);
subplot(2,2,1);plot(abs(Y1));axis([0,1000,0,1.0]); title('滤波前信号频谱');
subplot(2,2,2);plot(abs(X));axis([0,1000,0,1.0]); title('滤波后信号频谱'); subplot(2,2,3);plot(z1); title('滤波前信号波形'); subplot(2,2,4);plot(x); title('滤波后信号波形'); sound(x,fs,bits);
图形分析如图18、图19：
滤波前信号频谱100Magnitude (dB)滤波后信号频谱10.80.60.40.210.800.60.40.2-100-20000.10.20.30.40.50.60.70.8Normalized Frequency (?? rad/sample)0.9100500滤波前信号波形100000500滤波后信号波形10000Phase (degrees)10.50-0.500.10.20.30.40.50.60.70.8Normalized Frequency (?? rad/sample)0.9110.50-0.5024x 1064-2000-4000-6000-8000-1-1024x 1064图18 FIR数字低通滤波器幅频-相位特性曲线 图 19 滤波前后信号频谱和波形对比
2、窗函数设计高通 程序设计如下： clear;close all
[z1,fs,bits]=wavread('E:\\耿博.wav') y1=z1(1:8192); Y1=fft(y1);
fp=2800;fc=3000;As=100;Ap=1;Fs=8000; wc=2*pi*fc/Fs; wp=2*pi*fp/Fs; wdel=wc-wp; beta=0.112*(As-8.7); N=ceil((As-8)/2.285/wdel); wn= kaiser(N,beta); ws=(wp+wc)/2/pi;
26