IIR数字滤波器的实现（C语言）

经典滤波器和数字滤波器

　　一般滤波器可以分为经典滤波器和数字滤波器。

经典滤波器：假定输入信号中的有用成分和希望去除的成分各自占有不同的频带。如果信号和噪声的频谱相互重迭，经典滤波器无能为力。比如 FIR 和 IIR 滤波器等。　　
现代滤波器：从含有噪声的时间序列中估计出信号的某些特征或信号本身。现代滤波器将信号和噪声都视为随机信号。包括 Wiener Filter、Kalman Filter、线性预测器、自适应滤波器等

Z变换和差分方程

　　在连续系统中采用拉普拉斯变换求解微分方程，并直接定义了传递函数，成为研究系统的基本工具。在采样系统中，连续变量变成了离散量，将Laplace变换用于离散量中，就得到了Z变换。和拉氏变换一样，Z变换可用来求解差分方程，定义Z传递函数成为分析研究采样系统的基本工具。

　　对于一般常用的信号序列，可以直接查表找出其Z变换。相应地，也可由信号序列的Z变换查出原信号序列，从而使求取信号序列的Z变换较为简便易行。

　　Z变换有许多重要的性质和定理：

线性定理

　　设a，a₁，a₂为任意常数，连续时间函数f(t)，f₁(t)，f₂(t)的Z变换分别为F(z)，F₁(z)，F₂(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}$$

　　其中A₁=0.5，A₂=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})}$$

　　其中c_i为零点而d_i为极点。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中较大者的个数）；缺点：系数a_k，b_k对滤波器性能的控制关系不直接，因此调整不方便；具体实现滤波器时a_k，b_k的量化误差将使滤波器的频率响应产生很大的改变，甚至影响系统的稳定性。直接型结构一般用于实现低阶系统。

　　级联型：将系统传递函数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，因此被滤除，同时滤波也产生了一定的滞后和失真。

　　知道了差分方程的形式并通过MATLAB得到滤波器系数后很容易写出相应的代码来实现数字滤波，另外还有一个网站能根据设计指标直接生成C代码：http://www-users.cs.york.ac.uk/~fisher/mkfilter/

　　根据前面的设计指标，在网页上填入相应参数后提交，会得到下面的C语言代码。简单修改后就可以使用：

#define NZEROS 2
#define NPOLES 2
#define GAIN   1.482463775e+01

static float xv[NZEROS+1], yv[NPOLES+1];

static void filterloop()
{ 
　　for (;;)
   {
　　　　 xv[0] = xv[1]; xv[1] = xv[2]; 
        xv[2] = next input value / GAIN;
        yv[0] = yv[1]; yv[1] = yv[2]; 
        yv[2] =   (xv[0] + xv[2]) + 2 * xv[1] + ( -0.4128015981 * yv[0]) + (  1.1429805025 * yv[1]);
        next output value = yv[2];
    }
}