convolution model.cpp

石油地震勘探褶积的一个简单程序显示

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define PI 3.1415926
#define f0 40

void main()
{
	int   i,j,k,num1,num2;   
    double   dt,*R,wavelet[81],reflect[460],*s,*x; 
	FILE *fp,*reflection,*record;
	/***************************Ricker子波离散化******************************/ 
	dt=0.001;                                                  /*离散时间间隔*/
	R=wavelet;
	k=-40;
	for(i=0;k<=40;k++)                                    
	{
		wavelet[i]=(1-2*(PI*f0*k*dt)*(PI*f0*k*dt))*exp(-(PI*f0*k*dt)*(PI*f0*k*dt));
		i++;
	}
	/****************************输出数据到文件*******************************/
	if((fp=fopen("D:\\code\\data\\X_wavelet.txt","w+"))==NULL)
	{
		printf("cann't open file\n");
		return;
	}
	for(i=0;i<=80;i++)
		if (fprintf(fp,"%f\n",wavelet[i])==0)
			printf("file write error\n");		
	fclose(fp);
	/****************************输入反射系数序列******************************/
	num1=81;
	num2=231;
	for(i=0;i<=230;i++)
		reflect[i]=0;
	reflect[40]=0.2353;
	reflect[75]=0.1871;
	reflect[110]=-0.3731;
	reflect[140]=0.4247;
	reflect[175]=-0.2455;
	reflect[230]=0.1871;
	s=(double   *)malloc((num1+num2-1)*sizeof(double   ));
	x=(double   *)malloc(num2*sizeof(double   ));
	/**************************输出反射系数离散序列到文件**********************/
	if((reflection=fopen("D:\\code\\data\\reflection.txt","w+"))==NULL)
	{
		printf("cann't open file\n");
		return;
	}
	for(i=0;i<=230;i++)
		if (fprintf(reflection,"%f\n",reflect[i])==0)
			printf("file write error\n");		
	fclose(reflection);
	/*********************wavelet[]与reflect[]褶积的过程***********************/	
	 if(num1>=num2)                   /*先考虑wavelet数组长度大于reflect的情况*/              
	{
		 for(i=0;i<num1+num2-1;i++)
		 {
			 if(i<num2)
			{
				 for(s[i]=0,j=0;j<=i;j++)
					 s[i]+=wavelet[j]*reflect[i-j];
			 }
			 if(num2<=i&&i<num1)
			 {
				 for(s[i]=0,j=i-num2+1;j<=i;j++)
					 s[i]+=wavelet[j]*reflect[i-j];
			 }
			if(i>=num1)
			{
				for(s[i]=0,j=i-num2+1;j<num1;j++)
					s[i]+=wavelet[j]*reflect[i-j];
			}
		 }
	 }
	else                           /*再考虑reflect数组长度大于wavelet的情况*/  
	{
		for(i=0;i<num1+num2-1;i++)
		{
			if(i<num1)
			{
				for(s[i]=0,j=0;j<=i;j++)
					s[i]+=wavelet[j]*reflect[i-j];
			}
			if(num1<=i&&i<num2-1)
			{
				for(s[i]=0,j=0;j<num1;j++)
					s[i]+=wavelet[j]*reflect[i-j];
			}
			if(num2-1<=i)
			{
				for(s[i]=0,j=i-num2+1;j<num1;j++)
					s[i]+=wavelet[j]*reflect[i-j];
			}
		}
	}

    for(j=0,i=40;i<(num1+num2-1-40);i++,j++)                                /*输出褶积结果*/  
	    x[j]=s[i];
    for(j=0;j<num2;j++)                                /*输出褶积结果*/  
	    printf("x[%d]=%f\n ",j,x[j]);
    /*******************************输出结果离散序列***************************/
    if((record=fopen("D:\\code\\data\\record.txt","w+"))==NULL)
	{
		printf("cann't open file\n");
		return;
	}
	for(i=0;i<num2;i++)
		if (fprintf(record,"%f\n",x[i])==0)
			printf("file write error\n");		
	fclose(record);
}