第8题.cpp

最小平方反褶积程序,原代码

// 维纳滤波基本函数

  #include "stdlib.h"
  #include "math.h"
  #include "stdio.h"

/* 1. 莱文逊递推算法（t 双精度实型一维数组，长度n，存放Toeplitz阵中的自相关函数值。
                      n 整型变量。
                      b 数组，右端的常数向量，存放互相关函数值。
                      x 数组，方程组的解，即滤波因子。）
   本函数返回值＝－1，表示失败；＝1，成功。*/

int lvs(double t[], int n,double b[],double x[])
{
		int i,j,k;
		double a,beta,q,c,h,*y,*s;
	    s=(double *)malloc(n*sizeof(double));
		y=(double *)malloc(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[i]=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[i]=x[i]+h*s[i];
				 y[i]=s[i];
			 }
	         x[k]=h*y[k];
		}
	     free(s);free(y);
	     return(1);
} 

//褶积滤波 (x是原信号数组指针，n是原信号个数
          // h是滤波因子数组指针，m是滤波因子个数
         // y是输出信号数组指针，l是输出信号个数)

void con(double *x,int n,double *h,int m,double *y,int l)
{
	int i,j;
	for(i=0;i<l;i++)
	{
    	y[i]=0.0;
		for(j=0;j<n;j++)
		{
			if(i-j>=0&&i-j<m)
			y[i]+=x[j]*h[i-j];
		}
	}
}

//相关计算（计算x[n]与h[m]的互相关，结果为y[n])

void cor(double *x,int n,double *h,int m,double *y)
{
	int i,j;
	for(i=0;i<n;i++)
	{
		y[i]=0.0;
		for(j=0;j<m;j++)
		{
			if(i<=j)
				y[i]+=x[j]*h[j-i];
		}
	}
}
void main()
{
	int i,n=200,m=12,m1=60,l,l1;
	l=n+m-1;l1=l+m1-1;
	float t;
	double *h0,*h,*bx,*y,*y1,*x,*xx,b[12]={17.5,12.5,7.0,0.0,-3.5,-7.0,-3.0,-1.0,0.0,1.4,3.5,2.0};
	FILE *fp,*p,*p1,*p2,*p3,*p4;
	fp=fopen("D:\\胡学宝\\程序\\数据\\INPUT8.dat","r");
    p=fopen("D:\\胡学宝\\程序\\第8题\\第8题\\原反射系数(200).dat","w+");
	p1=fopen("D:\\胡学宝\\程序\\第8题\\第8题\\合成地震记录(211).dat","w+");
	p2=fopen("D:\\胡学宝\\程序\\第8题\\第8题\\自相关结果(211).dat","w+");
	p3=fopen("D:\\胡学宝\\程序\\第8题\\第8题\\反滤波因子(60).dat","w+");
	p4=fopen("D:\\胡学宝\\程序\\第8题\\第8题\\反射系数(222).dat","w+");
//动态申请存储空间
	h0=(double *)calloc(n,sizeof(double));
	h=(double *)calloc(m1,sizeof(double));
	bx=(double *)calloc(m1,sizeof(double));
	y=(double *)calloc(m1,sizeof(double));
	x=(double *)calloc(l,sizeof(double));
	xx=(double *)calloc(l,sizeof(double));
	y1=(double *)calloc(l1,sizeof(double));

//读数据
	
	for(i=0;i<n;i++)
	{
	   fscanf(fp,"%f",&t);
	   h0[i]=(double)t;
	}
	fclose(fp);
	for(i=0;i<m1;i++)
	{
		if(i==0) bx[i]=1.0;
		else bx[i]=0.0;
	}
//输出反射系数
	for(i=0;i<n;i++)
		fprintf(p,"%f\n",h0[i]);
	fclose(p);
//调用褶积函数，计算合成地震记录，并输出
	con(b,m,h0,n,x,l);
	for(i=0;i<l;i++)
		fprintf(p1,"%f\n",x[i]);
	fclose(p1);
//调用相关函数，计算自相关rxx[l]，并输出
	cor(x,l,x,l,xx);
	for(i=0;i<l;i++)
		fprintf(p2,"%f\n",xx[i]);
	fclose(p2);
//调用莱文森算法，计算反滤波因子，并输出
	i=lvs(xx,m1,bx,y);
	printf("%d\n",i);
	for(i=0;i<m1;i++)
		fprintf(p3,"%f\n",y[i]);
	fclose(p3);
//调用褶积函数，计算反射系数
	con(y,m1,x,l,y1,l1);
	for(i=0;i<l1;i++)
		fprintf(p4,"%f\n",y1[i]);
	fclose(p4);
}