通过状态变量计算IIR 滤波器的输出

通过状态变量计算IIR 系统的输出教材第四章已经介绍了IIR 系统的直接型结构。

这里再简单地回顾这种结构，然后阐述编程方法。

数字系统的系统函数可以写为]11][[)()()(100∑∑=-=--==Mk kk Mk k k z a z b H z X z Y z H )()(210z H z H H = （1）式中，子系统∑=-=Mk kk zb z H 01)( （2）由系统函数的零点决定。

而子系统 ∑=--=Mk kk z a Z H 1211)( （3）则由系统函数的极点决定。

式（5）中，)(1z H 写在)(2Z H 的左边，表示先实现零点，后实现极点。

这样的结构称为IIR 系统的直接形式Ⅰ。

图1示出这种结构），（33==M N ，框中标出用于编程的状态变量，详见下文。

图1 IIR 系统的直接形式Ⅰ），（33==M N对于线性时不变系统来说，上述两个子系统的次序可以调换，即)()()(120z H z H H z H = （4） 这时，先实现系统的极点，后实现零点。

对应于式（4）的结构示于图2。

合并图中的延迟支路即得图3。

这样的结构称为IIR 系统的直接形式Ⅱ，又名典范形式、正准型。

IIR 直接形式通常是指直接Ⅱ型。

直接Ⅱ型有以下特点：（1） 所需要的延时单元最少。

（2） 系统调整不方便（调整任何一个系数都会影响全部零、极点） （3） 受有限字长影响大。

图2 与图1等价的结构图3 IIR 系统的直接形式Ⅱ），（33==M N根据给定的系统结构和输入序列，怎样列出差分方程，求出输出序列？ 下面给出的程序有助于读者掌握这种方法。

这个程序计算图1所示IIR 直接形式Ⅰ的输出序列，演示状态变量的用法。

通过编程来实现数字网络时，最好在网络中设定状态变量。

所谓状态变量就是状态节点变量。

一个网络有许多节点，但不都是状态变量节点。

状态变量是一组根据已知输入和初始状态，可以计算出输出信号及系统内部任意节点变量值的最少节点变量。

图1设定7个状态变量，即)2(~)0(),3(~)0(v v u u 。

下面给出用C 语言编写的函数dir1，。

调用这个函数可以实现图1所示IIR 直接形式Ⅰ。

函数dir1的参数意义如下：Points : 输入参数，输入输出序列点数 Order : 输入参数，系统阶数H0 : 输入参数，系统函数的常数因子N[] : 输入参数，系统函数的分子多项式系数组 D[] : 输入参数，系统函数的分母多项式系数组 x[] : 输入参数，输入序列y[] : 输出参数，输出序列 函数dir1的运行过程如下：（1） 将状态变量数组 ] [u 和 ] [v 初始化为0。

（2）对输入序列每一个样点a. 求输出y[i] = N[0] * x[i] * H0b. 累加系统函数的分子多项式对输出y[i]的贡献。

c. 累加求系统函数的分母多项式对输出y[i]的贡献。

d.刷新状态变量数组（请注意刷新顺序：下标值最大的状态变量先被刷新）// 求直接 1 型 IIR 系统的输出void dir1(int points, int order, double H0, double N[], double D[], double x[], double y[]){int i, j, k;double u[20], v[20];for(i = 0; i<=order; i++) // 状态变量数组初始化 {u[i] = 0; v[i] = 0; }for(i = 0; i < points; i++) // 计算每一个输出点 {y[i] = N[0] * x[i] * H0; for(j = 1; j<= order; j++) y[i] = y[i] + u[j] * N[j]; for(k = 1; k<= order; k++) y[i] = y[i] - v[k] * D[k]; // 状态变量数组刷新 for(j = order; j>=2; j--) {u[j] = u[j - 1];v[j] = v[j - 1];}u[1] = x[i] * H0; v[1] = y[i]; } }下面介绍IIR 系统的级联实现形式。

在这种形式中，系统函数被分解为零、极点形式，即∏∏∑∑=-=-=-=---=-=N k kMk kM k kk Mk kk z p z zH z a zb z H 1111010)1()1(1)( （5）实际上，对IIR 滤波器来说，其系统函数的分子分母多项式是同阶的，即N M =。

式（5）中，}{k z 是零点数组，}{k p 是极点数组。

系统的零、极点可能是复数。

但由于}{k a 和}{k b 是实数数组，所以复数零、极点必将是共轭成对的。

合并共轭成对的零点或极点，可得一个实系数的二次多项式。

还可以将实值零、极点合并。

合并后，式（5）可以改写为∏∏==----=++++=Pk k Pk k k k k z H H z z z z H z H 101221122110)(11)(ααββ （6a ） 即IIR 系统是由P 个二阶节级联而成。

根据教材第8章的分析，有些带通、带阻IIR 滤波器是由Q 个四阶节级联而成的。

这样的系统函数可表为∏∏==--------=++++++++=Qk k Qk kk k k k k k k H H z z z z z z z z H z H 1014433221144332211011)(ααααββββ （6b ） 对于任何一种情况，IIR 系统的级联结构如图4 所示。

其中的基本单元是图5所示的二阶节或四阶节。

图4 IIR 系统的级联结构图5 IIR 系统的级联单元【例1】 图6 示出一个切比雪夫Ⅰ型低通IIR 滤波器的系统函数的分子、分母二次多项式因子。

这个滤波器的阶数为5，第1个二阶节其实是一阶节。

把它作为二阶节处理后，这个滤波器就有三个级联的二阶节。

图6 一个切比雪夫Ⅰ型低通IIR 滤波器的3个二阶节【例2】 光盘所附VB 可视化程序B49从文件BP401.dat （该文件附于本文末尾）读出一个IIR 带通滤波器的数据，用以演示IIR 系统级联型结构的实现。

这个滤波器有4个级联四阶节（4=M ）。

图7示出这个滤波器的结构。

图中，画出第k 个四阶节，并标出编程所用的状态变量。

给出级联结构的所有参数以及输入序列)(n x ，求输出序列)(n y 。

解 通过状态变量来实现IIR 级联系统是很方便的。

如图16所示，第k 节的状态变量0k w 和1k w 分别代表该节的输入)(n x k 和输出)(n y k 。

对于第k 节，应根据已知的0k w 求出1k w ，将1k w 赋给下一节的输入端的状态变量0)1(+k w ，然后开始计算下一节。

对于系统实时性的要求是：在系统下一个输入样点到来之前，应求出了当前输入样点的输出。

程序B49 用VB6,0 写成。

这里，用C 语言译出，供读者参考。

主程序main 调用函数 SOS4以计算级联IIR 系统的输出序列。

图7 一个带通IIR 滤波器的级联结构子程序 SOS4 定义各个参数如下：secNum ： 输入参数，级联二阶节数目points ： 输入参数，输入、输出序列点数 H0 ： 输入参数，系统函数的常数因子a[ ] ： 输入参数，二维数组，代表二阶节的反馈参数b[ ] ： 输入参数，二维数组，代表二阶节的正向传输参数 inSig[ ] ： 输入参数，输入序列 outSig[ ] ： 输出参数，输出序列函数SOS4的运行过程如下：（1） 将所有四阶节的状态变量数组初始化为0。

（2） 输入样点序号1+←i i ；四阶节序号0=j 。

（3） 对输入序列每一个样点a. 求第i 节的输出，以输出点状态变量赋给第1+i 节的输入点状态变量b. 四阶节序号刷新，即 1+←j jc. j 是否等于四阶节数目secNum ？如果不相等，转去步骤a ，计算下一节，如果相等，转去步骤（4）。

（2） 输入样点序号刷新，即 1+←i i 。

i 是否等于要求的计算点数？如果等于，计算完毕，否则，转去步骤（3），处理新样点。

程序是否正确，关键在于各个四阶节的状态变量是否被正确地刷新。

刷新要点如下： ● 程序运行前，所有状态变量都初始化为0。

● 某个四阶节开始处理前，应根据流图刷新状态变量2k w 。

● 某个四阶节处理完毕时，应依次刷新状态变量36~k k w w ，正确次序是 56k k w w ←，45k k w w ←，34k k w w ←，23k k w w ←。

然后把本节的输出点状 态变量赋给下一节的输入点状态变量。

#include <math.h> #include <iostream.h// 子程序 SOS4计算级联IIR 系统的输出序列void SOS4(int points, int secNum, double H0, double* ptra, double* ptrb, double inSig[], double outSig[]) {int i, j; double xx[10];double w[30][10]; // 状态变量数组, 第一维代表节序号，第二维代表节内变量序号 double a[10][5], b[10][5]; // 四阶节系数数组, 第一维代表节序号，第二维代表节内变量序号for(i=0; i<secNum; i++) for(j=0; j<5; j++) a[i][j] = *ptra++;for(i=0; i<secNum; i++)for(j=0; j<5; j++)b[i][j] = *ptrb++;// 初始化所有级联节的状态变量for(i=0;i<=30;i++)for(j=0;j<=6;j++)w[i][j] = 0;// 计算系统的输出序列for(i=0; i<points; i++){// 第0 节的输入是系统的输入乘以常数因子H0w[0][0] = inSig[i] * H0;for(j=0;j<secNum;j++){// 计算第j 节的内部状态变量w[j][2]，其中的状态变量w[j][0] 代表该节的输入w[j][2] = w[j][0] - w[j][3] * a[j][1] - w[j][4] * a[j][2] - w[j][5] * a[j][3] - w[j][6] * a[j][4];// 计算第j 节的内部状态变量w[j][1]，此状态变量代表该节的输出w[j][1] = w[j][2] * b[j][0] + w[j][3] * b[j][1] + w[j][4] * b[j][2] + w[j][5] * b[j][3] + w[j][6] * b[j][4];// 下一节的输入是本节的输出w[j + 1][0] = w[j][1];// 刷新状态变量w[j][6] = w[j][5]; w[j][5] = w[j][4]; w[j][4] = w[j][3]; w[j][3] = w[j][2];}// 求出输出序列的一点（最后一节的输出即是系统的输出）outSig[i] = w[secNum - 1][1];}}// 主函数void main(){int i, j;int points; // 点数int secNum; // 四阶节数目double *ptra, *ptrb; // 四阶节系数指针double twoPi = 2 * 3.14159265, H0, f1 = 50, fsamp = 600;double inSig[1024], outSig[1024]; // 输入输出序列double a[10][5], b[10][5]; // 四阶节系数数组// 本程序有一个二阶节和三个四阶节，现把二阶节也作为四阶节处理，故有四个四阶节secNum = 4;// 四阶节分子多项式系数b[0][4] = 0.0; b[0][3] = 0.0; b[0][2] = -1.0;b[0][1] = 0.0; b[0][0] = 1.0;//--------b[1][4] = 1.4384994; b[1][3] = -3.2274574; b[1][2] = 3.7271215;b[1][1] = -3.2274574; b[1][0] = 1.4384994;//--------b[2][4] = 1.1477683; b[2][3] = -3.2274574; b[2][2] = 4.3085839;b[2][1] = -3.2274574; b[2][0] = 1.1477683;//--------b[3][4] = 1.1018413; b[3][3] = -3.2274574; b[3][2] = 4.4004378;b[3][1] = -3.2274574; b[3][0] = 1.1018413;//---------// 四阶节分母多项式系数a[0][4] = 0.0; a[0][3] = 0.0; a[0][2] = .8734064;a[0][1] = -1.5115848; a[0][0] = 1.0;//---------a[1][4] = .7988188; a[1][3] = -2.7057797; a[1][2] = 4.0699851;a[1][1] = -3.0304316; a[1][0] = 1;//---------a[2][4] = .8777742; a[2][3] = -2.8534455; a[2][2] = 4.1628304;a[2][1] = -3.0506849; a[2][0] = 1;//---------a[3][4] = .9601206; a[3][3] = -3.0212184; a[3][2] = 4.2878291;a[3][1] = -3.0855728; a[3][0] = 1;// 四阶节常数因子H0 = 8.4885823E-05;points = 512; // 序列点数// 设定输入序列for(i=0; i<points; i++)inSig[i] = 0.5 * sin(i * twoPi * f1 / fsamp + twoPi / 4);ptra = &a[0][0]; ptrb = &b[0][0];// 计算输出序列SOS4(points, secNum, H0, ptra, ptrb, inSig, outSig);// 打印部分输出序列for(i=110;i<=128;i++)cout << i << " " << outSig[i] << endl;} // 主函数到此终止// -------------------------------------------------------------------------------------------------------- // 带通滤波器文件BP401.dat/*"Elliptic Digital Filter ""********************************""Band Pass"" Sampling Frequency = 600 Hz"" Points in Frequency Domain = 1024" " Pass Band :"" Fp1 = 45 Hz"" Fp2 = 80 Hz"" Apass = 0.2 dB"" Stop Band :"" Fs1 = 35 Hz"" Fs2 = 95 Hz"" Astop = 70 dB""----------------------------------" "Parameters of the Filter :"" 归一化转移函数 H(s) :"" "" 阶数 N = 7"" "" 常数 HO = 1.17844839134034E-04" " "" *****************"" 节号 k = 0----"" 分子多项式系数 "" 0.000000 0.000000 1.000000 <- 零次" " ........."" 分母多项式系数 "" 0.000000 1.000000 0.067574 <- 零次" " *****************"" 节号 k = 1----"" 分子多项式系数 "" 1.000000 0.000000 0.438499 <- 零次" " ........."" 分母多项式系数 "" 1.000000 0.113194 0.012099 <- 零次" " *****************"" 节号 k = 2----"" 分子多项式系数 "" 1.000000 0.000000 0.147768 <- 零次" " ........."" 分母多项式系数 "" 1.000000 0.066814 0.026474 <- 零次" " *****************"" 节号 k = 3----"" 分子多项式系数 "" 1.000000 0.000000 0.101841 <- 零次" " ........."" 分母多项式系数 "" 1.000000 0.021076 0.035930 <- 零次"" *****************"""" 非归一化转移函数 H(s) :"" "" 阶数 N = 7"" "" 常数 HO = 4.89536968039929E-04"" "" *****************"" 节号 k = 0----"" 分子多项式系数 "" 0.000000 0.000000 1.000000 <- 零次"" ........."" 分母多项式系数 "" 0.000000 1.000000 0.280708 <- 零次"" *****************"" 节号 k = 1----"" 分子多项式系数 "" 1.000000 0.000000 7.566915 <- 零次"" ........."" 分母多项式系数 "" 1.000000 0.470216 0.208777 <- 零次"" *****************"" 节号 k = 2----"" 分子多项式系数 "" 1.000000 0.000000 2.549946 <- 零次"" ........."" 分母多项式系数 "" 1.000000 0.277551 0.456845 <- 零次"" *****************"" 节号 k = 3----"" 分子多项式系数 "" 1.000000 0.000000 1.757413 <- 零次"" ........."" 分母多项式系数 "" 1.000000 0.087553 0.620016 <- 零次"" *****************""""-----"" 以下是便于读出的系统函数数据。