用于滤波，rls自适应滤波
#include
#include
#include
#define N 110000
float x[N];
float x1[110050];
float xout[N];
float w[50];
float p[50][50];
float u[50];
float kn[50];
float kn1[50];
float dn;
float e;
float outmiddle;
//阶数50//
//****************************************初始化程序用于计算w[][].p[]**********************************//
float initial ()
{
char i1,j1;
for(i1=0;i1<50;i1++)
{
w[i1]=0;
}
for(i1=0;i1<50;i1++)
{
for(j1=0;j1<50;j1++)
{
p[i1][j1]=0;
}
}
for(i1=0;i1<50;i1++)
{
p[i1][i1]=9999999.9;
}
return(0);
}
//****************************************************************************************************************//
//****************************************e[n]的求解**************************************************************//
float en(unsigned long int datanumber)
{
/*float result=0;
int i2,k2;
int j2;
j2=datanumber;
k2=datanumber+5;
e=0;
for(i2=0;i2<50;i2++)//计算w()*u()
{
e=e+w[i2]*u[i2];
}
dn=0;
for(i2=j2;i2<k2;i2++)
{
dn=dn+x[i2];
}
dn=dn/5;//给定dn为0.2
e=dn-e;
return(0);*/
}
//****************************************************************************************************************//
//****************************************k[n]的求解**************************************************************//
void knkn()
{
int i3,j3;
float knmiddle=0;
float middle[50];
float pluse=0;
for(i3=0;i3<50;i3++)
{
for(j3=0;j3<50;j3++)
{
knmiddle=knmiddle+p[i3][j3]*u[j3];
}
kn[i3]=knmiddle;
knmiddle=0;
}//计算p(n-1)*u(n)
for(i3=0;i3<50;i3++)
{
middle[i3]=kn[i3];
}
for(i3=0;i3<50;i3++)
{
pluse=pluse+middle[i3]*u[i3];
}//计算uh(n)*p(n-1)*u(n)
pluse=pluse+0.7;//衰减因数0.7
for(i3=0;i3<50;i3++)
{
kn[i3]=kn[i3]/pluse;
}
pluse=0;
}
//****************************************************************************************************************//
//*******************************************用于pn[]的更新*******************************************************//
float newpn()
{
float newpn[50][50];
float newpn1[50][50];
float newpn3[50][50];
float newpn2=0;
int i4,j4,k4;
for(i4=0;i4<50;i4++)
{
for(j4=0;j4<50;j4++)
{
newpn[i4][j4]=kn[i4]*u[j4];
newpn3[i4][j4]=p[i4][j4];
}
}//计算k(n)*u(n)
for(i4=0;i4<50;i4++)
{
for(j4=0;j4<50;j4++)
{
for(k4=0;k4<50;k4++)
{
newpn2=newpn2+newpn[i4][k4]*newpn3[k4][j4];
}
newpn1[i4][j4]=newpn2;
newpn2=0;
}
}//计算k()*u()*p()
for(i4=0;i4<50;i4++)
{
for(j4=0;j4<50;j4++)
{
p[i4][j4]=p[i4][j4]-newpn1[i4][j4];
p[i4][j4]=p[i4][j4]/0.7;//衰减因子是0.7
}
}
return(0);
}
//*****************************************************************************************************************//
//*************************************用于wn[]的更新**************************************************************//
float newwn()
{
char i5,j5;
for(i5=0;i5<50;i5++)
{
kn1[i5]=kn[i5];
}
for(i5=0;i5<50;i5++)
{
kn1[i5]=kn1[i5]*e;
w[i5]=w[i5]+kn1[i5];
}
return(0);
}
main()
{ //, f[M][N], g[M][N], a[M][M], P[N], K[M];
//float x1[200];
//float Knum = 0, Kden = 0, sum = 0, eP;
// int i, j, k, m, n;
float dn;
long int i=0, j=0,k=0,m=0,n1,n2,n3;
float result=0;
int i2,k2;
int j2;
char controller=1;
FILE * fpr, * fp_out;
int max[2][2]={1,2,3,4};
int max1[2]={1,2};
int max2[2]={2,3};
int maxmiddle;
//a[0][1] = 1;
if((fpr = fopen("wavewavedata.bin","rb")) == NULL)
{printf("cannot open file\n");
exit(0); }
for(i=0;i<N;i++)
{
fread(&x[i], 4, 1, fpr);
//printf("%f\n",x[i]);
}
/*for(i=0;i<49;i++)
{
x1[i]=x[110000-i-1];
}
for(i=49;i<110049;i++)
{
x1[i]=x[i-49];
}//得到要处理的110049个数据*/
fclose(fpr);
for(i=0;i<N;i++)
{
xout[i]=x[i];
}
if((fp_out = fopen("bbboutdata.bin","w")) == NULL)
{ printf("cannot open file\n");
return;
}
else printf("success fp_out = %d\n",fp_out);
for(j=0; j<N; j++)
{
fwrite(&xout[j],4, 1, fp_out);
//fread(&x1[i][j],4, 1, fpr);
//printf("%f\n",xout[j]);
}
fclose(fp_out);
getchar();
initial ();
for(m=0;m<N;m++)
{
//m=0;
n1=m;
n2=m+50;
k=49;
for(i=n1;i<n2;i++) //取出要处理的50个数据放于u[n]
{
u[k]=x1[i];
k--;
}
//printf("%f\n",u[0]);
while(controller)
{
//en(n3);
j2=n2-3;
k2=n2+2;
e=0;
for(i2=0;i2<50;i2++)//计算w()*u()
{
e=e+w[i2]*u[i2];
}
dn=0;
for(i2=j2;i2<k2;i2++)
{
dn=dn+x1[i2];
}
dn=dn/5;//给定dn为0.2
e=dn-e;
if(e<0.000000001)
{
controller=0;
/*printf("success!",m);
for(i=0;i<50;i++)
{
printf("%f\n",w[i]);
}
printf("success!",m);
for(i=0;i<50;i++)
{
printf("%f\n",kn[i]);
}*/
}
else
{
knkn();
newpn();
newwn();
}
}
/*for(j=0;j<50;j++)
{
printf("%f\n",w[j]);
}*/
outmiddle=0;
for(i=0;i<50;i++)
{
outmiddle=outmiddle+w[i]*u[i];
}
xout[m]=outmiddle;
//printf("output",m);
//printf("%f\n",xout[m]);
}
/*for(i=0;i<N;i++)
{
//xout[i]=x[i];
printf("%f\n",xout[i]);
}*/
/*if((fp_out = fopen("bboutdata.bin","w")) == NULL)
{ printf("cannot open file\n");
return;
}
else printf("success fp_out = %d\n",fp_out);
for(j=0; j<N; j++)
{
fwrite(&xout[j],4, 1, fp_out);
//fread(&x1[i][j],4, 1, fpr);
//printf("%f\n",xout[j]);
}
fclose(fp_out);
getchar();*/
}
扫一扫
250
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
