地震子波波形显示及一维地震合成记录制作

//最小相位子波
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#define PI 3.1415962

void waveletfun(int nt_wavelet,float dt,float *wavelet,float fm,float m1m2);
void reflectfun(int nt_record,float *Vp,float *des,float *reflect);
void convolutionfun(int nt_record, int nt_wavelet, float *record, float *reflect, float *wavelet);

// 主程序部分
void main()
{
	// 变量的定义
	int nt_wavelet=150;
	float *wavelet;
	float dt=0.002;
	int it;
	int j;
	float fm=10;
	float m1m2=2.0;
	int nt_record = 600;
	float *reflect;
	float *des;
	float *Vp;
	float *record=0;

	// 定义文件
	FILE *waveletfp;
	waveletfp = fopen ("minphase_wavelet.txt","w");
	FILE *reflectfp;
	reflectfp = fopen ("minphase_reflect.txt","w");
	FILE *recordfp;
	recordfp = fopen ("minphase_record.txt","w");
	
	//分配内存
	wavelet = (float *)malloc(sizeof(float)*nt_wavelet);
	reflect = (float *)malloc(sizeof(float)* nt_record);
	Vp = (float *)malloc(sizeof(float)* nt_record);
	des = (float *)malloc(sizeof(float)* nt_record);
	record = (float *)malloc(sizeof(float)* (nt_wavelet+nt_record-1));

	// 子波求取
	waveletfun(nt_wavelet, dt, wavelet, fm, m1m2);
	for(it=0; it<nt_wavelet; it++)
	{		
		fprintf( waveletfp,"%f\n",wavelet[it] );
	
	}

	// 反射率计算
	reflectfun( nt_record, Vp, des, reflect );
	for(it=0; it<nt_record; it++)
	{
		fprintf( reflectfp,"%f\n",reflect[it] );
	}
	
	// 褶积形成地震记录
	convolutionfun( nt_record, nt_wavelet, record, reflect, wavelet);
	for(j=0; j<nt_record+nt_wavelet-1; j++)
	{
		fprintf( recordfp, "%f\n", record[j]);
	}
}

//子波求取函数
void waveletfun(int nt_wavelet,float dt,float *wavelet,float fm,float m1m2)
{
	int it;
	float t;
	for(it=0; it<nt_wavelet; it++)
	{
		t=it*dt;
		wavelet[it] = exp(-2*fm*fm*t*t*log(m1m2))*sin(2*PI*fm*t);
	}   
}

//反射率计算函数
void reflectfun(int nt_record,float *Vp,float *des,float *reflect)
{
	// 速度密度参数赋值
	int it;
	for (it=0; it<nt_record; it++)
	{
		if(it<nt_record/5)
		{
			Vp[it]=3000;
			des[it]=2000;
		}
		else if(it>4*nt_record/5)
		{
			Vp[it]=3333;
			des[it]=2222;
		}
		else if(it>=nt_record/5 && it<2*nt_record/5)
		{
			Vp[it]=3300;
			des[it]=2200;
		}
		else if( it>=nt_record/3 &&  it<2*nt_record/3)
		{
			Vp[it]=3220;
			des[it]=2050;
		}
		else
		{
			Vp[it]=3120;
			des[it]=2020;
		}
	}
	for (it=0; it<nt_record; it++)
	{
	    reflect[it]=(Vp[it+1]*des[it+1]-Vp[it]*des[it])/(Vp[it+1]*des[it+1]+Vp[it]*des[it]);	
	}
	reflect[nt_record-1]=0.0;
}

//褶积函数
void convolutionfun(int nt_record, int nt_wavelet, float *record, float *reflect, float *wavelet)
{
	int it;
	int j;
	int a;
	int len=nt_record+nt_wavelet-1;
	float *p=new float[len];
	for(it=0; it<nt_wavelet; it++)
	{
		p[it]=wavelet[it]; 
	}
	for(it=nt_wavelet; it<len; it++)
	{
		p[it]=0.0;
	}
	for(j=0; j<len; j++)
	{
		record[j]=0.0;	
		if(j<nt_record)
		{
			a=j;
		}
		else
		{
			a=nt_record-1;
		}
		for(it=0; it<=a; it++)
		{
			record[j]+=reflect[it]*p[j-it];
		}
		
	}
}