通过状态变量计算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
#include 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 a[i][j] = *ptra++;
for(i=0; i] = inSig[i] * H0;
for(j=0;j1] = 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; ii * 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

*/