****学院
在设计IIR滤波器时,常用的方法是利用模拟滤波器来设计数字滤波器。广泛采取这种方法的因素有:1,模拟滤波器设计技术已非常成熟;2,可得闭合形式的解;3,关于模拟滤波器设计有完整的设计公式和图表可以利用很查阅。
为实现从模拟滤波器到数字滤波器的转换,需要从系统的描述方法来考虑转单位采样响应h(n)换问题,无论是模拟滤波器还是数字滤波器,描述系统的基本方法都有四种。
滤波器描述系统的方法
因此,IIR滤波器的设计方法是首先将数字滤波器的技术指标转化为对应模拟滤波器的技术指标,然后设计满足技术指标的模拟滤波器Ha(s),然后将设计出的模拟滤波器Ha(s)转换为满足技术指标的数字滤波器H(z)。将Ha(s)转换成H(z)的最终目的,是希望数字滤波器的频率响应H(ejω)尽量接近模拟滤波器Ha(jΩ)。将系统函数H(z)从s平面转换到z平面的方法有很多种,但工程上常用的有两种:一种是使数字滤波器的h(n)近似于模拟滤波器的ha(t),可导出脉冲响应不变法;另一种使数字滤波器的差分方程近似于模拟滤波器的微分方程,由此可导出双线性变换法。
3.10双线性变换法的基本原理
脉冲响应不变法使得数字滤波器在时域上能够较好的模仿模拟滤波器,但是由于从s平
sT面到z平面的映射z?e具有多值性,使得设计出来的数字滤波器不可避免的出现频谱混
迭现象。为了克服脉冲响应不变法可能产生的频谱混跌效应的缺点,我们使用一种新的变换——双线性变换。双线性变换法可认为是基于对微分方程的积分,利用对积分的数值逼近的道德。
仿真滤波器的传递函数H(s)为
H(s)??cskk?0Nkk?1Mk,M?N
k?dsH(s)??N将展开为部份分式的形式,并假设无重复几点,则
Ak s?sk?1pk那么,对于上述函数所表达的数字信号处理系统来讲,其仿真输入x(t)和模拟输出y(t)有
如下关系
y?(t)?spy(t)?Ax(t)
利用差分方程来代替导数,即
y?(t)?同时令
y(n)?y(n?1)
Ty(t)?1?y(n)?y(n?1)? 21x(t)??x(n)?x(n?1)?
2第5 页 共 17 页
这样,便可将上面的微分方程写为对应的差分方程形式
****学院
s1A?y(n)?y(n?1)??p?y(n)?y(n?1)???x(n)?x(n?1)? T22两边分别取z变换,可得
H(z)?Y(z)A? X(z)21?z?1??spT1?z?1这样,通过上述过程,就可得到双线性变换中的基本关系,如下所示
21?z?1s?? T1?z?12?s z?T2?sT所谓的双线性变换,仅是指变换公式中s与z的关系无论是分子部份还是分母部份都是线性的。
3.11用双线性变换法设计IIR数字滤波器的步骤
MATLAB中设计IIR数字滤波器的具体步骤如下:
(1) 把给出的数字滤波器的性能指标转换为模拟低通滤波器的性能指标; (2) 根据转换后的性能指标,通过滤波器结束选择函数,来确定滤波器的最小阶数n和固有频率wn;
(3) 由最小阶数n得到低通滤波器原型;
(4) 由固有频率wn把模拟低通滤波器转换为模拟低通、高通、带通或带阻滤波器; (5) 运用双线性变换法把模拟滤波器转换成数字滤波器。
3.12程序源代码和运行结果
3.12.1低通滤波器 clear
wp=100*2*pi; %通带截止频率 ws=150*2*pi; %阻带截止频率 rp=0.5; %通带衰减 rs=30; %阻带衰减 fs=2000; %采样频率
[n,wc]=cheb2ord(wp,ws,rp,rs,'s') %计算阶数,与截止频率 [z,p,k]=cheb2ap(n,rs); %建立切比雪夫2型数字滤波器 %零极点转换到空间状态表达式
[a,b,c,d]=zp2ss(z,p,k); %零极点转换到空间状态表达式 [at1,bt1,ct1,dt1]=lp2lp(a,b,c,d,wc); %低通转换到高通
[at2,bt2,ct2,dt2]=bilinear(at1,bt1,ct1,dt1,fs); %双线性变换 [num,den]=ss2tf(at2,bt2,ct2,dt2) %空间状态表达式转换到传递函数 %绘制幅频、相频图(频率响应特性图)
第6 页 共 17 页
****学院
figure(1);
freqz(num,den,128,fs); grid on;
title('幅频、相频图'); %绘制脉冲响应特性图 figure(2);
impz(num,den,128,fs); grid;
title('脉冲响应特性图 ') %滤波检验 figure(3);
t=0:0.0005:0.2;
x=sin(2*pi*50*t)+sin(2*pi*200*t); y=filter(num,den,x); plot(t,x,':',t,y,'-'); grid;
legend('X Signal','Y Signal'); title('滤波检验') 运行结果: n =6
wc =879.2559 num =
0.0287 -0.1085 0.2038 -0.2447 0.2038 -0.1085 0.0287 den =
1.0000 -4.4499 8.4145 -8.6176 5.0302 -1.5841 0.2103
幅频、相频图50Magnitude (dB)0-50-1000100200300400500600Frequency (Hz)7008009001000100Phase (degrees)0-100-200-300-4000100200300400500600Frequency (Hz)7008009001000
图3-2
第7 页 共 17 页
****学院
脉冲响应特性图 0.120.10.080.06 Amplitude0.040.020-0.02-0.04-0.0600.010.020.030.04nT (seconds)0.050.06
图3-3
滤波检验21.510.50-0.5-1-1.5-2X SignalY Signal00.020.040.060.080.10.120.140.160.180.2
图3-4
3.12.2高通滤波器
wp=100*2*pi; %通带截止频率 ws=150*2*pi; %阻带截止频率 rp=0.5; %通带衰减 rs=30; %阻带衰减 fs=2000; %采样频率
[n,wc]=cheb2ord(wp,ws,rp,rs,'s'); %计算阶数,与截止频率 [z,p,k]=cheb2ap(n,rs); %建立切比雪夫2型数字滤波器 [a,b,c,d]=zp2ss(z,p,k); %零极点转换到空间状态表达式
第8 页 共 17 页
****学院
[at1,bt1,ct1,dt1]=lp2hp(a,b,c,d,wc) %低通转换到高通
[at2,bt2,ct2,dt2]=bilinear(at1,bt1,ct1,dt1,fs) %双线性变换
[num,den]=ss2tf(at2,bt2,ct2,dt2); %空间状态表达式转换到传递函数 %绘制幅频、相频图(频率响应特性图) figure(1);
freqz(num,den,128,fs); grid on;
title('幅频、相频图') %绘制脉冲响应特性图 figure(2);
impz(num,den,128,fs); grid;
title('脉冲响应特性图 ') %滤波检验 figure(3);
t=0:0.0005:0.1;
x=sin(2*pi*50*t)+sin(2*pi*200*t); y=filter(num,den,x); plot(t,x,':',t,y,'-');
legend('X Signal','Y Signal'); title('滤波检验') 运行结果: n =6
wc =879.2559
num =0.5365 -3.0688 7.4574 -9.8501 7.4574 -3.0688 0.5365 den =1.0000 -4.5215 8.8272 -9.4731 5.8745 -1.9914 0.2878
幅频、相频图0Magnitude (dB)-20-40-600100200300400500600Frequency (Hz)7008009001000400Phase (degrees)3002001000-1000100200300400500600Frequency (Hz)7008009001000
图3-5
第9 页 共 17 页