sy2.cpp

系统辨识高效程序

#include <iostream.h>
#include <math.h>
#include <stdio.h>
#include <conio.h>
#define N 534
#define L 500
#define Nm 5
#define a 1
#define T0 1
#define sigma 1.0
#define Nv 12
#define mu 1

void mulmat(double* aa,double* b,double* c,int m,int n,int l)
{  int i,j,k;
   for(i=0;i<m;i++)
      for(j=0;j<l;j++)
	aa[i*l+j]=0;
   for(i=0;i<m;i++)
      for(j=0;j<l;j++)
       for(k=0;k<n;k++)
	  aa[i*l+j]+=b[i*n+k]*c[k*l+j];
}
void addmat(double*aa,double* b,double* c,int m,int n)
{  int i,j;
   for(i=0;i<m;i++)
      for(j=0;j<n;j++)
	aa[i*n+j]=0;
   for(i=0;i<m;i++)
      for(j=0;j<n;j++)
	  aa[i*n+j]=b[i*n+j]+c[i*n+j];

}
void kmulmat(double* aa,double* b,double k,int m,int n)
{  int i,j;
double xx,yy;
   for(i=0;i<m;i++)
      for(j=0;j<n;j++)
	aa[i*n+j]=0;
   for(i=0;i<m;i++)
      for(j=0;j<n;j++)
	 aa[i*n+j]=b[i*n+j]*k;
}
void submat(double* aa,double* b,double* c,int m,int n)
{  int i,j,k;
   double d[10*10];
   kmulmat(d,c,-1,m,n);
   addmat(aa,b,d,m,n);
}

void M_serial(double* u)
    { int i,k,m[Nm]={1,0,1,0,1};
//      FILE *fp;
//      fp=fopen("\\data.txt","a");
//      fprintf(fp,"%s","\n%the output of M_serial(u):\n");
      for(i=0;i<L;i++)
       { m[0]=(m[4]+m[1])%2;
	 for(k=Nm-1;k>0;k--) m[k]=m[k-1];
	 if(m[0]==0)   u[i]=1*a;
	 if(m[0]==1)   u[i]=-1*a;
//	 fprintf(fp,"%1.0f",u[i]);
//	 fputc(' ',fp);
//	 if(i%7==6) fprintf(fp,"%s","...\n");
       }
//      fclose(fp);
    }
void rand_N(double* u)
    { double kesi,sum;
      long x0,M;
      int i,k,A;
      M=32768;
      A=179;
      x0=11;
//      FILE *fp;
//      fp=fopen("\\data.txt","a");
//      fprintf(fp,"%s","\n%the output of rand_N(v):\n");
      for(k=0;k<L;k++)
	{ sum=0;
	  for(i=0;i<Nv;i++)
	   { x0=A*x0;
	     x0%=M;
	     kesi=double(x0)/M;
	     sum+=kesi;
	   }
      u[k]=0+sigma*(sum-6.0);
//      fprintf(fp,"%5.8f",u[k]);
//      fputc(' ',fp);
//      if(k%7==6)      fprintf(fp,"%s","...\n");
	}
//      fclose(fp);
    }
void gsout(double* u,double* z,double* v)
     {  double y[N],e[N];
	int k;
//	FILE *fp;
//	fp=fopen("\\data.txt","a");
//	fprintf(fp,"%s","\n%the output of gsout(G):\n");
	y[0]=0;e[0]=v[0];
	y[1]=u[0];e[1]=1.5*e[0]+v[1];
	for(k=2;k<L;k++)
	   {  y[k]=1.5*y[k-1]-0.7*y[k-2]+u[k-1]+0.5*u[k-2];
	      e[k]=1.5*e[k-1]-0.7*e[k-2]+v[k];
	   }
	for(k=0;k<L;k++)
	   {  z[k]=y[k]+e[k];
//	      fprintf(fp,"%5.8f",z[k]);
//	      fputc(' ',fp);
//	      if(k%7==6)  fprintf(fp,"%s","...\n");
	   }
//	fclose(fp);
     }
void makeh(double* h,double* z,double* u,int n,int k)
    { int i;
      for(i=0;i<n;i++)
	if(k-i-1>=0)
	{  h[i]=-z[k-i-1];
	   h[i+n]=u[k-i-1];
	}
	else
	{  h[i]=0;
	   h[i+n]=0;
	}
    }

void RFF(double* z,double* u,int n,double* sita,double* J)
    {int i,j,k;
     double P[400],K[20],h[20],eye[400],t1[400],t2[400],temp;
     for(i=0;i<2*n;i++)
	{ sita[i]=0.00001;
	  h[i]=0;
	  for(j=0;j<2*n;j++)
	      if(i==j)
	      {eye[i*2*n+j]=1;
	       P[i*2*n+j]=1000000;/*P(0)*/
	      }
	      else
	      {eye[i*2*n+j]=0;
	       P[i*2*n+j]=0;/*P(0)*/
	      }
	}
    (*J)=0;
     for(k=0;k<L;k++)
	{makeh(h,z,u,n,k);
	 mulmat(t1,h,sita,1,2*n,2*n);
	 (*J)+=(z[k]-(*t1))*(z[k]-(*t1));}
     makeh(h,z,u,n,0);
     for(k=1;k<L;k++)
	{  makeh(h,z,u,n,k);
	   /*-----J-----*/
	   mulmat(t1,h,sita,1,2*n,1);
	   temp=z[k]-(*t1);
	   mulmat(t1,h,P,1,2*n,2*n);
	   mulmat(t2,t1,h,1,2*n,1);
	   (*J)+=mu*temp*temp/(*t2+mu);
	   /*-----K-----*/
	   mulmat(t1,h,P,1,2*n,2*n);
	   mulmat(t2,t1,h,1,2*n,1);
	   temp=1.0/(*t2+mu);
	   mulmat(t1,P,h,2*n,2*n,1);
	   kmulmat(K,t1,temp,2*n,1);
	   /*-----P-----*/
	   mulmat(t1,K,h,2*n,1,2*n);
	   submat(t2,eye,t1,2*n,2*n);
	   mulmat(t1,t2,P,2*n,2*n,2*n);
	   kmulmat(P,t1,1.0/mu,2*n,2*n);
	   /*----sita----*/
	   mulmat(t1,h,sita,1,2*n,1);
	   temp=z[k]-(*t1);
	   kmulmat(t1,K,temp,2*n,1);
	   addmat(t2,sita,t1,2*n,1);
	   for(i=0;i<2*n;i++) sita[i]=t2[i];
	}
    }
void jamtosignal()
    {   double eita,sgmy,a0,a1,a2,b0,b1,b2;
	a0=1;	a1=-1.5;a2=0.7;
	b0=0;	b1=1;	b2=0.5;
	sgmy=pow(b0-b1,2)/(a0-a1);
	sgmy+=pow(b1*b1-pow(b2*a1/a0,2),2)/(a0*a0-pow(a2*a2/a0,2));
	sgmy+=b2*b2/a0;
	sgmy/=a0;
	eita=sqrt(sigma*sigma/sgmy);
	FILE *fp;
	fp=fopen("\\data.txt","a");
	fprintf(fp,"%s","\n%the output of jamtosignal(eita):\n");
	fprintf(fp,"%5.8f",eita);
	fclose(fp);
    }

void index(double* sita,int n)
    {int i;
     double k,ke,t1,t2,delta1,delta2,deltak,sita1[20]={-1.5,0.7,1,0.5,0,0,0,0,0,0,0,0};
     delta1=0;t1=0;t2=0;
     for(i=0;i<2*n;i++)
     {delta1+=pow((sita[i]-sita1[i])/sita1[i],2);
      t1+=pow(sita[i]-sita1[i],2);
      t2+=pow(sita1[i],2);
     }
     delta1=sqrt(delta1);
     delta2=sqrt(t1/t2);
     k=(1+0.5)/(1-1.5+0.7);
     t1=0;t2=0;
     for(i=0;i<n;i++)
     { t1+=sita[i+n];
       t2+=sita[i];
     }
     ke=t1/(1+t2);
     deltak=sqrt(fabs(k-ke)/k);
     FILE *fp;
     fp=fopen("\\data.txt","a");
     fprintf(fp,"%s","\n%the output of index(delta1):\n");
     fprintf(fp,"%5.8f",delta1);
     fprintf(fp,"%s","\n%the output of index(delta2):\n");
     fprintf(fp,"%5.8f",delta2);
     fprintf(fp,"%s","\n%the output of index(deltak):\n");
     fprintf(fp,"%5.8f",deltak);
     fprintf(fp,"%s","\n%the output of index(ke):\n");
     fprintf(fp,"%5.8f",ke);

     fclose(fp);
    }

void main()
    {   double u[N],z[N],v[N],tf,t,J1,J2;
	double sita[20]={0,0,0,0,0,0,0,0,0,0,0},tempsita[20],ramida;
	int i,n;
	FILE *fp;
	fp=fopen("\\data.txt","w+");
	fprintf(fp,"%s","\n%the output of pragramme:\n");
	fprintf(fp,"%s","\n%周超 200328014604111\n");
	fclose(fp);
	switch(L){
	case 100:tf=3.09;break;
	case 300:tf=3.03;break;
	case 500:tf=3.01;break;
	default: break;}
	rand_N(v);
	M_serial(u);
	gsout(u,z,v);
	J1=J2=1e10;n=1;
	t=tf+1;
	while(n<10 && t>tf)
	{  J1=J2;
	   for(i=0;i<2*n;i++) tempsita[i]=sita[i];
	   RFF(z,u,n,sita,&J2);
	   t=(J1-J2)/J2*(L-2*n)/2;
	   ramida=sqrt(J2/(L-2*n));
	   fp=fopen("\\data.txt","a");
	   fprintf(fp,"%s","\n%the output of RFF(J):\n");
	   fprintf(fp,"%s%d%s%5.8f%s%5.8f","  n=",n, "   J=",J2,"  t=",t);
	   fprintf(fp,"%s%5.8f","  ranmida=",ramida);
	   fprintf(fp,"%s","\n%the output of RFF(sita):\n");
	   for(i=0;i<2*n;i++)
	     { fprintf(fp,"%5.8f",sita[i]);
	       fputc(' ',fp); }
	    fputc('\n',fp);
	   fclose(fp);
	   n++;
	}
	jamtosignal();
	index(tempsita,2);
//	getch();

    }