所以H(s)和H(z)的关系为:
H(z)z?esT1?T?H(s?jk?) (2-22)
sk????式(2-22)是将模拟滤波器的系统函数H(s)作周期延拓,在经过式(2-21)的映射变换,映射到Z平面上,从而得到数字滤波器的系统函数H(z)。数字和模拟频率满足:???T的关系。经过z?esT的映射,s平面的左半平面映射为Z平面的单位圆内,因此,一个模拟滤波器的稳定性和因果性经过映射可以得到数字滤波器的稳定性和因果性。
经上面的分析,通过模拟滤波器的系统传递函数G(s)可求得数字滤波器系统传递函数H(z),它的设计具体步骤可以归纳如下:
a.利用???T,使数字滤波的指标?p和?s转换为模拟滤波器指标
?p和?s。
b.根据指标?p和?s来设计模拟滤波H(s)。
c.通过使用部分分式展开法,把H(s)展开为H(s)??d.最后把模拟的极点Pk转换为数字的极点eNAk(2-23) s?Pk?1kNSkT,得到
Ak H(z)??SkT?1 (2-24) Zk?11?e根据上面的理论,在MATLAB环境下用函数实现脉冲响应不变法设计数字低通滤波器。其函数为[bz,az]=impinvar(b,a,T),其中,b表示数字滤波器自变量为Z?1的分子多项式,a表示数字滤波器自变的量Z?1的分母多项式, T表示采样变换参量。
低通,采样频率为 1Hz,通带临界频率 fp =0.2Hz,通带内衰减
27
小于1dB(αp=1);阻带临界频率 fs=0.3Hz,阻带内衰减大于 25dB(αs=25)。设计一个数字滤波器满足以上参数。
FS=1
[n,Wn]=buttord(0.2*2*pi,0.3*2*pi,1,25,'s'); %临界频率采用角频率表示
[b,a]=butter(n,Wn,'s'); %freqs(b,a) %设计模拟的
[bz,az]=impinvar(b,a,FS); %映射为数字的 freqz(bz,az,512,FS) 仿真的结果图:
图 2-12 脉冲响应不变法频率响应(基于MATLAB实现) 对图2-12 分析,脉冲响应不变法的优点:①是频率坐标实线性的即???T,②数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应,时域逼近好。缺点:若抽样频率不高或其他原因将产生混叠失真,不能重现原来的模拟滤波器频率响应。故,脉冲响应不变法不适合带通、低通、高通、带阻等滤波器的设计。
28
2.双线性变换法
采用双线性不变法可以克服脉冲响应不变法产生频率响应的混叠失真,是因为双线性变换法的s平面与z平面是一一映射的关系,消除了多值变换性如下图2-13所示:
图 2-13 双线性变换法的映射关系
为了将S平面的j?轴压缩到s1平面j?1轴上的??/T到?/T一段上,可以通过正切变换实现:
??T???c.tan?1? (2-25)
?2?C为任意常数。令z得z??es1T21?z?1,c为2/T,得s?同样对z求解?1T1?z1?(T2)s这样的变换称为双线性变换,在MATLAB环境下用函
1?(T2)s数对双线性变换的实现如下。
①双线性变换实现 ButterWorth 低通 编程代码: %采样频率 10Hz,通带截止频率 fp=3Hz,阻带截止频率 fs=4Hz
%通带衰减小于 1dB,阻带衰减大于 20dB。
%使用双线性变换法由模拟滤波器原型设计数字滤波器 T=0.1; FS=1/T;
29
fp=3;fs=4; wp=fp/FS*2*pi; ws=fs/FS*2*pi; Rp = 1; %通带衰减 As = 20; %阻带衰减 % 频率预畸
OmegaP = (2/T)*tan(wp/2); % Prewarp Prototype Passband freq OmegaS = (2/T)*tan(ws/2); % Prewarp Prototype Stopband freq %设计 butterworth 低通滤波器原型
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS)));
OmegaC = OmegaP/((10^(Rp/10)-1)^(1/(2*N))); [z,p,k] = buttap(N); %获取零极点参数 p = p*OmegaC; k = k*OmegaC^N; B = real(poly(z)); b0 = k; cs = k*B; ds = real(poly(p)); % 双线性变换
[b,a] = bilinear(cs,ds,FS); % 绘制结果
30
freqz(b,a,512,FS); 仿真出来的结果图如2-13
图 2-13双线性变换实现 ButterWorth 低通 ②双线性变换法实现 Chebyshev 低通(I 型)
采样频率为10Hz,设计一个数字低通滤波器,要求其通带临界频率 fp= 3Hz ,通带αp= 1dB ,阻带临界频率fs= 4Hz ,阻带内衰减大于 15dB
程序代码: T=0.1; FS=1/T; fp=3;fs=4; Rp = 1; As = 15; % 频率预畸 wp=fp/FS*2*pi; ws=fs/FS*2*pi; OmegaP = (2/T)*tan(wp/2);
31