📄 aicdlg.cpp
字号:
}
else
{
cs.Format(cs+"\nb%d=%.3f",i+1,sita1[j+2*i]);
}
}
}
m_cStaticTextREFm.SetDrawText(cs);
cs.Format("辨识参数：");
for(j=0;j<2;j++)
{
for(i=0;i<na0;i++)
{
if(j==0)
{
cs.Format(cs+"\na%d=%.6f",i+1,sita0[j+2*i]);
}
else
{
cs.Format(cs+"\nb%d=%.6f",i+1,sita0[j+2*i]);
}
}
}
m_cStaticTextREFi.SetDrawText(cs);
cs.Format("辨识误差：");
for(j=0;j<2;j++)
{
for(i=0;i<na0;i++)
{
if(j==0)
{
cs.Format(cs+"\na%d=%.6f",i+1,fabs(sita0[j+2*i]-sita1[j+2*i]));
}
else
{
cs.Format(cs+"\nb%d=%.6f",i+1,fabs(sita0[j+2*i]-sita1[j+2*i]));
}
}
}
m_cStaticTextREFe.SetDrawText(cs);
delete m_tempInt;
m_tempInt=new int[na];
for(i=0;i<na;i++)
{
m_tempInt[i]=int(AIC[i]*20000);
}
m_cStaticGraph.SetDrawData(na,m_tempInt,1,1);
m_cStaticGraph.RedrawWindow();
RedrawWindow();
}
void CAICDlg::brmul(double *a, double *b, int m, int n, int k, double *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;
}
void CAICDlg::brinv(double *a, int n)
{
int *is,*js,i,j,k,l,u,v;
double d,p;
is=(int *)malloc(n*sizeof(int));
js=(int *)malloc(n*sizeof(int));
for (k=0; k<=n-1; k++)
{ d=0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{ l=i*n+j; p=fabs(a[l]);
if (p>d) { d=p; is[k]=i; js[k]=j;}
}
if (d+1.0==1.0)
{ free(is); free(js);
return;
}
if (is[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=is[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (js[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+js[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
l=k*n+k;
a[l]=1.0/a[l];
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=k*n+j; a[u]=a[u]*a[l];}
for (i=0; i<=n-1; i++)
if (i!=k)
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
}
for (i=0; i<=n-1; i++)
if (i!=k)
{ u=i*n+k; a[u]=-a[u]*a[l];}
}
for (k=n-1; k>=0; k--)
{ if (js[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=js[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (is[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+is[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
}
free(is); free(js);
return;
}
void CAICDlg::junyun(int n, double x[], double y[])
{
int i,b;
unsigned long c;
c=4294967;
for (i=1;i<=n+1;i++)
{
b=(int)(A*x[i-1])/c;
x[i]=(A*x[i-1])-b*c;
y[i-1]=x[i-1]/c;
}
return;
}
void CAICDlg::WhiteNoise(int n, double mea, double var)
{
// FILE *fp,*fp3;
FILE *fp4;
int N=n+1;
int i;
double *inta1,*inta2;
double *x1,*x2,*y1,*y2;
inta1=(double*)malloc(N*sizeof(double));
inta2=(double*)malloc(N*sizeof(double));
x1=(double*)malloc((N+1)*sizeof(double));
x2=(double*)malloc((N+1)*sizeof(double));
y1=(double*)malloc(N*sizeof(double));
y2=(double*)malloc(N*sizeof(double));
x1[0]=11.0;
x2[0]=13.0;
junyun(N-1,x2,y2);
junyun(N-1,x1,y1);
// if((fp=fopen("result.txt","w"))==NULL)
// {
// MessageBox("cannot open the file :result.txt");
// return;
// }
// if((fp3=fopen("WhiteNoise1.txt","w"))==NULL)
// {
// MessageBox("cannot open the file :WhiteNoise1.dat\n");
// return;
// }
if((fp4=fopen("WhiteNoise2.txt","w"))==NULL)
{
MessageBox("cannot open the file :WhiteNoise2.dat\n");
return;
}
for (i=0;i<N-1;i++)
{
inta1[i]=sqrt(-2*log(y1[i]))*cos(2*PAI*y2[i])*sqrt(var)+mea;
inta2[i]=sqrt(-2*log(y1[i]))*sin(2*PAI*y2[i])*sqrt(var)+mea;
// if(inta1[i]>0) fprintf(fp3," %3.12f ",inta1[i]);
// else fprintf(fp3,"%3.12f ",inta1[i]);
// if (((i+1)%10)==0) fprintf(fp3,"\n");
if(inta2[i]>0) fprintf(fp4," %3.12f ",inta2[i]);
else fprintf(fp4,"%3.12f ",inta2[i]);
if (((i+1)%10)==0) fprintf(fp4,"\n");
}
/* fprintf(fp,"\n第一组正态分布白噪声函数：\n");
for (i=0;i<N-1;i++)
{
if(inta1[i]>0) fprintf(fp," %3.12f ",inta1[i]);
else fprintf(fp,"%3.12f ",inta1[i]);
if (((i+1)%10)==0) fprintf(fp,"\n");
}
fprintf(fp,"\n第二组正态分布白噪声函数：\n");
for (i=0;i<N-1;i++)
{
if(inta2[i]>0) fprintf(fp," %3.12f ",inta2[i]);
else fprintf(fp,"%3.12f ",inta2[i]);
if (((i+1)%10)==0) fprintf(fp,"\n");
}
*/
// fclose(fp);
// fclose(fp3);
fclose(fp4);
free(inta1);free(inta2);free(x1);free(x2);free(y1);free(y2);
return;
}
void CAICDlg::MSerials(int n)
{
FILE *fp;
int i,j,Num;
int m_nAmpMserial=1;
Num=(int)pow(2,n)-1;
int *yita,*u;
u=(int *)malloc(Num*sizeof(int));
yita=(int *)malloc(Num*sizeof(int));
if((fp=fopen("MSerials.txt","w"))==NULL)
{
MessageBox("ERROR");
return;
}
j=0;
for(i=0;i<n;i++)
{
yita[i]=1;
u[i]=yita[i]*m_nAmpMserial;
fprintf(fp,"%3d ",u[i]);
j++;
}
for(i=n;i<Num;i++)
{
switch( n )
{
case 2:
yita[i]=yita[i-n]^yita[i-1];
break;
case 3 :
yita[i]=yita[i-n]^yita[i-2];
break;
case 4 :
yita[i]=yita[i-n]^yita[i-3];
break;
case 5:
yita[i]=yita[i-n]^yita[i-3];
break;
case 6 :
yita[i]=yita[i-n]^yita[i-5];
break;
case 7 :
yita[i]=yita[i-n]^yita[i-4];
break;
case 8:
yita[i]=yita[i-n]^yita[i-7]^yita[i-6]^yita[i-1];
break;
case 9 :
yita[i]=yita[i-n]^yita[i-5];
break;
case 10 :
yita[i]=yita[i-n]^yita[i-7];
break;
case 11:
yita[i]=yita[i-n]^yita[i-9];
break;
case 12 :
yita[i]=yita[i-n]^yita[i-8]^yita[i-5]^yita[i-1];
break;
case 13 :
yita[i]=yita[i-n]^yita[i-8]^yita[i-6]^yita[i-1];
break;
}
u[i]=(2*yita[i]-1)*m_nAmpMserial;
fprintf(fp,"%3d ",u[i]);
j++;
if(j==20)
{
fprintf(fp,"\n");
j=0;
}
}
free(yita);free(u);
fclose(fp);
}
void CAICDlg::OnButtonOk()
{
// TODO: Add your control notification handler code here
UpdateData(TRUE);
if(m_nN>m_nMLen)
{
CString cs;
cs.Format("数据长度的值必须小于M序列长度，请输入一个小于%d的数值!",m_nMLen);
MessageBox(cs);
m_ctrN.SetFocus();
return;
}
MSerials(m_nMn);
WhiteNoise(m_nN, m_dNmea, m_dNvar);
switch(m_nIDMethods)
{
case 0:
aic(m_nN,m_nMLoc);
break;
case 1:
Ftest(m_nN,m_nMLoc);
break;
case 2:
OnlyAic(m_nN,m_nMLoc);
break;
}
}
void CAICDlg::OnChangeEDITMn()
{
// TODO: If this is a RICHEDIT control, the control will not
// send this notification unless you override the CDialog::OnInitDialog()
// function and call CRichEditCtrl().SetEventMask()
// with the ENM_CHANGE flag ORed into the mask.
// TODO: Add your control notification handler code here
UpdateData(TRUE);
m_nMLen=(int)pow(2,m_nMn)-1;
UpdateData(FALSE);
m_ctrMLoc.SetBuddy(&m_editMLoc);
if(m_nMLen<m_nN)
{
MessageBox("M序列长度小于白噪声长度，请重新输入!");
m_ctrMLoc.SetPos(1);
m_ctrMLoc.SetRange(1,1);
return;
}
m_ctrMLoc.SetRange(1,m_nMLen-m_nN);
}
void CAICDlg::OnChangeEDITDataN()
{
// TODO: If this is a RICHEDIT control, the control will not
// send this notification unless you override the CDialog::OnInitDialog()
// function and call CRichEditCtrl().SetEventMask()
// with the ENM_CHANGE flag ORed into the mask.
// TODO: Add your control notification handler code here
UpdateData(TRUE);
m_ctrMLoc.SetBuddy(&m_editMLoc);
if(m_nMLen<m_nN)
{
MessageBox("M序列长度小于数据组数，请重新输入!");
m_ctrMLoc.SetPos(1);
m_ctrMLoc.SetRange(1,1);
return;
}
m_ctrMLoc.SetPos(1);
m_ctrMLoc.SetRange(1,m_nMLen-m_nN);
}
void CAICDlg::WriteData2File(char *filename, double *a, int Length, int column)
{
FILE *fp;
int i,j;
char fileName[30];
strcpy(fileName,filename);
strcat(fileName,".txt");
if((fp=fopen(fileName,"w"))==NULL) return;
for(i=0;i<Length;i++)
{
for(j=0;j<column;j++)
{
fprintf(fp,"%12.6f ",a[i*column+j]);
}
fprintf(fp,"\n");
}
fclose(fp);
}
void CAICDlg::Ftest(int n, int lu)
{
//本程序最大能算到10阶。
int N=n-10;
FILE *fp1,*fp2,*fp4;
int i,j,na,nb,na0,nb0;
double t[10],v[10],sigama[1];
double *z,*u,*e,*X1T,*Y,*Ye,*e1,*X1,*X2T,*X2,*XTX_inv;
double sita0[20],sita1[20],sita2[2],*sita_temp;
double *X1TX1,*X1TY,X2TX2[4],*X2TX1,*X1TX2;
double *B1,*AA,B2[4],BB[4],*AX2T,*BX2T;
double *CC,*temp1,*temp2;
na=1;nb=1;t[0]=0;
z=(double*)calloc(N+10,sizeof(double));
u=(double*)calloc(N+10,sizeof(double));
e=(double*)calloc(N+10,sizeof(double));
X1T=(double*)calloc(20*N,sizeof(double));
X2T=(double*)calloc(2*N,sizeof(double));
X2=(double*)calloc(N*2,sizeof(double));
Y=(double*)calloc(N,sizeof(double));
Ye=(double*)calloc(N,sizeof(double));
e1=(double*)calloc(N,sizeof(double));
BX2T=(double*)calloc(2*N,sizeof(double));
if((fp1=fopen("MSerials.txt","r"))==NULL)return;
if((fp2=fopen("WhiteNoise2.txt","r"))==NULL)return;
//读入输入u[]和白噪声e[]
for(i=0;i<N+10;i++)
{
fscanf(fp1,"%lf ",&u[i]);
}
for(i=0;i<N+10;i++)
{
fscanf(fp2,"%lf ",&e[i]);
}
fclose(fp1);
fclose(fp2);
//用下式计算理论输出z[]
//z(k)-1.185z(k-1)+0.814z(k-2)-0.518z(k-3)+0.349z(k-4)-0.117z(k-5)
//=1.08u(k-1)-0.745u(k-2)+0.475u(k-3)-0.253u(k-4)+0.123u(k-5)+e(k)
z[0]=e[0];
z[1]=1.185*z[0]+1.08*u[0]+e[1];
z[2]=1.185*z[1]-0.814*z[0]+1.08*u[1]-0.745*u[0]+e[2];
z[3]=1.185*z[2]-0.814*z[1]+0.518*z[0]+1.08*u[2]-0.745*u[1]+0.475*u[0]+e[3];
z[4]=1.185*z[3]-0.814*z[2]+0.518*z[1]-0.349*z[0]+1.08*u[3]-0.745*u[2]+0.475*u[1]-0.253*u[0]+e[4];
for(i=5;i<N+10;i++)
{
z[i]=1.185*z[i-1]-0.814*z[i-2]+0.518*z[i-3]-0.349*z[i-4]+0.117*z[i-5]\
+1.08*u[i-1]-0.745*u[i-2]+0.475*u[i-3]-0.253*u[i-4]+0.123*u[i-5]+e[i];
}
WriteData2File("u_f",u,N,1);
WriteData2File("e_f",e,N,1);
WriteData2File("z_f",z,N,1);
///////////////
for(i=0;i<N;i++)
{
Y[i]=z[10+i];
}
//na=nb=1;
//赋值给X1T矩阵
for(i=0;i<2;i++)
{
for(j=0;j<N;j++)
{
if(i==0)
{
X1T[i*N+j]=-z[9+j];
}
else
{
X1T[i*N+j]=u[9+j];
}
}
}
//将X1T矩阵转置得X1矩阵。
X1=(double*)calloc(N*2,sizeof(double));
for(i=0;i<2;i++)
{
for(j=0;j<N;j++)
{
X1[j*2+i]=X1T[i*N+j];
}
}
//计算X1TX1，并求逆。
X1TX1=(double*)calloc(2*2,sizeof(double));
brmul(X1T,X1,2,N,2,X1TX1);
brinv(X1TX1,2);
//求参数估计值sita
X1TY=(double*)calloc(2,sizeof(double));
brmul(X1T,Y,2,N,1,X1TY);
brmul(X1TX1,X1TY,2,2,1,sita1);
free(X1TY);
//求预测误差矩阵e1。
brmul(X1,sita1,N,2,1,Ye);
for(i=0;i<N;i++)
{
e1[i]=Y[i]-Ye[i];
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -