1、系统辨识实验报告学院：信息科学与技术学院专业：自动化日期：2016/4/26目录实验 14一实验内容及要求： .4二实验原理： .4三软件设计思想： .4四程序结构框图： .5五运行示意图： .5实验 28一实验内容及要求： .8二实验原理： .8三软件设计思想： .9四程序设计框图： .9五程序运行流程图： .10实验 312一实验内容及要求： .12二实验原理： .12三程序数据流程图： .12四实验运行结果： .13实验 414一实验内容及要求： .14二实验原理： .14三数据递推关系图： .14四实验运行结果： .15心得体会 .16附录（实验代码） .171. LabWork1 .
2、172. LabWork2 .213. LabWork3 .234. LabWork4 .26实验 1一实验内容及要求：1.编出矩阵 A 与 B 相乘得到的矩阵 R 的运算计算机程序要求：（1）A 和 B 的维数及数值可通过键盘及数据文件输入（2）计算结果 R 可由屏幕及文件输出2.将 1 改写为子程序3查找有关的资料，读懂及调通矩阵求逆程序，并改写为子程序。二实验原理：1. 两个矩阵 A、B 相乘得到 C 矩阵，首先要满足的条件是 A 的列数与 B 行数相等，否则不能相乘。当满足条件后，根据 C（i,k）= 可以(,)*(,)求得 C 矩阵。2. 当求矩阵的逆时，首先要判断其是否为方阵，若是
3、则可以对其进行下一步的操作。本次实验中求逆主要是通过构造一个增广矩阵（FangZ | E）矩阵的初等行变换得到（E | FZNi）的这样的一个矩阵就可以求得矩阵的逆。若矩阵 FangZ 不是满秩矩阵时，FangZ 没有 FZNi 。通过这样的求逆方式，避免了大方阵的求取行列式运算。三软件设计思想：1. 确定该软件的功能主要有：键盘输入两个矩阵然后相乘；文本 data 输入两个矩阵将结果放在文本 result 中；键盘输入一个方阵求得其逆矩阵。其中前两个的矩阵相乘运算部分设置为一个函数 Mul。2. 在 main 函数中提供两个关于矩阵的选择：multiplitation；invertion。其
4、相对应的子函数为 MulOp(a,b,c)，Inv()。3. 在 MulOp(a,b,c)子函数中，有两种输入矩阵的方式：way1，way2。相对应的功能为键盘输入，文本输入。并且两者在处理矩阵时，都调用了 Mul 函数。4. 在 Inv()子函数中，输入和显示原矩阵，和其相应的逆矩阵。调用qiuni（double FangZM,double FZNiM,int n）子函数，可以都得到原矩阵的逆矩阵。但当原矩阵不可逆时，系统输出为”The array is not invertible!“。四程序结构框图：五运行示意图：1. main 函数的主界面：2. MulOp 子函数的界面：main 函
5、数Inv 子函数矩阵求逆way2 文本输入矩阵MulOp 子函数矩阵乘积way1 键盘输入矩阵qiuni 子函数Mul 矩阵乘积子函数退出程序3. Inv 子函数的界面：4. 通过键盘操作计算两个矩阵的乘积：5. 求方阵的逆矩阵：实验 2一实验内容及要求：编写并调试动态模型仿真程序：模型：y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)已知：白噪声v(k)数据文件为 DV，数据长度为 L=500要求：（1） 产生长度为 L 的 M 序列数据文件 DU（2） 产生长度为 L 的模型输出数据文件 DY二实验原理：由于在现实中，白噪声序列很难求，所以寻找到
6、 M 序列在一定程度上可以代替白色噪声序列。由 L=500，所以 n=9。根据 M 序列的特征方程： ()=0+1+22+可知 9 阶移位寄存器的多项式为 ，及可得()=9+4+1c=0,0,0,1,0,0,0,0,19 级线性移位寄存器：图中 Ci 表示反馈的两种可能连接方式，Ci=1 表示连线接通，第 9-i 级加入反馈中；Ci=0 表示连线断开，第 9-i 级未参加反馈。系统产生 M 序列的结构流程图：8 7 1 0C0+ +C1 C8输出三软件设计思想：1. 该软件的主要功能是：产生 M 序列赋给 u(k)保存在 DU.txt 文件中；由 u(k)和v(k)求得 y(k)保存在 DY.
7、txt 文件中。2. 在 main 函数中给出 3 个选择：求 u(k)；求 y(k)；退出程序。其相对应的函数名称为 gener，ouput，exit。3. 在 gener 子函数中产生 M 序列 u(k)保存到 DU.txt 文本文件中。4. 在 output 子函数中，通过对 input 子函数（读入 v(k),u(k)的数据） 、deal 子函数(由公式 y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)求 y(k)的调用来达到生成 y(k)序列并保存到 DY.txt 文本文件中。输出 M 序列M 序列 DU(i)=a(i)寄存器向前移 1 位
8、=1=0( 2)移位寄存器c=0,0,0,1,0,0,0,0,1M 序列的长度为 L=500初始化寄存器a=1,1,0,0,0,1,0,1,0i（E|FZNi） */for(row=0;row=1;row-)for(i=row-1;i=0;i-) san=cirow*1.0/crowrow;for(col=row;col2*n;col+)cicol-=crowcol*san;/* 得到标准的（E|A） */for(row=0;rown;row+)for(col=0;coln;col+)crowcol+n/=crowrow;for(i=0;in;i+)for(j=0;jn;j+)FZNiij=c
9、ij+n;2. LabWork2/* 处理 v(k),u(k)的数据得输出 y(k) */void deal(double DVM,double DUM,double DYM)int k;DY0=0;DY1=0;for(k=2;k500;k+)DYk=1.5*DYk-1-0.7*DYk-2+DUk-1+0.5*DUk-2+DVk;/* 产生 m 序列 u(k)输出到文件 */void gener()int a9=1,1,0,0,0,1,0,1,0,c9=0,0,0,1,0,0,0,0,1,aa9;int DUM,temp,sum,i,j;DU0=a8;FILE *fp;system(“cls“
10、);if(fp=fopen(“DU.txt“,“w“)=NULL)printf(“cannot open the file!n“);exit(0);for(i=1;iM;i+)sum=0;for(j=0;j9;j+)sum+=aj*cj;aa0=sum%2;for(j=1;j9;j+)aaj=aj-1;for(j=0;j9;j+)aj=aaj;DUi=a8;for(i=0;iM;i+)fprintf(fp,“%d “,DUi);fclose(fp);printf(“nnnttPlease refer to the file(DU.txt)!nn“);system(“pause“);3. Lab
11、Work3#define M 500#define N 4void MUL1(double TXM,double XN,double ji1N)int i,j,k;for(i=0;iN;i+)for(k=0;kN;k+)for(j=0;jM-3;j+)ji1ik+=TXij*Xjk;void MUL2(double niN,double TXM,double ji2M)int i,j,k;for(i=0;iN;i+)for(k=0;kM-3;k+)for(j=0;jN;j+)ji2ik+=niij*TXjk;void MUL3(double ji2M,double YM,double xish
12、uN)int i,j;for(i=0;iN;i+)for(j=0;jM-3;j+)xishui+=ji2ij*Yj;void main()system(“cls“);system(“color 75“);int i,j;double DUM,DYM;double XMN=0,temp,YM,TXNM,ji1NN=0,niNN=0,ji2NM=0,xishuN=0;read(DU,DY);/* 得到大矩阵 X( 498 * 4 ) */for(i=0;iM-3;i+)for(j=0;j2;j+)Xij=-DYi+j+1;for(j=0;j2;j+)Xij+2=DUi+j+1;temp=0;tem
13、p=Xi0;Xi0=Xi1;Xi1=temp;temp=0;temp=Xi2;Xi2=Xi3;Xi3=temp;/* 得到 Y（498 * 1）=DY(1) DY(498)的列矩阵 */for(i=0;iM-3;i+)Yi=DYi+3;/* 得到 X（498 * 4）的转置矩阵 TX（4 * 498） */for(i=0;iN;i+)for(j=0;jM-3;j+)TXij=Xji; MUL1(TX,X,ji1); /* 得到 TX*X 的乘积 ji1 */qiuni(ji1,ni,N); /* 得到 ji1 的逆矩阵 ni */MUL2(ni,TX,ji2); /* 得到 ni*TX 的乘积
14、 ji2 */MUL3(ji2,Y,xishu); /* 得到 ji2*Y 的乘积 xishu */printf(“tthe xishu are:nn“);for(i=0;iN;i+)printf(“nnt %.4f n“,xishui);printf(“nn“);4. LabWork4#define M 500#define N 4double sitaNN;void MUL(double aN,double bN,double cN,int n1,int n2,int n4)int i,j,k;for(i=0;in1;i+)for(k=0;kn4;k+)for(j=0;jn2;j+)cik
15、+=aij*bjk;void trans(double dN,double eN,int n5,int n6)int i,j;for(i=0;in6;i+)for(j=0;jn5;j+)eij=dji;void consult()system(“cls“);int i,j,k;double DUM,DYM;double txNN=0,pNN,y,xNN,temp,gama;double ji1NN=0,ji2NN=0;double ji3NN=0,ji4NN=0;double ji5NN=0,ji6NN=0,ji7NN=0;read(DU,DY);/* 对 P（ 4 * 4 ）和 sita 进
16、行初始化 */for(i=0;iN;i+)for(j=0;jN;j+)if(i=j) pij=10e9;else pij=0;for(i=0;iN;i+)sitai0=0.0;/* 大循环开始 */for(i=2;iM;i+)/* 求 y(i+1) */y=DYi;/* 求 x(i+1)的(4*1) 的列阵 */for(k=0;k1;k+)for(j=0;j2;j+)xjk=-DYi+j-2;for(j=0;j2;j+)xj+2k=DUi+j-2;temp=x0k;x0k=x1k;x1k=temp;temp=x2k;x2k=x3k;x3k=temp;/* 求 gama(i+1)的值 gama=
17、1/(1+tx(i+1)*p(i)*x(i+1) */trans(x,tx,4,1); /* x(4 1) 到 tx(1 4) */MUL(tx,p,ji1,1,4,4); /* tx(1 4)*p(4 4)=ji1(1 4) */MUL(ji1,x,ji2,1,4,1); /* ji1(1 4)*x(4 1)=ji2(1 1) */gama=1.0/(1+ji200);/* 求 sita(i+1)的（4*1）的列阵 sita(i+1)=sita(i)+gama(i+1)*p(i)*x(i+1)*(y(i+1)-tx(i+1)*sita(i) */MUL(tx,sita,ji3,1,4,1);
18、 /* tx(1 4)*sita(4 1)=ji3(1 1) */MUL(p,x,ji4,4,4,1); /* p(4 4)*x(4 1)=ji4(4 1) */for(j=0;jN;j+)ji4j0=ji4j0*gama*(y-ji300);for(j=0;jN;j+)sitaj0=sitaj0+ji4j0;/*求 p(i+1)的（4*4 ）的矩阵 p(i+1)=p(i)-gama(i+1)*p(i)*x(i+1)*tx(i+1)*p(i) */MUL(p,x,ji5,4,4,1); /* p(4 4)*x(4 1)=ji5(4 1) */MUL(ji5,tx,ji6,4,1,4); /*
19、ji5(4 1)*tx(1 4)=ji6(4 4) */MUL(ji6,p,ji7,4,4,4); /* ji6(4 4)*p(4 4)=ji7(4 4) */for(j=0;jN;j+)for(k=0;kN;k+)pjk=pjk-gama*ji7jk;/* 使得计算过程当中使用的数组都清零 */for(j=0;jN;j+)for(k=0;kN;k+)ji1jk=0;ji2jk=0;ji3jk=0;ji4jk=0;ji5jk=0;ji6jk=0;ji7jk=0;xjk=0;txjk=0;/* 结果显示 */printf(“tthe xishu are:“);for(i=0;iN;i+)printf(“nnntt %.4f n“,sitai0);printf(“nn“);
