递推增广最小二乘法(RELS)程序

1、8. 递推增广最小二乘法（RELS）程序该程序的输入信号为幅值为1的M序列，噪声源为均值为0，方差为0.1的高斯分布白噪声。该算法把噪声模型的参数归结入系统的参数向量中。程序的大体流程是：读入数据产生数据给P，W等向量置初值在大循环中由递推增广最小二乘法计算参数向量q 值当误差足够小时结束递推。逐次递推的q 值存入文件“kuodazuixiaoercheng.txt”中。递推公式如下： 其中：程序中模型多项式A，B，C的阶次na,nb,nc以及大循环的次数N均采用宏定义赋值，此处给值分别为na=2,nb=2,nc=2, N=850，根据实际情况和需要可改变这4个变量的值，在宏定义中修改即可。源

2、程序：#include <iostream.h>#include <fstream.h>/#include <stdlib.h>#include <math.h>#include<stdio.h>int brmul(double a,double b,int m,int n,int k,double c) /int m,n,k; /double a,b,c; int i,j,l,u; for (i=0; i<=m-1; i+) for (j=0; j<=k-1; j+) u=i*k+j; cu=0.0; for (l=0;

3、 l<=n-1; l+) cu=cu+ai*n+l*bl*k+j; return 0; #define na 2#define nb 2#define nc 2#define N 850void main()int k,i,j;double ana=-1.5,0.7,bnb=1.0,0.5,cnc=-1.0,0.2,z1000,u1000,e1000;ifstream fip1("m.txt");for(i=0;i<=N+na;i+)fip1>>ui;ifstream fip2("white.txt");for(i=0;i<

4、=N+na;i+)fip2>>ei;for(i=0;i<1000;i+)zi=0;/cout<<e0;ofstream fop("kuodazuixiaoercheng.txt");for(k=na;k<=N+na;k+)for(i=0;i<na;i+)zk=zk-ai*zk-i-1;for(i=0;i<nb;i+)zk=zk+bi*uk-i-1;for(i=0;i<nc;i+)zk=zk+ci*ek-i-1;zk=zk+ek;/cout<<zk;double ctana+nb+nc,wna+nb+nc;do

5、uble P(na+nb+nc)*(na+nb+nc),eaN+na+1;double a1=pow(10,3);for(i=0;i<na+nb+nc;i+)ctai=0;for(i=0;i<N+na+1;i+)eai=0;for(i=0;i<na;i+)wi=-zna-1-i;for(i=na;i<na+nb;i+)wi=una-1-(i-na);for(i=na+nb;i<na+nb+nc;i+)wi=eana-1-(i-na-nb);/cout<<wi<<endl;for(i=0;i<(na+nb+nc);i+)for(j=0;

6、j<(na+nb+nc);j+)if(i=j)Pi*(na+nb+nc)+j=a1*a1;elsePi*(na+nb+nc)+j=0;double Kna+nb+nc;double c2na+nb+nc,c11;for(k=0;k<N;k+)brmul(P,w,(na+nb+nc),(na+nb+nc),1,c2);brmul(w,c2,1,(na+nb+nc),1,c1);/cout<<c10;for(i=0;i<na+nb+nc;i+)Ki=c2i/(c10+1);/cout<<c2i<<endl;/cout<<Ki<

7、<" "/cout<<endl<<endl;double d1;brmul(w,cta,1,(na+nb+nc),1,d);for(i=0;i<(na+nb+nc);i+)ctai=ctai+Ki*(zk+na-d0);/cout<<ctai;double d1(na+nb+nc)*(na+nb+nc),d2(na+nb+nc)*(na+nb+nc);double P1(na+nb+nc)*(na+nb+nc);brmul(K,w,(na+nb+nc),1,(na+nb+nc),d1);for(i=0;i<na+nb+n

8、c;i+)for(j=0;j<na+nb+nc;j+)if(i=j)d2i*(na+nb+nc)+j=1-d1i*(na+nb+nc)+j;elsed2i*(na+nb+nc)+j=-d1i*(na+nb+nc)+j;brmul(d2,P,(na+nb+nc),(na+nb+nc),(na+nb+nc),P1);for(i=0;i<(na+nb+nc)*(na+nb+nc);i+)/for(j=0;j<na+nb+nc;j+)Pi=P1i;double f1;brmul(w,cta,1,(na+nb+nc),1,f);eak+na=zk+na-f0;for(i=0;i<

9、na;i+)wi=-zk+na-i;for(i=na;i<(na+nb);i+)wi=uk+na-(i-na);for(i=na+nb;i<(na+nb+nc);i+)wi=eak+na-(i-na-nb);for(i=0;i<na;i+) fop <<"a真实值:"<<ai<<"辩识值"<<ctai<<endl; for(i=na;i<(na+nb);i+) fop <<"b真实值:"<<bi-na<<"

10、辩识值"<<ctai<<endl; for(i=na+nb;i<(na+nb+nc);i+) fop <<"c真实值:"<<ci-na<<"辩识值"<<ctai<<endl;for(i=0;i<na;i+) cout<<ctai<<endl; /fop <<"a真实值:"<<ai<<"辩识值"<<ctai<<endl; for(i=na;i<(na+nb);i+) cout<<ctai<<endl; /fop <<"b真实值:"<<bi-na<<"辩识值"<<ctai<