?? guji.cpp
字號:
/*******************************************************************/
/*-----------遞推最小二乘法估計及模型階次檢驗實驗------------------*/
/*******************************************************************/
#include <dos.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>
const int P=4,Np=(int)pow(2,P)-1,Len=600,realN=2; //P為M序列的階次，Np循環(huán)周期,realN理想階次
static int nL=2,L=150+nL*150,Nbeg=1,Nend=4,resultN;
static double *u=new double[Len]; //M序列
static double *z=new double[Len]; //有噪聲的輸出量
static double *v=new double[Len]; //白噪聲
static double *e=new double[Len]; //有色噪聲
static double *y=new double[Len]; //無噪聲的輸出量
static double *UU=new double[500*2*realN];
static double Sigma=0.5,a=1,miu=0.998; //Sigma噪聲標準差，a為M序列的幅度，miu遺忘因子
static double ab[4]={-1.5,0.7,1,0.5}; //真實參數
static double J[10],Lamda[10],Kj[10]; //損失函數、噪聲標準差、模型靜態(tài)增益估計
static double ZXB,deta1,deta2,detak; // 噪信比、參數估計平方相對偏差、參數估計平方根偏差、靜態(tài)增益估計相對偏差
static double TN[10],Q[4]; //定階法估計、二階辨識結果、
const double ta[3]={3.09,3.03,3.01}; //風險水平0.05時L＝100，300，500的閥值
/*********************以下是與矩陣相關的一些運算*************************/
/*----------------------構造一個單位陣----------------------------------*/
void CreateIMatrix(int n,double *im)
{
for(int i=0;i<n*n;i++)
{
if(int(i/n)==i%n)
*(im+i)=1;
else *(im+i)=0;
}
}
/*-----------------------矩陣加法-----------------------------------------*/
void AddMatrix(double *jia1,double *jia2,double *he,int row,int vol)
{
int n=row*vol;
for(int i=0;i<n;i++)
{
*(he+i)=*(jia1+i)+(*(jia2+i));
}
}
/*-----------------------矩陣數乘-----------------------------------------*/
void MulNumMatrix(double *one,double num,double *ji,int row,int vol)
{
int n=row*vol;
for(int i=0;i<n;i++)
{
*(ji+i)=num*(*(one+i));
}
}
/*------------------------矩陣乘法----------------------------------------*/
void MulMatrix(double *one,int row,int rv,int vol,double *two,double *ji)
{
int i,j,k;
double temp;
for(i=0;i<row;i++)
{
for(j=0;j<vol;j++)
{
temp=0;
for(k=0;k<rv;k++)
{
temp+=(*(one+rv*i+k))*(*(two+k*vol+j));
}
*(ji+i*vol+j)=temp;
}
}
}
/*********************以下是模型輸出************************/
/*-------------------利用位運算生成M序列-------------------*/
void CreateU(int L)
{
int i;
unsigned m=2; //給初始值u[0]=1
for(i=0;i<L;i++)
{
m=m>>1; //右移位
if(m&1)
u[i]=-a; //邏輯為1，幅度－a
else
u[i]=a; //邏輯為0，幅度a
if((m&9)==1||(m&9)==8)
m=m|16; //循環(huán)
}
}
/*--------------------生成白噪聲v-----------------------------*/
void CreateV(int L)
{
double M=32768,A=179,x0=11,ksai,x11;//乘同余法參數
int i,k;
x11=x0;
for(k=0;k<L;k++)
{
ksai=0;
for(i=1;i<12;i++)
{
x11=A*x11;
x11=fmod(x11,M); //生成u[0,1]均勻分布隨機數
ksai=ksai+float(x11/M);
}
v[k]=Sigma*(ksai-6.0); //生成正態(tài)分布的白噪聲
}
}
/*--------------------生成有色噪聲e-----------------------------*/
void CreateE(int L)
{
int k;
for(k=0;k<L;k++)
{
if(k==0)e[k]=v[k]; //零階系統
else
if(k==1)
e[k]=v[k]-ab[0]*e[k-1]; //一階系統
else
e[k]=v[k]-ab[0]*e[k-1]-ab[1]*e[k-2];//二階系統
}
}
/*----------------------傳遞函數的輸出Z(k)-------------------*/
void CountZ(int L,double *ab)
{
int k;
y[0]=0;
y[1]=u[0];
for(k=2;k<L;k++)
y[k]=-ab[0]*y[k-1]-ab[1]*y[k-2]+ab[2]*u[k-1]+ab[3]*u[k-2];//從方框圖中得出的模型表示
for(k=0;k<L;k++)
z[k]=y[k]+e[k]; //含有有色噪聲的輸出量
}
/**************************以下是RLS同F－test子程序***********************/
/*-----------------------遞推最小二乘法----------------------------------*/
void Fzhen(int n)
{
double *thita=new double[2*n];//辨識出的參數
double *hk=new double[2*n]; //已知數據
double *P1=new double[4*n*n],*P2=new double[4*n*n];//P1、P2分別為前一時刻和當前時刻的P矩陣
double *K=new double[2*n]; //增益矩陣K
double te,tm,*ti=new double[4*n*n],*temp;
int i,k;
J[n]=0; //損失函數
/*------------初始化thita和P-------------------*/
CreateIMatrix(2*n,P1);
MulNumMatrix(P1,1E+6,P1,2*n,2*n); //為收斂快，P＝10(6)*I
for(i=0;i<2*n;i++)
{
thita[i]=0.001; //thita的初始值很小
UU[i]=thita[i];
}
/*----------------遞推核心算法------------------*/
for(k=n;k<L;k++)
{
for(i=0;i<n;i++)
{
hk[i]=-z[k-i-1];
hk[i+n]=u[k-i-1]; //hk為已知的輸入輸出量
}
//計算K(k)
MulMatrix(hk,1,2*n,2*n,P1,ti);
MulMatrix(ti,1,2*n,1,hk,K);
tm=1/(K[0]+miu);
MulMatrix(P1,2*n,2*n,1,hk,ti);
MulNumMatrix(ti,tm,K,2*n,1);
//計算P(k)
MulMatrix(K,2*n,1,2*n,hk,ti);
MulNumMatrix(ti,-1,ti,2*n,2*n);
for(i=0;i<2*n;i++)
ti[i*2*n+i]+=1;
MulNumMatrix(ti,1/miu,ti,2*n,2*n);
MulMatrix(ti,2*n,2*n,2*n,P1,P2);
//計算thita(k)
MulMatrix(hk,1,2*n,1,thita,ti);
te=z[k]-ti[0];
MulNumMatrix(K,te,ti,2*n,1);
AddMatrix(thita,ti,thita,2*n,1);
//計算J(k)
J[n]=(J[n]+tm*te*te);
temp=P1;
P1=P2;
P2=temp;
}
if(n==realN)
{
for(i=0;i<2*n;i++)
Q[i]=thita[i]; //存儲理想階次下的辨識結果
}
Lamda[n]=sqrt(J[n]/(L-2*n));//估計噪聲標準差
/*--------------估計模型靜態(tài)增益----------------*/
te=0;tm=0;
for(i=0;i<n;i++)
{
te+=thita[i];
tm+=thita[i+n];
}
if(te==-1)
Kj[n]=1E+10;
else
Kj[n]=tm/(1+te);
/*------------保存遞推過程中的辨識參數數據-----------*/
FILE *file=fopen("最終結果.txt","w");
if((file=fopen("最終結果.txt","w"))==NULL)
{
printf("cannot open the file.\n");
return;
}
fprintf(file,"%f\t\t%f\t\t%f\t\t%f\t\t\n",Q[0],Q[1],Q[2],Q[3]);
fclose(file);
}
/*--------------------F-Test定階--------------------------------*/
bool AcceptN(int n)
{
if(n<1||n>8)return false;
TN[n]=(J[n]/J[n+1]-1)*(L-2*n-2)/2;//統計量
if(TN[n]<=ta[nL])return true; //接受H：n+1>=n0；階次估計值為n+1
else return false; //拒絕H：n+1>=n0
}
/*----------------------過程輸出方差的計算------------------------------*/
void fangcha(double *a,double *b,double *at,double *bt,int n)
{
int i,k=n-1;
for(i=0;i<n;i++) //具有周期性
{
at[k*n+i]=a[i];
bt[k*n+i]=b[i];
}
for(k=n-2;k>=0;k--)
{
for(i=0;i<=k;i++)
{
at[k*n+i]=at[(k+1)*n+i]-(at[(k+1)*n+k+1]*at[(k+1)*n+k+1-i])/at[(k+1)*n];
bt[k*n+i]=bt[(k+1)*n+i]-(bt[(k+1)*n+k+1]*at[(k+1)*n+k+1-i])/at[(k+1)*n];
}
}
}
/*--------------------------性能計算---------------------------*/
void CountG()
{
/*------------計算噪聲方差----------------*/
double Deltae=0,e0=0;
int i;
for(i=0;i<L;i++)
e0+=e[i];
e0/=L; //有色噪聲均值
for(i=0;i<L;i++)
Deltae+=pow((e[i]-e0),2); //有色噪聲方差
Deltae/=(L-1);
/*-----------計算過程輸出方差--------------*/
double Deltay=0,a[3],b[3];
a[2]=ab[1];
a[1]=ab[0];
a[0]=1;
b[0]=0;
b[1]=ab[2];
b[2]=ab[3];
int n=3;
double *at=new double[n*n],*bt=new double[n*n];
fangcha(a,b,at,bt,n);
for(i=0;i<n;i++)
Deltay+=pow(bt[i*n+i],2)/at[i*n+i];
Deltay=Deltay/a[0];
/*-----------計算噪信比---------------*/
ZXB=sqrt(Deltae/Deltay);
printf("系統輸出噪信比:%f\n",ZXB);
/*----------參數估計平方相對偏差------*/
deta1=0;
for(i=0;i<2*realN;i++)
deta1+=pow((ab[i]-Q[i])/ab[i],2);
deta1=sqrt(deta1);
printf("參數估計平方相對偏差:%f\n",deta1);
/*---------參數估計平方根偏差-----------*/
double temp=0;
deta2=0;
for(i=0;i<2*realN;i++)
{
deta2+=pow((ab[i]-Q[i]),2);
temp+=pow(ab[i],2);
}
deta2=sqrt(deta2/temp);
printf("參數估計平方根偏差:%f\n",deta2);
/*---------靜態(tài)增益估計相對偏差-----------*/
temp=0,Kj[0]=0;
for(i=0;i<realN;i++)
{
temp+=ab[i]; //ai之和
Kj[0]+=ab[i+realN]; //bi之和
}
if(temp==-1)
Kj[0]=1E+10;
else
Kj[0]=Kj[0]/(1+temp);
detak=sqrt(fabs((Kj[0]-Kj[realN])/Kj[0]));
printf("靜態(tài)增益估計相對偏差:%f\n",detak);
return;
}
/*-------------------------主函數-----------------------------*/
void main()
{
loop:
printf("----------遞推最小二乘法及模型階次檢驗實驗-------------\n");
printf("請輸入標準差值(范圍0～1):\n");
scanf("%lf",&Sigma);
getchar();
if(Sigma<0||Sigma>1) //標準差在[0,1]之間
{
printf("error.\n");
goto loop; //重新輸入sigma
}
printf("遺忘因子：\n");
scanf("%f",&miu);
getchar();
printf("數據長度：\n");
scanf("%d",&L);
getchar();
printf("開始階數：\n");
scanf("%d",&Nbeg);
getchar();
printf("終止階數：\n");
scanf("%d",&Nend);
getchar();
if(Nbeg>=Nend)resultN=-2;
CreateU(L); //生成M序列
CreateV(L); //生成輸出白噪聲v
CreateE(L); //生成有色噪聲e
CountZ(L,ab);//傳遞函數的輸出Z(k)
int i=Nbeg;
Fzhen(i);
for(i++;i<=Nend;i++) //循環(huán)確定理想階次
{
Fzhen(i);
if(AcceptN(i-1))
break;
}
printf("理想階次：%d\n",realN);
if(Nbeg>realN||i<realN)
Fzhen(realN);
printf("辨識結果:\n");
for(i=0;i<2*realN;i++)
{
printf("%f\t",Q[i]);
}
printf("\n");
printf("噪聲標準差估計：%f\n",Lamda[realN]);
printf("靜態(tài)增益估計：%f\n",Kj[realN]);
CountG();
if(i>Nend)resultN=-1;
else resultN=i-1;
printf("\n");
goto loop;
return;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -