在文章如何快速设计一个FIR滤波器(一)以及如何快速设计一个FIR滤波器(二)等文章中,我们讨论了如何设计FIR(Finite Impulse Response Filter),FIR有很多优点,比如可以获得线性相位,不存在稳定性问题等,典型的FIR具有如下形式:
当
时幅值响应可以变得更为陡峭和理想,这样做的缺点就是需要采集更多的点进行计算,有时候这么做不太适用。因此工程上很多时候采用的是IIR(Infinite Impulse Response Filter),IIR和FIR最大的区别就是输出不仅仅取决于输入,还取决于输出,IIR的标准形式如下:
IIR的设计和FIR有较大的区别,接下来我们就简单讨论一下,先从一个最常见的传递函数说起。
一、从一个最常见的传递函数说起
惯性环节是我们碰到最多的传递函数了,其基本形式如下:
这个函数有什么特点呢?——无外乎就是从幅值响应和相位响应两个角度来看了。我们只考虑稳态时系统的响应,也就是令
,于是传递函数形式可以改写为:
我们来简单的分析一下这个传递函数,很显然:
当
时,
;
当
时,
;
当
时,
;
也就是说,
在低频是幅值响应大(没衰减,为1),频率增加到
时幅值响应下降到
,当频率增加到
时,幅值响应降到了几乎为零,因此这个滤波器应该是一个低通滤波器。下面我们从数学的角度严谨的看一下,
幅值响应为:
相位响应为:
画个图更直观:
可见,这确是一个低通滤波器。我们知道截止频率的定义是当输出幅值响应下降到输入幅值的-3dB (
),也就是0.707(也就是
)时对应的频率。因此,此处
就是截止频率。
二、如何设计一个高通模拟滤波器
对于低通滤波器,令截止频率
,形式为:
,我们把这个称之为标准化的低通滤波器。那假如我们现在想设计一个高通滤波器
该怎么办呢?
同样令
,我们的目标是:
当
时,
;
当
时,
;
当
时,
;
很简单,将
即可,则传递函数变形为:
我们画个图看一下:
哈哈,果然是高通滤波器,而且截止频率就是
。
三、如何设计一个带通模拟滤波器
那假如我们现在想设计一个带通滤波器该怎么办呢?
也就是我们想要:
当
时,
;
当
时,
;
当
时,
;
当
时,
;
、
分别为带通的高和低端的频率。我们经过几次尝试以后呢发现,可以将标准形式的低通滤波器中
进行如下替换就可以实现:
,其中
,
。
则变换后的传递函数变为:
于是就得到了一个带通滤波器。
四、如何设计一个带阻模拟滤波器
那假如我们现在想设计一个带阻滤波器该怎么办呢?
同样令
,也就是我们想要:
当
时,
;
当
时,
;
当
时,
;
当
时,
;
同样,我们将标准形式的
进行变换,
,则可以得到:
画个图看看:
可见,确实得到了一个带阻滤波器。
好了,需要停下来总结一下,我们有了一个标准的低通模拟滤波器
,通过适当的变形,我们可以得到相应的高通、带通及带阻模拟滤波器。但是我们现在是要设计数字滤波器啊,这不是跑题了吗?——也不见得,按照我们的直觉,模拟滤波器和数字滤波器应该存在某种联系,如果我们找到了这个联系,就可以通过模拟滤波器来得到数字滤波器,问题就解决了。那这个联系是什么呢?——双线性变换。
五、什么是双线性变换
在如何理解离散傅里叶变换及Z变换一文中我们介绍了什么是z变换,以及z和s的关系:
,其中
为采样周期。这就是模拟(s)和数字(z)的联系,理论上可以直接用的,指数函数用起来比较复杂,不方便,我们需要一个更简单的形式。
稍微变一下形:
我们知道,函数
的泰勒展开为:
利用上式分别对分子和分母进行泰勒展开,并只取前两项:
稍微变一下形:
这就是双线性变换了,为啥这个名字呢?——因为它是用两个s域的线性函数的比值来逼近z的。我们稍微做一个推演,看看双线性变换代表什么。首先看一下z域的单位圆:
我们知道:在s域,s=0时对应的是零频率,对应到z域就是z=1;
在s域,s=∞时对应的是无穷大频率f=∞。通过非线性变换,我们知道此时z=-1,对应的频率就是
;
因此,双线性变换本质就是就是将s域(模拟量)无穷大的频率映射到z域的
,因为根据香浓采样定理,数字信号的包含的最大分量的频率就是一半的采样频率,否则就会产生混叠。
可见,对于双线性变换而言,在s域的频率和z域对应的频率不同,发生了一定的弯曲,也就意味着截止频率在s域和在z域是不一样的,现在我们需要找到这种关系。
在z域,假设截止频率为
,则
根据双线性变换计算得到的s为
此时s对应模拟的频率为:
所以
即
也就是说对于双线性变而言,模拟的截止频率和数字的截止频率是不同的,不同的原因是因为双线性变换是近似变换,不是准确换算。
值得一提的是,当
时,即采样频率远大于截止频率,可以得到
也就是模拟的截止频率和数字的截止频率差不多,为简单起见,也可以直接用数字的截止频率代替模拟的截止频率。
六、如何设计数字滤波器
有了前面的铺垫,我们就可以进行数字的IIR滤波器设计了,基本有如下五个步骤:选择一个归一化的模拟滤波器
;
确定数字滤波器的截止频率,也就是响应为-3dB(幅值衰减为
)是对应的频率;
利用公式
计算对应的模拟截止频率;
选择合适的变换,得到
比如对于低通滤波器:
其中
。
在
中,用
进行替换,得到数字化的滤波器
。
举个例子,现在假设要设计一个低通滤波器,截止频率为
,采样频率为
。
step1:选择归一化的滤波器
;
step2:数字截止频率为10Hz;
step3:对应的模拟截止频率为:
;
step4:将
中s进行替换:
,
。得到:
,其中
;
step5: 用
进行替换,得到
于是就得到了我们想要的滤波器。
七、如何利用MATLAB进行IIR设计
前面手动设计IIR滤波器是为了告诉大家基本的设计原理,实际工程中我们一般采用计算机来进行设计,快速、准确、可视化程度高,下面我们就来看看如何利用MATLAB来设计IIR滤波器。
最常用的函数就是butter函数,具体语法为:
[a b]= butter(N,Wn,'low');
a和b就是就是IIR滤波器分子和分母对应的系数:
为截止频率,low代表低通滤波器。
我们来验证一下我们前面设计的IIR滤波器对不对。根据定义
,MATLAB中输入[a b]= butter(1,0.4,'low'),其计算结果为:
a =[0.4208 0.4208],b =[1.0000 -0.1584]
可见,是一致的(因为手算位数少,会有一定的计算误差)。
可能有的童鞋就纳闷了,在设计FIR时,用的函数是fir1、fir2等,一看就是FIR滤波器,为啥到IIR函数的名字就叫butter了,难道第一个设计IIR的人喜欢吃黄油?——哈哈,当然不是,其实butter是butterworth(巴特沃斯)的简写,那butterworth又是什么呢?——看拼写像是一个人名,没错,这就是一个人名。那为啥用这个名字呢?
还记得我们一开始选取的标准低通滤波器的函数吗?
稳态时可用
,则
即
这其实是一类函数的特殊形式:
其中
为滤波器的阶次,
为截止频率,当N=1时时就是我们选择的那个标准的一节惯性环节,感兴趣的童鞋可以试一下这类函数的性质。这类函数最早是由一个叫butterworth的人提出的,因此,基于这类函数而衍生的滤波器就是butterworth滤波器,MATLAB就用butter函数来致敬butterworth。
当然基准函数的获取,除了butterworth的方法,还有其他的方法,他们的特点分别为:Butterworth(巴特沃斯)滤波器:滤波器具有单调下降的幅频特性;
Chebyshev(切比雪夫)滤波器:幅频特性在铜带或阻带内有波动,可提高选择性;大约有3/4通带接近线性相位;
Bessel(贝塞尔)滤波器:通带内有较好的线性相位;
Ellipse(椭圆)滤波器:较好的线性相位;大约有1/2通带接近线性相位。
它们的幅值响应如下图所示。
除了butter函数外,当然我们不需要记住其他的命令,当需要其他类型的IIR时,打开MATLAB,在命令行输入:filterDesigner或者fdatool(老版MATLAB),你就能看到如下界面:
通过勾勾选选就能设计IIR滤波器哦!