ergoergo的专栏

自适应积分算法，适用于对函数的积分。对在不同区间变化趋势不同的函数进行积分。算法收敛迅速。占用的程序空间很小。对数组的操作采用了堆栈的思想。C语言实现，易于在DSP，ARM，单片机上移植

最小二乘法用于参数估计，用于检测和参数估计.c语言实现。用六个方程估计3个参数。

#include <stdio.h> #include <stdlib.h> #include <math.h> #define N 110000 double x[N]; double xout[N]; double f0[N]; double g0[N]; double f00[N]; double g00[N]; double p0; double km; double am[1000]; double amchange[1000]; int m; double error; double error1; double amtry[20]={-1.561735,1.436875,-1.228909,1.114787,-1.286210,1.066792,-0.818686,0.850963,-0.777519,0.624041,-0.474253,0.475683,-0.404298,0.217569,-0.185316,0.154146,-0.082471,-0.043769,0.069651}; //****************************************初始化程序用于计算PO，f0,g0以及确定阶数m**********************************// double initial () { double z; int i1,j1; for(i1=0;i1<N;i1++) { f0[i1]=x[i1]; g0[i1]=x[i1]; p0=p0+x[i1]*x[i1]; } p0=p0/N; return(0); } //****************************************************************************************************************// //**************************************计算KM********************************************************************// double kmnew() { double upkm=0; double downkm=0; int i2,j2; for(i2=m+1;i2<N;i2++) { downkm=downkm+f0[i2]*f0[i2]+g0[i2-1]*g0[i2-1]; upkm=upkm+f0[i2]*g0[i2-1]; } downkm=downkm/2; upkm=-upkm; km=upkm/downkm; downkm=0; upkm=0; return(0); } //*****************************************************************************************************************// //**********************************************计算新的f0和g0*****************************************************// void update() { int n; double kmkm; double middle1; double middle2; int i3,j3; kmkm=km; for(n=1+m;n<N;n++) { middle1=f0[n]+kmkm*g0[n-1]; middle2=kmkm*f0[n]+g0[n-1]; f00[n]=middle1; g00[n]=middle2; } for(n=1+m;n<N;n++) { f0[n]=f00[n]; g0[n]=g00[n]; } middle1=0; middle2=0; } //*******************************************************************************************************************// //***************************************计算新的am******************************************************************// void update2() { double kmkmkm; double middle3; int n2; int i4,j4; kmkmkm=km; for(i4=1; i4<=m-1; i4++) { am[i4] = amchange[i4] + kmkmkm*amchange[m-i4]; //a[0][0]初值应为多少？？？ printf("a = %f , ",am[i4]); } am[m] = kmkmkm; for(i4=1; i4<=m; i4++) { amchange[i4] = am[i4]; } } //********************************************************************************************************************// 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; double xnew; unsigned long int i=0, j=0,xmxm; int i4,j4,n3,n4,controller; FILE * fpr, * fp_out; //a[0][1] = 1; //将二进制音频数据读入数组x[N]中 if((fpr = fopen("magic5.bin","rb")) == NULL) {printf("cannot open file\n"); exit(0); } for(i=0;i<N;i++) { fread(&x[i], 8, 1, fpr); //printf("%f\n",x[i]); } fclose(fpr); /*for(i=0;i<N;i++) { xout[i]=x[i]; printf("%f\n",xout[i]); } if((fp_out = fopen("boutdata.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],8, 1, fp_out); //fread(&x1[i][j],4, 1, fpr); //printf("%f\n",xout[j]); } fclose(fp_out); getchar(); */ controller=1; initial (); for(m=1;m<20;m++) { kmnew(); update2(); update(); //while(controller) //{ error=(1-km*km)*p0; error1=p0; p0=error; error1=error-error1; if(error1>-0.00005) { //controller=0; printf("success! m = %d\n",m); } else { /*m=m+1; kmnew(); update(); update2();*/ } // } //printf("m=%d\n",m); } /*for(i=20;i<N;i++) { for(j=1;j<=19;j++) { xmxm=xmxm+amtry[j]*x[i-j]; } xout[i]=-xmxm; xmxm=0; //printf("%f\n",x[i]); }*/ for(i=0;i<N;i++) { xout[i]=x[i]-f0[i]; //printf("%f\n",xout[i]); } if((fp_out = fopen("boutdata.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],8, 1, fp_out); //fread(&x1[i][j],4, 1, fpr); //printf("%f\n",xout[j]); } fclose(fp_out); getchar(); }

用于滤波，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();*/ }