淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
1.0000 2.7098 4.6379 5.2252 4.3685 2.7207 Columns 7 through 12
1.2885 0.4561 0.1181 0.0211 0.0023 0.0001
0Magnitude (dB)-100-200-300050010001500Frequency (Hz)200025000Phase (degrees)-200-400-600-800-1000050010001500Frequency (Hz)20002500
图7巴特沃思数字低通滤波器幅频-相频特性
3.2 脉冲响应不变法
顾名思义，脉冲响应不变法就是要求数字滤波器的脉冲响应序列h(n)与模拟滤波器的脉冲响应ha(t)的采样值相等，即
h(n)=ha(t)的关系，我们知道
H(z)
z?eSTt?nT=ha(nT) （3-10）
式中，T为采样周期。根据模拟信号的拉普拉斯变换与离散序列的Z变换之间
=
1T?Ha(S?jk?ks) （3-11）
此式表明，ha(t)的拉普拉斯变换在s平面上沿虚轴，按照周期?s=2?/T延拓后，按式z=eST，进行Z变换，就可以将Ha(s)映射为H(z)。事实上，用脉冲响应不变法设计IIR滤波器，只适合于Ha(s)有单阶极点,且分母多项式的阶次高于分子多项式阶次的情况[5]。将Ha(s)用部分分式表示：
Ha(s)=LT[ha(t)]=?i?1NAis?si （3-12）
式中，LT[·]代表拉普拉斯变换，si为的单阶极点。将Ha(s)进行拉普拉斯反变换，即可得到
ha(t)=?AieSitu(t) （3-13）
i?1N12
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
N式中，u(t)是单位阶跃函数。则ha(t)的离散序列h(n)=ha(nT)=?AieSinTu(nT)
i?1对h(n)进行z变换之后，可以得到数字滤波器的系统函数H(z)
H(z)=?h(n)zn?0??n=
?1?ei?1NAiSiTz?1 （3-14）
对比Ha(s)与H(z)，我们会发现：s域中Ha(s)的极点是si，映射到z平面之后，其极点变成了eSiT，而系数没有发生变化，仍为Ai。因此，在设计IIR滤波器时，我们只要找出模拟滤波器系统函数Ha(s)的极点和系数Ai，通过脉冲响应不变法，代入H(z)的表达式中，即可求出H(z)，实现连续系统的离散化[8]。
但是脉冲响应不变法只适合于设计低通和带通滤波器，而不适合于设计高通和带阻滤波器。因为，如果模拟信号ha(t)的频带不是介于??/T之间，则会在
??/T的奇数倍附近产生频率混叠现象，映射到z平面后，则会在???附近产生频率混叠现象。从而使所设计的数字滤波器不同程度的偏离模拟滤波器在???附近的频率特性，严重时使数字滤波器不满足给定的技术指标。为此，希望设计的滤波器是带限滤波器，如果不是带限的，例如，高通滤波器、带阻滤波器，需要在高通滤波器、带阻滤波器之前加保护滤波器，滤出高于折叠频率?/T以上的频带，以免产生频率混叠现象[9]。但这样会增加系统的成本和复杂性。因此，高通与带阻滤波器不适合用这种方法。
将模拟滤波器转化为数字滤波器，牵涉到一个关键的问题，即寻找一种转换关系，将s平面上的Ha(s)转换成z平面上的H(z)。这里Ha(s)是模拟滤波器的传输函数，H(z)是数字滤波器的系统函数[10]。为了确保转换后的H(z)稳定且满足技术要求，转换关系要满足以下要求：
1)将因果稳定的模拟滤波器转换为数字滤波器后，仍然是因果稳定的。我们知道，当模拟滤波器的传输函数Ha(s)的极点全部位于s平面的左平面时，模拟滤波器才是因果稳定的；对于数字滤波器而言，因果稳定的条件是其传输函数H(z)的极点要全部位于单位圆内。因此，转换关系应是s平面的左半平面映射到z平面的单位圆内。
2)数字滤波器的频率响应与模拟滤波器的频率响应相对应，S平面的虚轴映射为Z平面的单位圆，而响应的频率之间是线性变换关系。
在MATLAB中，脉冲响应不变法的调用函数是impinvar，其调用格式为： a. [bz,az]=impinvar(b,a,fs) b. [bz,az]=impinvar(b,a)
13
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
c. [bz,az]=impinvar(b,a,fs,tol)
该函数的功能是将分子向量为b、分母向量为a的模拟滤波器，转换为分子向量为bz、分母向量为az的数字滤波器。fs为采样频率，单位为Hz，默认值为1Hz。tol指误差容限，表示转换后的离散系统函数是否有重复的极点。
例如，用脉冲响应不变法设计一个契比雪夫型数字低通滤波器，指标要求：通带截止频率?p=1000Hz，阻带截止频率?s=1200Hz，采样频率fs=5000Hz，通带衰减系数Rp=0.3dB，阻带衰减系数Rs=40dB。
程序如下：
wp=1000*2*pi;ws=1200*2*pi;fs=2500;Rp=0.3;Rs=40; [N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s'); %估计滤波器最小阶数 [z,p,k]=cheb1ap(N,Rp); %模拟滤波器函数引用 [A,B,C,D]=zp2ss(z,p,k); %返回状态转移矩阵形式 [AT,BT,CT,DT]=lp2lp(A,B,C,D,Wn); %频率转换
[b,a]=ss2tf(AT,BT,CT,DT); %返回传递函数形式 [bz,az]=impinvar(b,a,fs); %调用脉冲相应不变法 [H,W]=freqz(bz,az); %返回频率响应 plot(W*fs/(2*pi),abs(H)); %画图 grid;
xlabel('frequency/Hz'); ylabel('magnitude'); N,Wn
程序在MATLAB环境下的运行及结果如图8所示：
1.41.21magnitude0.80.60.40.200200400600800frequency/Hz100012001400
图8 契比雪夫型数字低通滤波器幅频特性曲线
N =11 Wn =6.2832e+003
14
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
3.3 完全设计函数法
此法是根据设计要求，直接调用函数来设计数字滤波器。所用到的函数有butter、cheby1、cheb2ord、ellipd以及besself等。Butter用来直接设计巴特沃思数字滤波器，cheby1用来直接设计切比雪夫?型滤波器，cheb2ord用来设计切比雪夫
?型滤波器，ellipd用来设计椭圆滤波器，besself用来设计贝塞尔滤波器。
例如，用完全设计函数法设计一个巴特沃思数字低通滤波器，技术指标要求为：wp=1000;ws=1200;Rp=0.3;Rs=40;fs=8000;
程序如下：
wp=1000;ws=1200;Rp=0.3;Rs=40;fs=8000;
[N,Wn]=buttord(wp/(fs/2),ws/(fs/2),Rp,Rs) %估计滤波器最小阶数 [b,a]=butter(N,Wn);
[H,W]=freqz(b,a); %返回频率响应 plot(W*fs/(2*pi),abs(H)); %画图 grid;
xlabel('Frequency/Hz'); ylabel('magnitude');
程序在MATLAB环境下的运行及结果如图9所示：
1.41.21magnitude0.80.60.40.2005001000150020002500Frequency/Hz300035004000
图9巴特沃思数字低通滤波器幅频特性曲线
N =29 Wn =0.2611
15
淮北煤炭师范学院学士学位论文 基于MATLAB的数字滤波器设计
3.4 语音滤波实例
根据语音信号的特点给出有关滤波器的性能指标：1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz， As=100dB，Ap=1dB；2）高通滤波器性能指标，fc=2800 Hz，fp=3000 Hz As=100dB，Ap=1dB；3）带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB。用IIR对其分析，最后，利用MATLAB中的函数freqz画出各滤波器的频率响应。
1．语音信号的采集：
1）利用windows下的录音机（开始—程序—附件—娱乐—录音机，文件—属性—立即转换—8000KHz，8位，单声道），录制一段自己的话音“信号”， 时间控制在1秒左右，然后将音频文件保存“E:\\耿博.wav”
2）在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。
[z1,fs,bits]=wavread('E:\\耿博.wav') 运行结果：fs =44100 bits =16
wavread函数调用格式：
y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。 [y,fs,nbits]=wavread(file)，采样值放在向量y中，fs表示采样频率（Hz），nbits表示采样位数。
y=wavread(file,N)，读取前N点的采样值放在向量y中。
y=wavread(file,[N1,N2])，读取从N1点到N2点的采样值放在向量y中。 2．语音信号的频谱分析 ①首先画出语音信号的时域波形 z1=wavread('E:\\耿博.wav'); plot(z1);图像输出如图10
②对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性
z1=wavread('E:\\耿博.wav'); y1=z1(1:8192); Y1=fft(y1); n=0:8191;
plot(n,Y1);图像输出如图11
16