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

资源描述

通过状态变量计算IIR 系统的输出教材第四章已经介绍了IIR系统的直接型结构。这里再简单地回顾这种结构，然后阐述编程方法。数字系统的系统函数可以写为（1）式中，子系统（2）由系统函数的零点决定。而子系统（3）则由系统函数的极点决定。式（5）中，写在的左边，表示先实现零点，后实现极点。这样的结构称为IIR系统的直接形式Ⅰ。图1示出这种结构，框中标出用于编程的状态变量，详见下文。图1 IIR系统的直接形式Ⅰ 对于线性时不变系统来说，上述两个子系统的次序可以调换，即（4）这时，先实现系统的极点，后实现零点。对应于式（4）的结构示于图2。合并图中的延迟支路即得图3。这样的结构称为IIR系统的直接形式Ⅱ，又名典范形式、正准型。IIR直接形式通常是指直接Ⅱ型。直接Ⅱ型有以下特点：（1）所需要的延时单元最少。（2）系统调整不方便（调整任何一个系数都会影响全部零、极点）（3）受有限字长影响大。图2 与图1等价的结构图3 IIR系统的直接形式Ⅱ 根据给定的系统结构和输入序列，怎样列出差分方程，求出输出序列？下面给出的程序有助于读者掌握这种方法。这个程序计算图1所示IIR直接形式Ⅰ的输出序列，演示状态变量的用法。通过编程来实现数字网络时，最好在网络中设定状态变量。所谓状态变量就是状态节点变量。一个网络有许多节点，但不都是状态变量节点。状态变量是一组根据已知输入和初始状态，可以计算出输出信号及系统内部任意节点变量值的最少节点变量。图1设定7个状态变量，即。下面给出用C语言编写的函数dir1，。调用这个函数可以实现图1所示IIR直接形式Ⅰ。函数dir1的参数意义如下： Points : 输入参数，输入输出序列点数 Order : 输入参数，系统阶数 H0 : 输入参数，系统函数的常数因子 N[] : 输入参数，系统函数的分子多项式系数组 D[] : 输入参数，系统函数的分母多项式系数组 x[] : 输入参数，输入序列 y[] : 输出参数，输出序列函数dir1的运行过程如下：（1）将状态变量数组和初始化为0。（2）对输入序列每一个样点 a. 求输出y[i] = N[0] * x[i] * H0 b. 累加系统函数的分子多项式对输出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系统的级联实现形式。在这种形式中，系统函数被分解为零、极点形式，即（5）实际上，对IIR 滤波器来说，其系统函数的分子分母多项式是同阶的，即。式（5）中，是零点数组，是极点数组。系统的零、极点可能是复数。但由于和是实数数组，所以复数零、极点必将是共轭成对的。合并共轭成对的零点或极点，可得一个实系数的二次多项式。还可以将实值零、极点合并。合并后，式（5）可以改写为（6a）即IIR系统是由个二阶节级联而成。根据教材第8章的分析，有些带通、带阻IIR滤波器是由个四阶节级联而成的。这样的系统函数可表为（6b）对于任何一种情况，IIR系统的级联结构如图4 所示。其中的基本单元是图5所示的二阶节或四阶节。图4 IIR系统的级联结构图5 IIR系统的级联单元【例1】图6 示出一个切比雪夫Ⅰ型低通IIR滤波器的系统函数的分子、分母二次多项式因子。这个滤波器的阶数为5，第1个二阶节其实是一阶节。把它作为二阶节处理后，这个滤波器就有三个级联的二阶节。图6 一个切比雪夫Ⅰ型低通IIR滤波器的3个二阶节【例2】光盘所附VB可视化程序B49从文件BP401.dat （该文件附于本文末尾）读出一个IIR带通滤波器的数据，用以演示IIR系统级联型结构的实现。这个滤波器有4个级联四阶节（）。图7示出这个滤波器的结构。图中，画出第个四阶节，并标出编程所用的状态变量。给出级联结构的所有参数以及输入序列，求输出序列。解通过状态变量来实现IIR级联系统是很方便的。如图16所示，第节的状态变量和分别代表该节的输入和输出。对于第节，应根据已知的求出，将赋给下一节的输入端的状态变量，然后开始计算下一节。对于系统实时性的要求是：在系统下一个输入样点到来之前，应求出了当前输入样点的输出。程序B49 用VB6,0 写成。这里，用C语言译出，供读者参考。主程序main 调用函数 SOS4以计算级联IIR系统的输出序列。图7 一个带通IIR滤波器的级联结构子程序 SOS4 定义各个参数如下： secNum ：输入参数，级联二阶节数目 points ：输入参数，输入、输出序列点数 H0 ：输入参数，系统函数的常数因子 a[ ] ：输入参数，二维数组，代表二阶节的反馈参数 b[ ] ：输入参数，二维数组，代表二阶节的正向传输参数 inSig[ ] ：输入参数，输入序列 outSig[ ] ：输出参数，输出序列函数SOS4的运行过程如下：（1）将所有四阶节的状态变量数组初始化为0。（2）输入样点序号；四阶节序号。（3）对输入序列每一个样点 a. 求第i节的输出，以输出点状态变量赋给第节的输入点状态变量 b. 四阶节序号刷新，即 c. 是否等于四阶节数目secNum？如果不相等，转去步骤a，计算下一节，如果相等，转去步骤（4）。（2）输入样点序号刷新，即。是否等于要求的计算点数？如果等于，计算完毕，否则，转去步骤（3），处理新样点。程序是否正确，关键在于各个四阶节的状态变量是否被正确地刷新。刷新要点如下： ● 程序运行前，所有状态变量都初始化为0。 ● 某个四阶节开始处理前，应根据流图刷新状态变量。 ● 某个四阶节处理完毕时，应依次刷新状态变量，正确次序是，，，。然后把本节的输出点状态变量赋给下一节的输入点状态变量。 #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 节的输入是系统的输入乘以常数因子 H0 w[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 <- 零次" " *****************" "" "-----" " 以下是便于读出的系统函数数据。各个级联节的分母多项式常数项被归一化为 1。相应地，应修改系统函数的常数因子。" 7 8.48858232100412E-05 0 0 -1 0 1 1.43849941664248 -3.22745743423451 3.72712153916394 -3.22745743423451 1.43849941664248 1.14776825159894 -3.22745743423451 4.30858386925103 -3.22745743423451 1.14776825159894 1.1018412986817 -3.22745743423451 4.4004377750855 -3.22745743423451 1.1018412986817 0 0 .873406389810402 -1.511584845034 1 .798818791539774 -2.70577969437508 4.06998514488549 -3.0304315878117 1 .877774162422211 -2.85344554747275 4.16283038596367 -3.05068489154574 1 .96012064406542 -3.02121835368777 4.28782914195782 -3.08557281557954 1 */

展开阅读全文