在信号系统领域,基于傅里叶变换的谱分析和滤波器设计是两个最基本的问题。在模拟系统中,常常用电子元件逼近所需的传递函数,而在数字系统中,则是用差分方程法去逼近所需传递函数。在硬件实现上,主要是用FPGA中的若干和积门来完成。
1.因果稳定系统
首先还是来看一下系统的性质。这里我们只讨论线性时不变系统(其概念之前讲过),这种系统的性质可完全由其冲激响应决定。在数字系统中,我们用冲激响应h(n)的z变换H(z)来表示系统函数,z=exp(jw)时,则表示其频率响应。
(z有一个收敛域)
正常情况下,我们首先要求系统是稳定的,这里有个重要的定理是:系统稳定的充分且必要条件是其冲激响应绝对可和,如下
再结合上面H(z)的形式,不难得出:H(z)的收敛域要包含z=1(注意这里的前后因果关系,初看优点别扭)。这样就可以用另一种说法来表述上面的定理:系统稳定的充分必要条件是其传递函数的收敛域包含单位圆z=1。
再来看收敛域和极点的关系。一般系统的极点都多余其零点,系统函数可表示为
即可以拆分成多个单极点函数的和。把它变为级数形式
要是级数收敛,n的范围至关重要。对于右边序列(典型的如因果序列),n=(K,+无穷),则收敛域为z>pk;反之,对于左边序列(典型的如逆因果序列),n=(-无穷,K),则收敛域为z<pk。
再结合之前的定理,可得出结论:对于因果系统,要稳定,所有极点都要在单位圆内;对于逆因果系统,要稳定,所有极点都要在单位圆外;对于双边系统,要稳定,往右信号造成的极点在单位圆内,往左信号造成的极点在单位圆外。
2.最小相位系统和全通系统
系统函数中z=exp(jw)是,得到的就是频率响应,其中相频响应也可得到:
现在不看某一个频率处的特定相位,而去考虑频率w在单位圆上转一圈,相位的变化。另外考虑到,单位圆外的零极点,在频率转一圈后,对相位的贡献为0(这不是说,相位随频率不变,而是可能先变大,后变小,反正最终和最初是一样的);而单位圆内的零极点,对相位的贡献为2pi,可以简单画图看出来。则整个频域内,相位变化可表示如下:
mi,ni分别是单位圆内的零极点的个数。
对于因果稳定系统,极点都在单位原内ni=N,则:
即相位永远是延时的,称为相位滞后系统。特别的,若所有零点都在单位圆内,mo=0,则在整个频域内,相位延时差最小,为0。这样的系统即称为最小相位延时系统。反之,所有零点都在单位圆外,称为最大相位延时系统。(注意,它们一定都是因果稳定系统)
对于逆因果稳定系统,极点都在单位圆外ni=0,则:
称为相位超前系统,也有最大最小之分。
最小相位系统在通信中有重要应用,它有一些主要性质:
- 在所有幅频响应相同的系统中,最小相位系统具有最小的相位延时。注意这里说的,不是之前讲的一个系统内,整个通频域内相位变化最小,而是指两个系统间,在任意频率处,它的相位延时都小于另一个非最小相位系统(这才是最小相位的真正含义吧);
- 在所有幅频响应相同的系统中,最小相位系统在任意频率处的群延时(即相频的导数的负值)最小;
- 在所有幅频响应相同的系统中,最小相位系统的冲激响应最集中于n=0处;
- 在所有幅频响应相同的系统中,只有唯一的最小相位系统;
- 任意一个非最小相位因果稳定系统,都可以看成是有一个最小相位系统和一个全通系统的级联。
证明不那么简单,可去参见ppt。
全通系统,顾名思义就是在整个通频带内,幅度响应为1。其特点就是零极点成复倒数对出现,即有一个零点r*ejy,一定有一个极点1/r*ejy(注意它们是复倒数),在图中就是相对于单位圆成对称镜像。例如
它的主要应用有3个:
- 上面说的,任意一个非最小相位系统,都可以看成是一个最小相位系统和一个全通系统的级联;
- 若设计出的系统是非稳定的,可以级联一个全通系统,消除单位圆外的极点,使得系统稳定,且幅频响应不变;
- 相位均衡器,不改变幅频特征的情况下,调整相频响应,使其在通频带内为线性相位。
3.FIR滤波器设计
FIR滤波器很简单,它实际上是一个全零点模型(MA滑动平均模型),滤波器系数只包含滑动平均的B,而自回归的A=1。滤波过程就是X和B的卷积,再联系滤波器的原理,容易知道,FIR滤波器的系数B就是滤波器的冲激响应h。
因此它的设计就简单了,只要设定滤波器的频率响应H,进行ifft后就得到冲激响应h,把h直接作为滤波器系数B对原信号滤波即可。
如上图,是设定的一个低通滤波器的频响,注意了,数字系统有周期延拓的特点。把它ifft得到的冲激响应h(t),注意了,t从负到正的,这里把它平移到正位置,这也是造成FIR滤波器零相位响应变为线性相位响应的原因。还有一点,原来h(t)不是因果的,平移之后就变成因果的了,这样理解对吗??
另外,所得h(t)的长度较长,可截取一段,降低阶数。
如上图所示可见,截取的越长,所得的频响越接近设定频响,第二组就是截取61点(fftl=512)的频响,几乎和设定频响一样。第一组是截取的21点的频响,也大差不差,不过滤波器要做21阶,用硬件做还是比较浪费资源的。在看其相频响应,如下图
下面来看一个例子,用上面设计的FIR滤波器,得到滤波系数h,直接用matlab函数y=filter(h,1,x)进行滤波,即以h为B,A为1进行滤波,得到的谱如图:
用FIR设计带通、高通等滤波器也十分简单,和上面方法一样,如下图:
综上,FIR滤波器设计及其简单,且具有线性相位,对各个频率分量的时间延时都相同。其缺点是,所需阶数一般较多,实现起来浪费硬件资源。
4.IIR滤波器
IIR滤波器设计起来比较复杂,用的是ARMA模型,即A、B都不为1,即系统的传递函数既有零点又有极点。那么,用它的冲激响应,就无法计算出滤波器系数,这是与FIR滤波器的关键区别。
因此也就不能用FIR那样的方法,直接设定频响,IFFT得到冲激响应,从而得到滤波器系数。那怎么办呢?借用模拟滤波器,因为给定一个频响,有很多模拟滤波器的传递函数模型,如巴特沃斯、切毕雪夫等,不过这些传递函数是拉普拉斯形式的,要把它们变成Z变换形式,从而就可得到滤波器系数了。
比如巴特沃斯滤波器模型,其幅度平方为
其参数主要有两个,wc为3dB频点,N为阶数。这两个参数可以通过滤波器指标计算得到,如在200Hz处下降1dB,在600Hz处下降10dB,则可以有两个方程计算出wc和N。
当然也可以简单地设定,比如有信号,带宽500Hz,采样频率2000Hz,我设计一巴特沃斯滤波器,N=2,wc=300。则可以分别计算得在各频率处,幅度的下降值
100Hz | 200Hz | 300Hz | 400Hz | 500Hz |
0.05dB | 0.26dB | 3dB | 6.2dB | 9.4dB |
下面来设计该滤波器来验证一下,由上面的幅度谱平方式,求出其极点并将其分母分解因式,所有的极点肯定是成共轭对称出现的,要使系统稳定,取左侧的极点。由于形式固定,各阶时的极点已做成表格,可以直接查得,不过是wc归一化的。如2阶时有:
更一般性的,最终得到多个单极点模型的和,并用s/wc代替原s,然后根据H(s)到H(z)的变化公式,有:
再把这个和式通分变为一项,计算还是比较复杂的,最终得到
即A=[1 -0.807 0.264],B=[0 0.423],进行滤波得到的波形及频谱如下:
可以看出,300Hz处确实为0.707(3dB),其它各频率处的衰减也都和上表中的吻合。由于阶数较低,滤波器的性能不是很好。若取4阶(计算就更加复杂了,不过有计算机程序专门来设计),其性能就可逼近前面的21阶FIR滤波器,且大大节约资源。
再看其相频响应,因为滤波器的冲激响应是无限的,无法直接得到其频响。这里只能通过滤波前后信号的相位谱来推断滤波器的相频响应了,如下图
可以看出其相位发生了畸变,即说明IIR滤波器没有线性相位,但一般只要保证,在通频带内,相位畸变在一定范围内即可,这还是能满足的,另外还可级联全通系统去均衡相位。
设计IIR带通、高通等滤波器,其原理是先在模拟域(s域)内变换,然后用上述的方法(冲激响应不变法,或其它方法),变为数字滤波器。关键就在于s域内怎么变换,以带通为例,得到截止频率Wc的低通滤波器的H(s)后,将其中的s进行一个替代
用jw代替s得到频率对应关系,并观察Wlp从-Wc到Wc的变化时,Wbp的值:
可见得到的新滤波器是带通的,以为中心频率,为带宽。变换到数字系统,方法和前面一样,计算很复杂,就不做了。其它的各种滤波器也可仿造此法,只是s域的变换关系不同,如高通
5.matlab中的相关函数
看了一下,最有用的从差分方程到零极点的计算函数,这一步虽然理论上很简单,但因式分解的数学计算过于复杂,matlab有函数可以轻松完成这点。如由差分式得到H(z)如下:
直接利用函数[r,p,c]=residuez(B,A),得到的p为极点,r为留数(即分子的系数,终于知道留数是什么了…),可得到零极点模型如下:
可以简单推算一下,确实是准确的。
另外,若给出的H(z)分子分母已是积式,如下:
用函数ploy(p)可以直接得到系数A,然后再用residuez求出和式的系数。
最后还有一个,有滤波器系数A,B可直接得到冲激响应,impz(B,A),如下:
impz([1 0.3 0.5],[1],8) impz([1 0.3 0.5],[1 -0.1],8)
可看出,A=1时,全零点模型,即FIR,冲激响应和滤波器系数一样;A!=1时,零极点模型,冲激响应是无限长的IIR,但很快趋于收敛。和理论是一致的。
6.响应的各种描述
以一个例题来说明,已知x(n)=0.25nu(n),求解差分方程y(n)-1.5y(n-1)+0.5y(n-2)=x(n),n>=0,初始条件y(-2)=10,y(-1)=4。
首先做z变换,这里就要注意的,要求n>=0的,但还包含有n=-1,-2,所以要分开
这样就完整地包含了初始状态,这本是一个基本知识,不过好像忘了,特此复习一下。
然后做分解
很显然,前一部分为零输入响应,后一部分为零状态响应。
稍作合并,重新分解得
前一部分为通解,即系统形式决定的,后一部分为特解,由输入信号的形式决定的。
再重新合并、分解得
前一部分为暂态响应,很快趋于0,后一部分为稳态响应,它是由在单位圆上的极点造成的,一直持续下去。