24小时热门版块排行榜

返回列表

baobiao007

木虫 (职业作家)

应助: 201 (大学生)
金币: 6482.7
帖子: 3050
在线: 1009.9小时
虫号: 505962

[资源] 【分享】最小平方反滤波（反褶积）【无重复】

//最小平方反滤波，学习反褶积必须掌握的第一个程序

#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 if(j-k>=0&&j-k 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 {
y=0.0;
for(j=0; j { k=i-j;
if(k>=0 && k 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 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 r2=r1[M/2+i];
//形成并解Toeplize方程
      atlvs(r2,N,dirac,y);
//反子波与子波卷积变回单位冲激序列，说明计算正确
conv(x,N,y,N,r1,M);
for(i=0; i       printf("r1[%d]=%f\n",i,r1);
}