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

8. 递推增广最小二乘法（RELS ）程序
该程序的输入信号为幅值为1的M 序列，噪声源为均值为0，方差为0.1的高斯分布白噪声。

该算法把噪声模型的参数归结入系统的参数向量中。

程序的大体流程是：读入数据—产生数据—给P ，W 等向量置初值—在大循环中由递推增广最小二乘法计算参数向量θ 值—当误差足够小时结束递推。

逐次递推的θ 值存入文件“kuodazuixiaoercheng.txt ”中。

递推公式如下：
]ˆ)()1([ˆˆT 011N a e N N N N N z θw K θθ-++=++
10T 001)]()()(1)[()(-++=N N N N N a e a e a e N w P w w P K
)()]([)1(T 01N N N a e N P w K I P +-=+
其中：
)1()([)(T 0+---=a a e n N z N z N w )1()(+-b n N u N u
)]1(ˆ)(ˆ+-c a a n N e N e
程序中模型多项式A ，B ，C 的阶次na,nb,nc 以及大循环的次数N 均采用宏定义赋值，此处给值分别为na=2,nb=2,nc=2, N =850，根据实际情况和需要可改变这4个变量的值，在宏定义中修改即可。

源程序：
#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; c[u]=0.0;
for (l=0; l<=n-1; l++)
c[u]=c[u]+a[i*n+l]*b[l*k+j];
}
return 0;
}
#define na 2
#define nb 2
#define nc 2
#define N 850
void main()
{
int k,i,j;
double a[na]={-1.5,0.7},b[nb]={1.0,0.5},c[nc]={-1.0,0.2},z[1000],u[1000],e[1000]; ifstream fip1("m.txt");
for(i=0;i<=N+na;i++)
fip1>>u[i];
ifstream fip2("white.txt");
for(i=0;i<=N+na;i++)
fip2>>e[i];
for(i=0;i<1000;i++)
z[i]=0;
//cout<<e[0];
ofstream fop("kuodazuixiaoercheng.txt");
for(k=na;k<=N+na;k++)
{for(i=0;i<na;i++)
z[k]=z[k]-a[i]*z[k-i-1];
for(i=0;i<nb;i++)
z[k]=z[k]+b[i]*u[k-i-1];
for(i=0;i<nc;i++)
z[k]=z[k]+c[i]*e[k-i-1];
z[k]=z[k]+e[k];
// cout<<z[k];
}
double cta[na+nb+nc],w[na+nb+nc];
double P[(na+nb+nc)*(na+nb+nc)],ea[N+na+1];
double a1=pow(10,3);
for(i=0;i<na+nb+nc;i++)
cta[i]=0;
for(i=0;i<N+na+1;i++)
ea[i]=0;
for(i=0;i<na;i++)
w[i]=-z[na-1-i];
for(i=na;i<na+nb;i++)
w[i]=u[na-1-(i-na)];
for(i=na+nb;i<na+nb+nc;i++)
{w[i]=ea[na-1-(i-na-nb)];
//cout<<w[i]<<endl;
}
for(i=0;i<(na+nb+nc);i++)
for(j=0;j<(na+nb+nc);j++)
{if(i==j)
P[i*(na+nb+nc)+j]=a1*a1;
else
P[i*(na+nb+nc)+j]=0;
}
double K[na+nb+nc];
double c2[na+nb+nc],c1[1];
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<<c1[0];
for(i=0;i<na+nb+nc;i++)
{K[i]=c2[i]/(c1[0]+1);
//cout<<c2[i]<<endl;
//cout<<K[i]<<" ";
}
//cout<<endl<<endl;
double d[1];
brmul(w,cta,1,(na+nb+nc),1,d);
for(i=0;i<(na+nb+nc);i++)
{cta[i]=cta[i]+K[i]*(z[k+na]-d[0]);
//cout<<cta[i];
}
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+nc;i++)
for(j=0;j<na+nb+nc;j++)
{if(i==j)
d2[i*(na+nb+nc)+j]=1-d1[i*(na+nb+nc)+j];
else
d2[i*(na+nb+nc)+j]=-d1[i*(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++)
P[i]=P1[i];
double f[1];
brmul(w,cta,1,(na+nb+nc),1,f);
ea[k+na]=z[k+na]-f[0];
for(i=0;i<na;i++)
w[i]=-z[k+na-i];
for(i=na;i<(na+nb);i++)
w[i]=u[k+na-(i-na)];
for(i=na+nb;i<(na+nb+nc);i++)
w[i]=ea[k+na-(i-na-nb)];
for(i=0;i<na;i++)
{
fop <<"a真实值:"<<a[i]<<"辩识值"<<cta[i]<<endl;
}
for(i=na;i<(na+nb);i++)
{
fop <<"b真实值:"<<b[i-na]<<"辩识值"<<cta[i]<<endl; }
for(i=na+nb;i<(na+nb+nc);i++)
{
fop <<"c真实值:"<<c[i-na]<<"辩识值"<<cta[i]<<endl; }
}
for(i=0;i<na;i++)
{cout<<cta[i]<<endl;
//fop <<"a真实值:"<<a[i]<<"辩识值"<<cta[i]<<endl;
}
for(i=na;i<(na+nb);i++)
{cout<<cta[i]<<endl;
//fop <<"b真实值:"<<b[i-na]<<"辩识值"<<cta[i]<<endl; }
for(i=na+nb;i<(na+nb+nc);i++)
{cout<<cta[i]<<endl;
//fop <<"c真实值:"<<c[i-na]<<"辩识值"<<cta[i]<<endl; }
}。