1、精選優(yōu)質(zhì)文檔-傾情為你奉上精選優(yōu)質(zhì)文檔-傾情為你奉上專(zhuān)心-專(zhuān)注-專(zhuān)業(yè)專(zhuān)心-專(zhuān)注-專(zhuān)業(yè)精選優(yōu)質(zhì)文檔-傾情為你奉上專(zhuān)心-專(zhuān)注-專(zhuān)業(yè)實(shí)驗(yàn)一 矩陣相乘與求逆實(shí)驗(yàn)要求請(qǐng)編出矩陣A和B相乘得到矩陣R的運(yùn)算計(jì)算機(jī)程序，要求A和B的維數(shù)及數(shù)值可以通過(guò)鍵盤(pán)及數(shù)據(jù)文件輸入，計(jì)算結(jié)果R可由屏幕及文件輸出。將題一的程序改為子程序。查找有關(guān)的資料，讀懂及調(diào)通矩陣求逆子程序，并改寫(xiě)為子程序。實(shí)驗(yàn)過(guò)程鍵盤(pán)輸入矩陣A和B，屏幕上顯示計(jì)算結(jié)果矩陣R。#include#include int main() int i,j,k,m,n,p,q; printf(請(qǐng)輸入A矩陣的行數(shù)：); scanf(%d,&m); printf(請(qǐng)輸

2、入A矩陣的列數(shù)：); scanf(%d,&n); printf(請(qǐng)輸入B矩陣的行數(shù)：); scanf(%d,&p); printf(請(qǐng)輸入B矩陣的列數(shù)：); scanf(%d,&q); if(n!=p) printf(請(qǐng)檢查A矩陣的列數(shù)與B矩陣的行數(shù)是否一致！); else double *matrixA,*matrixB,*matrixC; matrixA=(double *)malloc(sizeof(double *)*m); matrixA0=(double *)malloc(sizeof(double)*m*n); matrixB=(double *)malloc(sizeof(do

3、uble *)*p); matrixB0=(double *)malloc(sizeof(double)*p*q); matrixC=(double *)malloc(sizeof(double *)*m); matrixC0=(double *)malloc(sizeof(double)*m*q); for(i=1;im;i+) matrixAi=matrixAi-1+n; for(i=1;ip;i+) matrixBi=matrixBi-1+q; for(i=1;im;i+) matrixCi=matrixCi-1+q; printf(請(qǐng)輸入A矩陣的數(shù)：n); for(i=0;im;i+)

4、 for(j=0;jn;j+) scanf(%lf,&matrixAij); printf(n輸入的A矩陣是：n); for(i=0;im;i+) for(j=0;jn;j+) printf(%7.3lft,matrixAij); printf(n); printf(請(qǐng)輸入B矩陣的數(shù)：n); for(i=0;ip;i+) for(j=0;jq;j+) scanf(%lf,&matrixBij); printf(n輸入的B矩陣是：n); for(i=0;ip;i+) for(j=0;jq;j+) printf(%7.3lft,matrixBij); printf(n); for(i=0;im;i

5、+) for(j=0;jq;j+) matrixCij=0; for(i=0;im;i+) for(j=0;jq;j+) for(k=0;kn;k+) matrixCij=matrixCij+matrixAik*matrixBkj; printf(n矩陣A*B的結(jié)果是：n); for(i=0;im;i+) for(j=0;jq;j+) printf(%7.3lft,matrixCij); printf(n); free (matrixA); free (matrixB); free (matrixC); getchar(); getchar();運(yùn)行的結(jié)果如下：實(shí)現(xiàn)文本的矩陣輸入與計(jì)算結(jié)果的輸

6、出，并將相乘改為子程序。#includestdio.hint main() void matrix_mutipl(double *matrixA,double *matrixB,double *matrixC,int m,int n,int q); char ch;int m=0,n=0,p=0,q=0,i=0,j=0,k=0;char infileA80,infileB80,outfile80;FILE *fpin1,*fpin2,*fpin3;double *matrixA,*matrixB,*matrixC;printf(Enter the infileA name:n); scanf(

7、%s,infileA); if(fpin1=fopen(infileA,r)=NULL) printf(Cant open file %sn,infileA); exit(0); do ch=fgetc(fpin1); if(ch=t) n+; if(ch=n) m+; while(!feof(fpin1); n=(n+m)/m; printf(輸入的矩陣A行與列數(shù)分別為：%d,%dn,m,n); fclose(fpin1); fflush(stdin); matrixA=(double *)malloc(sizeof(double *)*m); matrixA0=(double *)mall

8、oc(sizeof(double)*m*n); for(i=1;im;i+) matrixAi=matrixAi-1+n; for(i=0;im;i+) for(j=0;jn;j+) matrixAij=0; if(fpin1=fopen(infileA,r)=NULL) printf(Cant open file %sn,infileA); exit(0); do ch=fgetc(fpin1); while(!feof(fpin1)&(ch9); if (feof(fpin1) exit(0); ungetc(ch,fpin1); do i=k/n; j=k%n; fscanf(fpin1

9、,%lf,&matrixAij); k+; while(!feof(fpin1)&km*n); fclose(fpin1); printf(輸入的A矩陣為：n); for(i=0;im;i+) for(j=0;jn;j+) printf(%lf,matrixAij); printf(t); printf(n); fflush(stdin); printf(Enter the infileB name:n); scanf(%s,infileB); if(fpin2=fopen(infileB,r)=NULL) printf(Cant open file %sn,infileB); exit(0)

10、; do ch=fgetc(fpin2); if(ch=t) q+; if(ch=n) p+; while(!feof(fpin2); q=(p+q)/p; printf(輸入的矩陣B行與列數(shù)分別為：%d,%dn,p,q); fclose(fpin2); fflush(stdin); matrixB=(double *)malloc(sizeof(double *)*p); matrixB0=(double *)malloc(sizeof(double)*p*q); for(i=1;ip;i+) matrixBi=matrixBi-1+q; for(i=0;ip;i+) for(j=0;jq;

11、j+) matrixBij=0; if(fpin2=fopen(infileB,r)=NULL) printf(Cant open file %sn,infileB); exit(0); do ch=fgetc(fpin2); while(!feof(fpin2)&(ch9); if (feof(fpin2) exit(0); ungetc(ch,fpin2); k=0; do i=k/q; j=k%q; fscanf(fpin2,%lf,&matrixBij); k+; while(!feof(fpin2)&kp*q); fclose(fpin2); printf(輸入的B矩陣為：n); f

12、or(i=0;ip;i+) for(j=0;jq;j+) printf(%lf,matrixBij); printf(t); printf(n); fflush(stdin); matrixC=(double *)malloc(sizeof(double *)*m); matrixC0=(double *)malloc(sizeof(double)*m*q); for(i=1;im;i+) matrixCi=matrixCi-1+q; for(i=0;im;i+) for(j=0;jq;j+) matrixCij=0; matrix_mutipl(double*)matrixA,(double

13、*)matrixB,(double*)matrixC,m,n,q); printf(矩陣A*B為：n); for(i=0;im;i+) for(j=0;jq;j+) printf(%lf,*(*(matrixC+i)+j); printf(t); printf(n); fflush(stdin); printf(Enter the outfile name:n); scanf(%s,outfile); if(fpin3=fopen(outfile,w)=NULL) printf(Cant open file %sn,outfile); exit(0); for(i=0;im;i+) for(j

14、=0;jq;j+) fprintf(fpin3,%lf,*(*(matrixC+i)+j); fprintf(fpin3,t); fprintf(fpin3,n); fflush(stdin); getchar(); void matrix_mutipl(double *matrixA,double *matrixB,double *matrixC,int m,int n,int q) int k=0,j=0,i=0; for(i=0;im;i+) for(j=0;jq;j+) for(k=0;kn;k+) *(*(matrixC+i)+j)=*(*(matrixC+i)+j)+*(*(mat

15、rixA+i)+k)*(*(*(matrixB+k)+j);運(yùn)行結(jié)果如下： 文本文件1.dat和2.dat需要放在程序的同一文件夾下，否則需要說(shuō)明路徑。矩陣求逆子程序#includestdio.h#includemalloc.h#includemath.h void main() int inv(double *p,int n); double a44=1,2,0,0,1,2,0,0,0,0,3,0,0,0,0,1,*ab; ab=a0; int n=4,i=0,j; i=inv(ab,n); if(i!=0) for(i=0;in;i+) putchar(n); for(j=0;jn;j+)

16、 printf(%f ,aij); getchar();int inv(double *p,int n) void swap(double *a,double *b); int *is,*js,i,j,k,l; for(i=0;in;i+) putchar(n); for(j=0;jn;j+) printf(%f ,*(p+i*n+j); puts(nnnn); double temp,fmax; is=(int *)malloc(n*sizeof(int); js=(int *)malloc(n*sizeof(int); for(k=0;kn;k+) fmax=0.0; for(i=k;in

17、;i+) for(j=k;jfmax) fmax=temp; isk=i;jsk=j; if(fmax+1.0)=1.0) free(is);free(js); printf(no inv); return(0); if(i=isk)!=k) for(j=0;jn;j+) swap(p+k*n+j,p+i*n+j); if(j=jsk)!=k) for(i=0;in;i+) swap(p+i*n+k,p+i*n+j); pk*n+k=1.0/pk*n+k; for(j=0;jn;j+) if(j!=k) pk*n+j*=pk*n+k; for(i=0;in;i+) if(i!=k) for(j

18、=0;jn;j+) if(j!=k) pi*n+j=pi*n+j-pi*n+k*pk*n+j; for(i=0;i=0;k-) if(j=jsk)!=k) for(i=0;in;i+) swap(p+j*n+i),(p+k*n+i); if(i=isk)!=k) for(j=0;jn;j+) swap(p+j*n+i),(p+j*n+k); free(is); free(js); return 1;void swap(double *a,double *b) double c; c=*a; *a=*b; *b=c;如上所示，本子程序先將輸入的矩陣進(jìn)行全對(duì)角占優(yōu)的方法，避免了矩陣處理過(guò)程中的“病

19、態(tài)”情況。實(shí)驗(yàn)總結(jié)在整個(gè)實(shí)驗(yàn)過(guò)程中，出現(xiàn)了三個(gè)問(wèn)題：文本輸入是實(shí)現(xiàn)矩陣相乘的重要點(diǎn)，由于本人在實(shí)驗(yàn)的文本文件中加了一些提示語(yǔ)言，因此需要判斷何時(shí)是數(shù)據(jù)，再進(jìn)行矩陣數(shù)據(jù)的讀入，因此增加了一個(gè)數(shù)據(jù)判斷過(guò)程。矩陣程序改為子程序時(shí)，需要實(shí)現(xiàn)地址之間的傳遞，我們不能簡(jiǎn)單地進(jìn)行matrixAij處理，應(yīng)該改為*(*(matrixA+i)+j)進(jìn)行處理，而且還要實(shí)現(xiàn)二維數(shù)組的動(dòng)態(tài)分配存儲(chǔ)空間，如下：double *matrixA;matrixA=(double *)malloc(sizeof(double *)*m); matrixA0=(double *)malloc(sizeof(double)*m*n

20、); for(i=1;im;i+) matrixAi=matrixAi-1+n; for(i=0;im;i+) for(j=0;jn;j+) matrixAij=0;最后的動(dòng)態(tài)數(shù)組賦初值也是很重要的，否則在數(shù)據(jù)處理過(guò)程中會(huì)出現(xiàn)錯(cuò)誤。矩陣求逆的子程序中，需要先對(duì)矩陣進(jìn)行全對(duì)角占優(yōu)的處理，能很好避免數(shù)據(jù)處理過(guò)程中一些“病態(tài)”的處理，而且能避免復(fù)雜的計(jì)算量問(wèn)題。實(shí)驗(yàn)二 M序列的生成與輸出計(jì)算實(shí)驗(yàn)要求編寫(xiě)并調(diào)試動(dòng)態(tài)模型仿真程序模型：已知白噪聲v（k）數(shù)據(jù)文件為DV，數(shù)據(jù)長(zhǎng)度L=500要求：（1）產(chǎn)生長(zhǎng)度為L(zhǎng)的M序列數(shù)據(jù)文件產(chǎn)生長(zhǎng)度為L(zhǎng)的模型輸出數(shù)據(jù)文件DY實(shí)驗(yàn)過(guò)程產(chǎn)生M序列 M序列是一種幅值為二值的偽

21、隨機(jī)序列，周期=（2n-1）（其中n為系統(tǒng)的階次，即辨識(shí)的參數(shù)的個(gè)數(shù)）。 有兩種產(chǎn)生方法： 1.采用隨機(jī)函數(shù)。即先隨機(jī)一個(gè)周期的M序列的值（-a或a），再據(jù)此更新其他的數(shù)據(jù)即可。#include#include#include#includemain()int i,k,n=3,m=2,p;int *a,*mfuc;FILE *fp;p=pow(m,n)-1;a=(int *)malloc(sizeof(int *)*p);mfuc=(int *)malloc(sizeof(int *)*500);srand( (unsigned)time( NULL ) );for(i=0;ip;i+) ai

22、=0;for(i=0;ip;i+) ai=0.5+(float)rand()/(RAND_MAX+1); ai=2*ai-1; printf(%dn,ai);printf(nnnnn);for(i=0;i500;i+) k=i%p; mfuci=ak; printf(%dn,mfuci);if(fp=fopen(m.dat,w)=NULL)printf(Cannot open this file!n);exit(0);for(i=0;i500;i+)fprintf(fp,%d,mfuci);fputc(n,fp);fclose(fp);fflush(stdin);getchar();M序列用模

23、二加和門(mén)電路生成。我們以4階的進(jìn)行說(shuō)明，首先給4個(gè)初始狀態(tài)，如0001，再下一位的數(shù)據(jù)為該數(shù)據(jù)的前面第3、4位數(shù)據(jù)的異或值，一次循環(huán)下去，就可以得到一個(gè)周期為15的M序列，因?yàn)榫敌枰獮?，我們這里命0為-1，實(shí)驗(yàn)程序如下：#includemain()int i=0;int a=1,b=1,c=-1,d=-1,e=0;int m=500;int *m1;m1=(int *)malloc(sizeof(int *)*m);FILE *fpin;printf(%dn,a);printf(%dn,b);printf(%dn,c);printf(%dn,d);m10=a;m11=b;m12=c;m13

24、=d; for(i=4;im;i+) e=a; a=b; b=c; c=d; if (e=a) d=-1; else d=1; printf(%dn,d); m1i=d; if(fpin=fopen(u.dat,w)=NULL) printf(Cannot open this file!n); exit(0); for(i=0;im;i+) fprintf(fpin,%d,m1i); fputc(n,fpin); fclose(fpin); free(m1); fflush(stdin); getchar(); 在上面程序運(yùn)行過(guò)后，會(huì)生成一個(gè)u.dat的文件，即為以后的輸入的數(shù)據(jù)。y值的輸出。

25、 輸入的值選用前面我們生成的M序列，這樣可以開(kāi)始計(jì)算對(duì)應(yīng)的輸出值，這里還有白噪聲的干擾，具體的編程如下：#includemain()int i;int num=500;FILE *fpu,*fpy,*fpdv;double u500,dv500,y500;double u_k1=0,u_k2=0,y_k=0,y_k1=0,y_k2=0; for(i=0;inum;i+) ui=0;dvi=0; if(fpu=fopen(u.dat,r)=NULL) printf(Cannot open u.dat file!n); exit(0); i=0; do fscanf(fpu,%lf,&ui); i

26、+; while(!feof(fpu)&i500); fclose(fpu); if(fpdv=fopen(DV,r)=NULL) printf(Cannot open DV file!n); exit(0); i=0; do fscanf(fpdv,%lf,&dvi); i+; while(!feof(fpdv)&i500); fclose(fpdv); y0=0;y1= 0; for(i=2;inum;i+) u_k1=ui-1; u_k2=ui-2; y_k=1.5*y_k1-0.7*y_k2+u_k1+0.5*u_k2+dvi-2; yi=y_k; y_k2=y_k1; y_k1=yi

27、; for(i=0;inum;i+) printf(%lfn,yi); if(fpy=fopen(y.dat,w)=NULL) printf(Cannot open y.dat file!n);exit(0); for(i=0;inum;i+) fprintf(fpy,%lf,yi); fputc(n,fpy); fclose(fpy); fflush(stdin);getchar(); 通過(guò)上面的程序，我們可以得到一個(gè)y.dat的數(shù)據(jù)輸出數(shù)據(jù)，方便我們以后調(diào)用。三、實(shí)驗(yàn)分析采用M序列來(lái)替代白噪聲，是因?yàn)樗鼪](méi)有白噪聲的缺點(diǎn)，如：要求無(wú)限長(zhǎng)時(shí)間等，但它是接近白噪聲的，是比較理想的辨識(shí)輸入的激勵(lì)信

28、號(hào)。數(shù)據(jù)的輸出值，我們是根據(jù)離散函數(shù)進(jìn)行計(jì)算的，是一種理論的數(shù)據(jù)輸出，我們?cè)诖思由习自肼暎阅７聦?shí)際的數(shù)據(jù)輸出。實(shí)驗(yàn)三 成批最小二乘法（BLS）實(shí)驗(yàn)要求 編寫(xiě)并調(diào)試動(dòng)態(tài)離散時(shí)間模型LS成批算法程序。 要求： 1.原始數(shù)據(jù)由DU和DY讀出； 2.調(diào)用求逆及相乘子程序； 3.顯示參數(shù)辨識(shí)結(jié)果。二、實(shí)驗(yàn)過(guò)程成批最小二乘法的計(jì)算流程系統(tǒng)模型：或?qū)懗桑浩渲?，u（k）我們選為M序列產(chǎn)生的數(shù)據(jù)，即u.dat里的數(shù)據(jù)，v（k）為白色隨機(jī)噪聲。即可寫(xiě)成以下形式：式中，對(duì)于k=n+1,n+2,.,n+m，以上方程可以構(gòu)成一個(gè)線性方程組?？梢詫?xiě)成：其中，參數(shù)估計(jì)其中u(k)應(yīng)?。?）隨機(jī)序列（白噪聲）；2）偽隨即序

29、列（M序列）。具體的程序如下：#includemain() void matrix_mutipl(float *matrixA,float *matrixB,float *matrixC,int m,int n,int q); float inv(float *p,int n); void swap(double *a,double *b); int i,j; float u500; float *y,*x,*xt,*c,*d,*o,*yt;y=(float *)malloc(sizeof(float *)*500);y0=(float *)malloc(sizeof(float)*500*1

30、);for(i=1;i500;i+) yi=yi-1+1; x=(float *)malloc(sizeof(float *)*498);x0=(float *)malloc(sizeof(float)*498*4);for(i=1;i498;i+) xi=xi-1+4; xt=(float *)malloc(sizeof(float *)*4);xt0=(float *)malloc(sizeof(float)*4*498);for(i=1;i4;i+) xti=xti-1+498; c=(float *)malloc(sizeof(float *)*4);c0=(float *)mallo

31、c(sizeof(float)*4*4);for(i=1;i4;i+) ci=ci-1+4; d=(float *)malloc(sizeof(float *)*4);d0=(float *)malloc(sizeof(float)*4*498);for(i=1;i4;i+) di=di-1+498; o=(float *)malloc(sizeof(float *)*4);o0=(float *)malloc(sizeof(float)*4*1);for(i=1;i4;i+) oi=oi-1+1; yt=(float *)malloc(sizeof(float *)*498);yt0=(fl

32、oat *)malloc(sizeof(float)*498*1);for(i=1;i498;i+) yti=yti-1+1; for(i=0;i500;i+) ui=0;for(i=0;i500;i+)yi1=0; for(i=0;i498;i+)for(j=0;j4;j+)xij=0;for(i=0;i4;i+)for(j=0;j498;j+)xtij=0;for(i=0;i4;i+)for(j=0;j4;j+)cij=0; for(i=0;i4;i+)for(j=0;j498;j+)dij=0; for(i=0;i4;i+)for(j=0;j1;j+)oij=0; for(i=0;i49

33、8;i+)for(j=0;j1;j+)ytij=0; FILE *fpdu,*fpdy;if(fpdu=fopen(u.dat,r)=NULL) printf(Cant open file u.datn); exit(0); i=0; do fscanf(fpdu,%f,&ui); i+; while(!feof(fpdu)&i500); fclose(fpdu); if(fpdy=fopen(y.dat,r)=NULL) printf(Cant open file y.datn); exit(0); i=0; do fscanf(fpdy,%g,&yi0); i+; while(!feof(

34、fpdy)&i500); fclose(fpdy);for(i=0;i498;i+)for(j=0;j1;j+)ytij=yi+2j; for(j=0;j2;j+) for(i=0;i498;i+)xij=(-1)*yi-j+10;for(j=2;j4;j+) for(i=0;i498;i+)xij=ui-j+3;for(i=0;i498;i+)for(j=0;j4;j+) xtji=xij; matrix_mutipl(float *)xt,(float *)x,(float *)c,4,498,4);i=inv(*c,4);if(i!=0) for(i=0;i4;i+) putchar(n

35、); for(j=0;j4;j+) printf(%f ,cij); matrix_mutipl(float *)c,(float *)xt,(float *)d,4,4,498);matrix_mutipl(float *)d,(float *)yt,(float *)o,4,498,1);printf(nn);for(i=0;i4;i+)printf(%gn,oi0);getchar();void matrix_mutipl(float *matrixA,float *matrixB,float *matrixC,int m,int n,int q) int k=0,j=0,i=0; fo

36、r(i=0;im;i+) for(j=0;jq;j+) for(k=0;kn;k+) *(*(matrixC+i)+j)=*(*(matrixC+i)+j)+*(*(matrixA+i)+k)*(*(*(matrixB+k)+j);float inv(float *p,int n) void swap(double *a,double *b); int *is,*js,i,j,k,l; for(i=0;in;i+) putchar(n); for(j=0;jn;j+) printf(%f ,*(p+i*n+j); puts(nnnn); double temp,fmax; is=(int *)

37、malloc(n*sizeof(int); js=(int *)malloc(n*sizeof(int); for(k=0;kn;k+) fmax=0.0; for(i=k;in;i+) for(j=k;jfmax) fmax=temp; isk=i;jsk=j; if(fmax+1.0)=1.0) free(is);free(js); printf(no inv); return(0); if(i=isk)!=k) for(j=0;jn;j+) swap(p+k*n+j,p+i*n+j);/交換指針 if(j=jsk)!=k) for(i=0;in;i+) swap(p+i*n+k,p+i*

38、n+j); /交換指針 pk*n+k=1.0/pk*n+k; for(j=0;jn;j+) if(j!=k) pk*n+j*=pk*n+k; for(i=0;in;i+) if(i!=k) for(j=0;jn;j+) if(j!=k) pi*n+j=pi*n+j-pi*n+k*pk*n+j; for(i=0;i=0;k-) if(j=jsk)!=k) for(i=0;in;i+) swap(p+j*n+i),(p+k*n+i); if(i=isk)!=k) for(j=0;jn;j+) swap(p+j*n+i),(p+j*n+k); free(is); free(js); return 1

39、;void swap(double *a,double *b) double c; c=*a; *a=*b; *b=c;運(yùn)行的實(shí)驗(yàn)結(jié)果如下：三、實(shí)驗(yàn)分析在運(yùn)算過(guò)程會(huì)出現(xiàn)求逆的過(guò)程，因此會(huì)不可避免地增加整個(gè)運(yùn)算量，尤其是當(dāng)整個(gè)系統(tǒng)的階次較高時(shí)候，會(huì)占用較大內(nèi)存。2）由于運(yùn)算的結(jié)果相乘等處理過(guò)程中，會(huì)不可避免地出現(xiàn)偏差，如上，我們得到的最大的數(shù)據(jù)誤差為4.8%，已經(jīng)能較好地辨識(shí)系統(tǒng)參數(shù)。3）成批最小二乘法不適合系統(tǒng)的在線辨識(shí)。實(shí)驗(yàn)四 遞推最小二乘法（RLS）實(shí)驗(yàn)要求編寫(xiě)并調(diào)試動(dòng)態(tài)離散時(shí)間模型LS遞推算法程序。要求：原始數(shù)據(jù)由DU和DY讀出；顯示辨識(shí)結(jié)果；設(shè)置選擇變量決定是否輸出中間結(jié)果。二、實(shí)驗(yàn)過(guò)

40、程 遞推最小二乘法的原理 由于成批最小二乘法存在大矩陣的求逆運(yùn)算，而且占用內(nèi)存較多，不適合在線辨識(shí)。我們出現(xiàn)了遞推最小二乘法，它的基本思想是新的估計(jì)值（k+1）=老的估計(jì)值（k)+修正項(xiàng)模型：其中的x（k)為2n*1維的，為2n*1維y（k+1)=測(cè)量值 標(biāo)量P（k)=信息陣 2n*2n對(duì)稱(chēng)方陣（k）=修正系數(shù) 標(biāo)量遞推的框圖如下：初值選取P（0）=aI，a為足夠大的正數(shù)（1051012）（0）=，為足夠小的數(shù)（0）輸入信號(hào)u(k)一般取M序列實(shí)驗(yàn)編程如下：#include#includemain()void matrix_mutipl(float *matrixA,float *matrix

41、B,float *matrixC,int m,int n,int q);int i,j,k,a=0;FILE *fpdu,*fpdy;float u500;float r=0,deta=0;float *y,*x,*o,*p,*xt,*m1,*m2,*m3,*m4,*m5;y=(float *)malloc(sizeof(float *)*500);y0=(float *)malloc(sizeof(float)*500*1);for(i=1;i500;i+) yi=yi-1+1; for(i=0;i500;i+)yi0=0; for(i=0;i500;i+) ui=0;x=(float *)

42、malloc(sizeof(float *)*4);x0=(float *)malloc(sizeof(float)*4*1);for(i=1;i4;i+) xi=xi-1+1; for(i=0;i4;i+) xi0=0;o=(float *)malloc(sizeof(float *)*4);o0=(float *)malloc(sizeof(float)*4*1);for(i=1;i4;i+) oi=oi-1+1; for(i=0;i4;i+) oi0=0;p=(float *)malloc(sizeof(float *)*4);p0=(float *)malloc(sizeof(floa

43、t)*4*4);for(i=1;i4;i+) pi=pi-1+4; for(i=0;i4;i+)for(j=0;j4;j+)if(i=j) pij=pow(10,6);elsepij=0;xt=(float *)malloc(sizeof(float *)*1);xt0=(float *)malloc(sizeof(float)*1*4);for(i=1;i1;i+) xti=xti-1+4; for(i=0;i4;i+)xt0i=0; m1=(float *)malloc(sizeof(float *)*1);m10=(float *)malloc(sizeof(float)*1*4);fo

44、r(i=1;i1;i+) m1i=m1i-1+4; for(i=0;i4;i+)m10i=0; m2=(float *)malloc(sizeof(float *)*1);m20=(float *)malloc(sizeof(float)*1*1);for(i=1;i1;i+) m2i=m2i-1+1; for(i=0;i1;i+)m20i=0;m3=(float *)malloc(sizeof(float *)*1);m30=(float *)malloc(sizeof(float)*1*1);for(i=1;i1;i+) m3i=m3i-1+1; for(i=0;i1;i+)m30i=0;

45、 m4=(float *)malloc(sizeof(float *)*4);m40=(float *)malloc(sizeof(float)*4*1);for(i=1;i4;i+) m4i=m4i-1+1; for(i=0;i4;i+)m4i0=0;m5=(float *)malloc(sizeof(float *)*4);m50=(float *)malloc(sizeof(float)*4*4);for(i=1;i4;i+) m5i=m5i-1+4; for(i=0;i4;i+) for(j=0;j4;j+)m5ij=0; if(fpdu=fopen(u.dat,r)=NULL) printf(Cant open file u.datn); exit(0); i=0; do fscanf(fpdu,%g,&ui); i+; while(!feof(fpdu)&i500); fclose(fpdu); if(fpd