| 查看: 1933 | 回复: 0 | ||
[资源]
【分享】最小平方反滤波（反褶积）【无重复】
|
|
//最小平方反滤波，学习反褶积必须掌握的第一个程序 #include #include #include /*解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 z=0.0; k=i-(n-1); for(j=0; j } } /* 线性卷积 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 y=0.0; for(j=0; j if(k>=0 && 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 } 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 //形成并解Toeplize方程 atlvs(r2,N,dirac,y); //反子波与子波卷积变回单位冲激序列，说明计算正确 conv(x,N,y,N,r1,M); for(i=0; i } |
回复此楼» 猜你喜欢
西南科技大学 调剂 绵阳
已经有10人回复
-
散金币，求基金中标！！
已经有47人回复
-
地质学论文润色/翻译怎么收费?
已经有81人回复
找到一些相关的精华帖子，希望有用哦~
非线性最小二乘拟合（初值问题）lsqcurvefit 函数的设定
已经有14人回复
- 求太阳能技术书籍+有效期至2010年5月1日
已经有11人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
