els.c

递推增广最小二乘法的VC程序代码

#include "math.h"
#include "stdio.h"	 
#include "brmul.c"
#include "brinv.c"
int main()
{
	FILE *fp1,*fp2,*fp;
	int u[600]; double v[600], e[600];
	int i, j, k, count;
	double y1[6][1],y[510],y2[1][1];
	double t1[6][1],t2[1][6],t3[1][1],t5[6][6],t6[6][6];
	double o[6][1], o1[500][6], w[6][1], w1[1][6];
	double p[6][6],	p1[6][6], km[6][1], er[6][1];
	double maxe;
	t3[0][0] = 0.0;
	y2[0][0] = 0.0;

	fp1 = fopen("data1.txt","r");
	fp2 = fopen("data2.txt","r");
	for( i = 0; i < 600; i++ ) 
	{
		fscanf(fp1,"%d",&u[i]);//输入信号M序列
		fscanf(fp2,"%lf",&v[i]);//均值为0,方差为0.1的白噪声
	}
	fclose(fp1);
	fclose(fp2);

	y[0] = v[0];
	y[1] = v[1];
	for( i = 2; i < 510; i++ )
		y[i] = 1.5*y[i-1]-0.7*y[i-2]+u[i-1]+0.5*u[i-2]-v[i-1]+0.2*v[i-2]+v[i];

	y1[0][0] = -1.5;
	y1[1][0] = 0.7;
	y1[2][0] = 1.0;
	y1[3][0] = 0.5;
	y1[4][0] = -1.0;
	y1[5][0] = 0.2;
	for( i = 0; i < 6; i++ )	o[i][0] = 0;
	for( i = 0; i < 6; i++ )
		for( j = 0; j < 6;j++ )
		{
			if( i == j )  p[i][j] = 1000000;
			else	p[i][j] = 0;
		}
	e[1] = e[0] = 0;
	for( k = 2; k < 500; k++ )
	{
		for( i = 0; i < 2; i++ ) w[i][0] = w1[0][i] = -y[1-i+k];
		for( i = 2; i < 4; i++ ) w[i][0] = w1[0][i] = u[3-i+k];
		for( i = 4; i < 6; i++ ) w[i][0] = w1[0][i] = e[3-i+k];

		brmul(p,w,6,6,1,t1);
		brmul(w1,p,1,6,6,t2);
		brmul(t2,w,1,6,1,t3);
		for( i = 0; i < 6; i++ ) km[i][0] = t1[i][0]/(1+t3[0][0]);
		brmul(w1,o,1,6,1,y2);
		for( i = 0; i < 6; i++ ) o[i][0] = o[i][0] + km[i][0]*(y[k+2] - y2[0][0]);
		brmul(km,w1,6,1,6,t6);
	   	for( i = 0; i < 6; i++ )
			for( j = 0; j < 6; j++ )
			{
				if( i == j )  p1[i][j] = 1 - t6[i][j];
				else	p1[i][j] = -t6[i][j];
			}
		brmul(p1,p,6,6,6,t5);
		for( i = 0; i < 6; i++ )
			for( j = 0; j < 6; j++ )	  p[i][j] = t5[i][j];
		
		brmul(w1,o,1,6,1,y2);
		e[k] = y[k+2] - y2[0][0];

		for ( i = 0; i < 6; i++ ) o1[k][i] = o[i][0]; 

		maxe = 0;
		for( i = 0; i < 6; i++ )  
		{	
			er[i][0] = fabs((o[i][0] - y1[i][0]));
			maxe = max( er[i][0], maxe );
		 }
		count = k;
		if( maxe < 0.01 ) break; 
	}
	fp = fopen("Data3.txt","w");
	fprintf(fp,"a0\t\ta1\t\tb0\t\tb1\t\tc1\t\tc2\n"); 
	for( i = 2; i < count+1; i++ )
	{
		fprintf(fp,"%f\t%f\t%f\t%f\t%f\t%f\n",o1[i][0],o1[i][1],o1[i][2],o1[i][3],o1[i][4],o1[i][5]); 
	}
	fclose(fp);

	return 0;  
}