lso.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];
	int i, j, k, count;
	double y1[4][1],y[510],y2[1][1];
	double t1[4][1],t2[1][4],t3[1][1],t5[4][4],t6[4][4];
	double o[4][1], o1[500][4], w[4][1], w1[1][4];
	double p[4][4], km[4][1], er[4][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];

	y1[0][0] = -1.5;
	y1[1][0] = 0.7;
	y1[2][0] = 1.0;
	y1[3][0] = 0.5;
	for( i = 0; i < 4; i++ )	o[i][0] = 0;
	for( i = 0; i < 4; i++ )
		for( j = 0; j < 4;j++ )
		{
			if( i == j )  p[i][j] = 1000000;
			else	p[i][j] = 0;
		}
	for( k = 0; 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];

		brmul(p,w,4,4,1,t1);
		brmul(w1,p,1,4,4,t2);
		brmul(t2,w,1,4,1,t3);
		for( i = 0; i < 4; i++ ) km[i][0] = t1[i][0]/(1+t3[0][0]);
		brmul(w1,o,1,4,1,y2);
		for( i = 0; i < 4; i++ ) o[i][0] = o[i][0] + km[i][0]*(y[k+2] - y2[0][0]);
		brmul(t1,w1,4,1,4,t5);
		brmul(t5,p,4,4,4,t6);
		for( i = 0; i < 4; i++ )
			for( j = 0; j < 4; j++ ) p[i][j] = p[i][j] - t6[i][j]/(1+t3[0][0]);	
		for ( i = 0; i < 4; i++ ) o1[k][i] = o[i][0]; 

		maxe = 0;
		for( i = 0; i < 4; i++ )  
		{	
			er[i][0] = fabs((o[i][0] - y1[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\n"); 
	for( i = 0; i < count+1; i++ )
	{
		fprintf(fp,"%f\t%f\t%f\t%f\n",o1[i][0],o1[i][1],o1[i][2],o1[i][3]); 
	}
	printf("the output result in data3.txt\n");
	fclose(fp);

	return 0;  
}