最小平方反滤波c程序小例

天涯

#include<stdio.h> #include<stdlib.h> #include<math.h> /*解Teoplize矩阵的莱文森递推算法 t[]--矩阵的第一行数据,n个 b[]--线性方程组最右边的一列数据 x[]--计算结果 */ int atlvs(double t[],int n,double b[],double x[]) { int i,j,k; double a,beta,q,c,h,*y,*s; s=(double *)calloc(n,sizeof(double)); y=(double *)calloc(n,sizeof(double)); a=t[0]; if (fabs(a)+1.0==1.0) { free(s); free(y); printf("fail\n"); return(-1);} y[0]=1.0; x[0]=b[0]/a; for (k=1; k<=n-1; k++) { beta=0.0; q=0.0; for (j=0; j<=k-1; j++) { beta=beta+y[j]*t[j+1]; q=q+x[j]*t[k-j]; } if (fabs(a)+1.0==1.0) { free(s); free(y); printf("fail\n"); return(-1);} c=-beta/a; s[0]=c*y[k-1]; y[k]=y[k-1]; if (k!=1) for (i=1; i<=k-1; i++) s=y[i-1]+c*y[k-i-1]; a=a+c*beta; if (fabs(a)+1.0==1.0) { free(s); free(y); printf("fail\n"); return(-1);} h=(b[k]-q)/a; for (i=0; i<=k-1; i++) { x=x+h*s; y=s;} x[k]=h*y[k]; } free(s); free(y); return(1); } //互相关 void corre(double x[],int m,double y[],int n,double z[],int l) { int i,j,k; for(i=0; i<l; i++) { z=0.0; k=i-(n-1); for(j=0; j<m; j++) if(j-k>=0&&j-k<n) z+=x[j]*y[j-k]; } } /* 线性卷积 y(n)=x(n)*h(n) m--length of x(n); n--length of h(n); l=m+n-1 length of y(n) */ void conv(double x[],int m,double h[],int n,double y[],int l) { int i,j,k; for(i=0; i<l; i++) { y=0.0; for(j=0; j<m; j++) { k=i-j; if(k>=0 && k<n) y += x[j]*h[k]; } } } //产生零相位子波,num--length,f0--frequency void WaveLet(double x[],int num,double f0) { double pi=3.1415926; double det=0.002;//2ms采样 double m=2.0;//余弦波衰减比例 double a; a=2*f0*f0*log(m); for(int i=0; i<num; i++) x=exp(-a*i*i*det*det)*cos(2*pi*f0*i*det); } void main() { const int N=129; const int M=2*N-1; double x[N];//子波 double y[N];//计算出的反子波 double dirac[N]={1.0};//单位冲激序列 double r1[M], r2[N]={0.0}; int i; WaveLet(x,N,40.0); //子波自相关 corre(x,N,x,N,r1,M); for(i=0; i<N; i++) r2=r1[M/2+i]; //形成并解Toeplize方程 atlvs(r2,N,dirac,y); //反子波与子波卷积变回单位冲激序列，说明计算正确 conv(x,N,y,N,r1,M); for(i=0; i<M; i++) printf("r1[%d]=%f\n",i,r1); }