ls.c

最小二乘一次完成算法

#include <stdio.h>
#include <math.h>
#include "brmul.c"
#include "brinv.c"
void main()
{  
	int i,j;
	double a1=-1.5,a2=0.7,b1=1.0,b2=0.5;
	double u[600]={0},v[600]={0};//u为输入，v为干扰
	double g[51][1]={0},y[500][1]={0};//g为系数,y2为输出y1[600],
	double w[500][51]={0},wt[51][500]={0};//w为数据阵,wt为w的转置
	double p[51][51]={0},q[51][500]={0};//p为wt*w,q为p的逆乘wt
	double b[51][1]={0};//所估计的参数
	double z[510][1]={0};//z为另一组输出
	double m[500][4],mt[4][500];//m为w，mt为m的转置
	double p2[4][4]={0},q2[4][500]={0};//p2 为mt*m,q2为p的转置
	double c[4][1]={0};
	double z2[500][1]={0};
	double J=0;//误差平方和
	FILE *fp1,*fp2,*fp3,*fp4,*fp5;
	fp1=fopen("m.txt","r");
	fp2=fopen("whitenoise.txt","r");
	fp3=fopen("估计的参数g.txt","w");
	fp4=fopen("真值.txt","w");
	fp5=fopen("估计的参数.txt","w");
	for(i=0;i<600;i++)
	{	 
		fscanf(fp1,"%lf",&u[i]);
		fscanf(fp2,"%lf",&v[i]);
	}
	fclose(fp1);
	fclose(fp2);
	for(i=0;i<51;i++)
	{   
		g[i][0]=0.25*exp(-0.5*i)-exp(-0.25*i)+0.75*exp(-i/6);
	    fprintf(fp4,"%lf     ",g[i][0]);
	}
	fclose(fp4);
	for(i=0;i<500;i++)
	{  
		for(j=0;j<51;j++)
		{   
			w[i][j]=u[i+50-j];  
			wt[j][i]=w[i][j];
		}
	}
	brmul(w,g,500,51,1,y);
	for(i=0;i<500;i++)
	  y[i][0]=y[i][0]+v[i];
	brmul(wt,w,51,500,51,p);
	brinv(p,51);
	brmul(p,wt,51,51,500,q);
	brmul(q,y,51,500,1,b);
	fprintf(fp3,"估计参数\n");
	for (j=0;j<51;j++)
	{   
		J=J+b[i][0]-g[i][0];
		fprintf(fp3,"%10f   ",b[j][0]);
	}
	fprintf(fp3,"\n");
	fprintf(fp3,"估计值与真值之差\n");
	for(i=0;i<51;i++)
		fprintf(fp3,"%10f    ",b[i][0]-g[i][0]);
	fprintf(fp3,"\n");
	fprintf(fp3,"误差平方和：");
	fprintf(fp3,"%lf   ",J);
	printf("误差平方和J=%lf   ",J);
	printf("\n");
	fclose(fp3);
	z[0][0]=v[0];
	z[1][0]=v[1];
	for(i=2;i<510;i++)
		z[i][0]=-a1*z[i-1][0]-a2*z[i-2][0]+b1*u[i-1]+b2*u[i-2]+v[i];

	for(i=0;i<500;i++)
		for (j=0;j<2;j++)
		{  
			m[i][j]=-z[i+1-j][0];
			m[i][j+2]=u[i+1-j];
		}
	for(i=0;i<4;i++)
		for(j=0;j<500;j++)
			mt[i][j]=m[j][i];
	for(i=0;i<500;i++)
		z2[i][0]=z[i+2][0];
	brmul(mt,m,4,500,4,p2);
	brinv(p2,4);
	brmul(p2,mt,4,4,500,q2);
	brmul(q2,z2,4,500,1,c);
	for(i=0;i<4;i++)
	{ 
		printf("%lf   ",c[i][0]);
		fprintf(fp5,"%lf\n",c[i][0]);
	}
}