经典滤波器和数字滤波器
⼀般滤波器可以分为经典滤波器和数字滤波器。
经典滤波器:假定输⼊信号中的有⽤成分和希望去除的成分各⾃占有不同的频带。如果信号和噪声的频谱相互重迭,经典滤波器⽆能为⼒。⽐如 FIR 和 IIR 滤波器等。
现代滤波器:从含有噪声的时间序列中估计出信号的某些特征或信号本⾝。现代滤波器将信号和噪声都视为随机信号。包括 Wiener Filter、Kalman Filter、线性预测器、⾃适应滤波器等
Z变换和差分⽅程
在连续系统中采⽤拉普拉斯变换求解微分⽅程,并直接定义了传递函数,成为研究系统的基本⼯具。在采样系统中,连续变量变成了离散量,将Laplace变换⽤于离散量中,就得到了Z变换。和拉⽒变换⼀样,Z变换可⽤来求解差分⽅程,定义Z传递函数成为分析研究采样系统的基本⼯具。 对于⼀般常⽤的信号序列,可以直接查表出其Z变换。相应地,也可由信号序列的Z变换查出原信号序列,从⽽使求取信号序列的Z变换较为简便易⾏。
Z变换有许多重要的性质和定理:
线性定理
设a,a1,a2为任意常数,连续时间函数f(t),f1(t),f2(t)的Z变换分别为F(z),F1(z),F2(z),则有$$\mathbf{Z}[af(t)]=aF(z)$$ $$ \mathbf{Z}[a_1f_1(t)+a_2f_2(t)]=a_1F_1(z)+a_2F_2(z)$$
滞后定理
设连续时间函数在t<0时,f(t)=0,且f(t)的Z变换为F(z),则有$$\mathbf{Z}[f(t-kT)]=z^{-k}F(z)$$
应⽤Z变换求解差分⽅程的⼀个例⼦:已知系统的差分⽅程表达式为$y(n)-0.9y(n-1)=0.05u(n)$,若边
界条件$y(-1)=1$,求系统的完全响应。
解:⽅程两端取z变换$$Y(z)-0.9[z^{-1}Y(z)+y(-1)]=0.05\frac{z}{z-1}\\Y(z)=\frac{0.05z^2}{(z-1)(z-0.9)}+\frac{0.9y(-1)z}{z-
0.9}$$
可得$$\frac{Y(z)}{z}=\frac{A_1z}{z-1}+\frac{A_2z}{z-0.9}$$
其中A1=0.5,A2=0.45,于是$y(n)=0.5+0.45 \times(0.9)^n \quad(n\geq0)$
IIR数字滤波器的差分⽅程和系统函数
IIR数字滤波器是⼀类递归型的线性时不变因果系统,其差分⽅程可以写为:$$y(n)=\sum_{i=0}^{M}a_ix(n-i)+\sum_{i=1}^{N}b_iy(n-i)$$
进⾏Z变换,可得:$$Y(z)=\sum_{i=0}^{M}a_iz^{-i}X(z)+\sum_{i=1}^{N}b_iz^{-i}Y(z)$$
于是得到IIR数字滤波器的系统函数:$$H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{i=0}^{M}a_iz^{-i}}{1-\sum_{i=1}^{N}b_iz^{-
i}}=a_0\frac{\prod_{i=1}^{M}(1-c_iz^{-1})}{\prod_{i=1}^{N}(1-d_iz^{-1})}$$
其中ci为零点⽽di为极点。H(z)的设计就是要确定系数、或者零极点,以使滤波器满⾜给定的性能指标。
IIR数字滤波器结构
数字滤波器的功能本质上是将⼀组输⼊数字序列通过⼀定的运算后转变为另⼀组输出数字序列。滤波器系统函数可以表达为多种不同的形式,每⼀种对应着不同的算法,也就对应着不同的实现结构。例如:$$H(z)=\frac{1}{1-0.3z^{-1}-0.4z^{-2}}$$
可以分解为:$$H(z)=\frac{1}{1-0.8z^{-1}}\cdot \frac{1}{1+0.5z^{-1}}$$
或$$H(z)=\frac{0.6154}{1-0.8z^{-1}}+\frac{0.3846}{1+0.5z^{-1}}$$
上述同⼀系统的三种不同描述形式就对应了不同的实现结构,或者说不同的滤波器结构可以实现相同的传递函数。IIR滤波器常见的结构形式有直接Ⅰ型、直接Ⅱ型(典范型)、级联型、并联型。通过差分⽅程能够画出包含反馈结构的数字⽹络称为直接型。
直接Ⅰ型滤波器的⽹络结构可以根据差分⽅程很直观地画出(The Direct-Form I structure implements the feed-forward terms first followed by the feedback terms):
可以看出直接Ⅰ型滤波器需要N+M个延时单元(N≥M)。直接Ⅱ型结构是对直接Ⅰ型的变型,将上⾯⽹络的零点与极点的级联次序互换,并将延时单元合并。实现N阶滤波器只需要N个延时单元(The Direct-Form II structure uses fewer delay blocks than Direct-Form I),故称为典范型。
直接Ⅱ型看上去不那么直观,可以通过下图进⾏理解。我们可以将整个滤波器系统看成A、B两个⼦系统串联⽽成,由于为线性系统因此交换顺序不影响最终输出结果,传递函数可写为:$$H(z)=B(z)\cdot \frac{1}{A(z)}=\frac{1}{A(z)}\cdot B(z)$$
直接型优点:直接型结构简单,⽤的延迟器较少(N和M中较⼤者的个数);缺点:系数ak,bk对滤波器性能的控制关系不直接,因此调整不⽅便;具体实现滤波器时ak,bk的量化误差将使滤波器的频率响应产⽣很⼤的改变,甚⾄影响系统的稳定性。直接型结构⼀般⽤于实现低阶系统。
级联型:将系统传递函数H(z)因式分解为多个⼆阶⼦系统,系统函数就可以表⽰为这些⼆阶⼦系统传递函数的乘积。实现时将每个⼆阶⼦系统⽤直接型实现,整个系统函数⽤⼆阶环节的级联实现。
并联型:与级联型类似,⽤部分分式展开法将系统函数表⽰为⼆阶⼦系统传递函数的和。每个⼆阶⼦系统仍然⽤直接型实现,整个系统函数⽤⼆阶环节的并联实现。
在IIR滤波器设计过程中,通常利⽤模拟滤波器来设计数字滤波器,要先根据滤波器的性能指标设计出相应的模拟滤波器的系统函数H(s),然后由H(s)经变换得到所需要的数字滤波器的系统函数H(z)。常⽤的变换⽅法有冲激响应不变法和双线性变换法。下⾯使⽤MATLAB等⼯具直接⽣成数字滤波器系数:
在MATLAB命令⾏中输⼊fdatool打开滤波器设计⼯具箱,为了便于分析,我们先从设计⼀个简单的2阶低通滤波器。Design Method⽤于选择IIR滤波器还是FIR滤波器,这⾥我们选择IIR滤波器,类型选择Butterworth,当然也可以选择其他类型,不同类型的频率响应不同,选择后默认的滤波器结构是直接II型。ResponseType⽤于选择低通、⾼通、带通、带阻等类型,选择低通滤波“Lowpass”。
Frequency Specifications⽤于设置采样频率以及截⽌频率,这⾥填⼊200以及20,即采样率为200Hz,20Hz以上的频率将被滤除掉。Fiter Order 选择滤波器阶数,为了简单起见,先选择⼀个2阶滤波器做实验。
参数设置好后点击Design filter按钮,将按要求设计滤波器。默认⽣成的IIR滤波器类型是Direct-Form II,Second-Order Sections(直接Ⅱ型,每个Section是⼀个⼆阶滤波器),在⼯具栏上点击Filter Coefficients图标或菜单栏上选择Analysis→Filter Coefficients可以查看⽣成的滤波器系数。
MATLAB中⼆阶滤波器差分⽅程公式如下(注意反馈项符号为负号):$$y[n]=b_0 \cdot x[n]+b_1 \cdot x[n-1] + b_2 \cdot x[n-2] -a_1 \cdot y[n-1] -a_2 \cdot y[n-2]$$
⾼阶IIR滤波器的实现是采⽤⼆阶滤波器级联的⽅式来实现的。默认情况下,Filter Coefficients把结果分成多个2阶Section显⽰,其中还有增益。增益的⽬的是为了保证计算的精度和系统的稳定性。选择[edit]→[Convert to Single Section],这时候系数变成我们熟悉的形式:
按照上⾯的公式,滤波器差分⽅程为:y[n] = 0.06745527*x[n] + 0.134910547*x[n-1] + 0.06745527*x[n-2] - (-
1.1429805025)*y[n-1] - (0.412801596)*y[n-2]
滤波器设计完成后还可以⽣成Simulink模型进⾏仿真:按照下图中数字标号进⾏,第⼀步点击左边Realize Model图标,第⼆步勾
选“Build model using basic elements”这⼀项,右边四个灰⾊的项将⾃动打钩,最后点击“Realize Model”,matlab将⾃动⽣成滤波器模型,在弹出的窗⼝中双击模型可以观察该模型的内部结构。
下⾯是按照设计要求⽣成的2阶滤波器直接Ⅰ型的结构:
Direct-Form I
下⾯是直接Ⅱ型的内部结构:
Direct-Form II
使⽤⽣成的滤波器搭建⼀个简答的测试模型:将两个幅值为1,频率分别为10Hz、50Hz的正弦波叠加,输⼊滤波器后观察滤波前后的波形。仿真时间设为1s,仿真参数中求解器类型设为固定步长,求
解器选择discrete(它适⽤于离散⽆连续状态的系统),步长设为
0.005s(200Hz)
点击Run按钮开始进⾏仿真:
打开⽰波器结果如下图所⽰:上⾯⼀栏是不同频率叠加的波形,下⾯是10Hz正弦波和滤波后得到的波形的对⽐。由于50Hz正弦波频率⾼于滤波器截⽌频率20Hz,因此被滤除,同时滤波也产⽣了⼀定的滞后和失真。
根据前⾯的设计指标,在⽹页上填⼊相应参数后提交,会得到下⾯的C语⾔代码。简单修改后就可以使⽤:
#define NZEROS 2