最小平方反滤波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;
if (fabs(a)+1.0==1.0)
{ free(s); free(y); printf("fail\n"); return(-1);}
y=1.0; x=b/a;
for (k=1; k<=n-1; k++)
{ beta=0.0; q=0.0;
for (j=0; j<=k-1; j++)
{ beta=beta+y*t;
q=q+x*t;
}
if (fabs(a)+1.0==1.0)
{ free(s); free(y); printf("fail\n"); return(-1);}
c=-beta/a; s=c*y; y=y;
if (k!=1)
for (i=1; i<=k-1; i++)
s=y+c*y;
a=a+c*beta;
if (fabs(a)+1.0==1.0)
{ free(s); free(y); printf("fail\n"); return(-1);}
h=(b-q)/a;
for (i=0; i<=k-1; i++)
{ x=x+h*s; y=s;}
x=h*y;
}
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*y;
}
}
/* 线性卷积
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*h;
}
}
}
//产生零相位子波,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;//子波
double y;//计算出的反子波
double dirac={1.0};//单位冲激序列
double r1, r2={0.0};
int i;
WaveLet(x,N,40.0);
//子波自相关
corre(x,N,x,N,r1,M);
for(i=0; i<N; i++)
r2=r1;
//形成并解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);
}
页:
[1]