任务一:查找地震子波概念并调查“零相位雷克子波”表达式。
地震子波是一段具有确定的起始时间、能量有限且有一定延续长度的信号,它是地震记录中的基本单元。
零相位雷克子波表达式:
任务二:设主频为30Hz,基于c语言实现“雷克子波”程序。
任务三:给定随机反射系数,并实现其c程序。
任务四:基于“雷克子波”以及随即反射系数,合成地震记录,实现其c程序。
任务五:修改零相位”雷克子波“表达式,得到90°、180°、270°相位的雷克子波表达式,并实现其c程序。
仔细分析将会发现课设的难度并不高,相对于已经上过c程序设计课程的同学来说,只要沉下心来,认真思考、并查阅一些文献资料,两天左右的时间完全可以调试出来。
下面是笔者自己编写的c程序,可能有所疏漏 欢迎大家批评指正。
#include<stdio.h>
#include<iostream>
#include<time.h>
#include<math.h>
#include<stdlib.h>
using namespace std;
#define PI 3.1415926
#define T 0.001
#define FM 30
void convolution(double a[],double b[],double c[],int M, int N);
int main(void)
{
double Ricker[200];
FILE*fpLei,*fpFan,*fpDi;
fpLei=fopen("雷克子波.cvs","w");
if(!fpLei)
{
cout<<"雷克子波文件打开失败!!!"<<endl;
exit(0);
}
fprintf(fpLei,"%s,%s\n","时间t","雷克子波");
for(int i=0;i<200;++i)
{
if(i==0)
cout<<"雷克子波文件数据开始导入!!!"<<endl;
Ricker[i]=(1-2*PI*(__int64)(i-50)*FM*T*PI*(i-50)*FM*T)*exp(-PI*(i-50)*FM*T*PI*(i-50)*FM*T);
fprintf(fpLei,"%.4f,%.4f\n",T*(__int64)(i-50),Ricker[i]);
if(i==199)
cout<<"雷克子波文件数据全部导入成功!!!"<<endl;
}
fclose(fpLei);
srand(unsigned(time(NULL)));
double Reflect[500];
fpFan=fopen("反射系数.cvs","w");
if(!fpFan)
{
cout<<"反射系数文件打开失败!!!"<<endl;
exit(0);
}
for(int i=0;i<500;i++)
{
if(i==0)
{
fprintf(fpFan,"%s","反射系数");
cout<<"反射系数文件开始读入!!!"<<endl;
}
int reflect=rand()%200+1;
if(reflect<=100)
Reflect[i]=reflect/100.0;
else
Reflect[i]=(reflect-200)/100.0;
fprintf(fpFan,"%lf\n",Reflect[i]);
if(i==499)
cout<<"反射系数数据文件读入完成!!!"<<endl;
}
fclose(fpFan);
double Di[599];
convolution(Ricker,Reflect,Di,100,500);
fpDi=fopen("地震子波.cvs","w+");
for(int i=0;i<599;++i)
{
if(i==0)
cout<<"地震子波文件开始读入!!!"<<endl;
fprintf(fpDi,"%lf\n",Di[i]) ;
if(i==598)
cout<<"地震子波文件读入完成!!!"<<endl;
}
fclose(fpDi);
double Ricker_90[199];
double Zheji_90[100];
for(int i=0;i<100;++i)
{
Zheji_90[i]=(i-50)/PI;
}
convolution(Ricker,Zheji_90,Ricker_90,100,100);
FILE*Lei_90,*Lei_180,*Lei_270;
Lei_90=fopen("雷克子波90.cvs","W");
if(!Lei_90)
{
cout<<"雷克子波90文件的读入失败!!!"<<endl;
exit(0);
}
for(int i=0;i<100;i++)
{
if(0==i)
{
cout<<"开始读入雷克子波90文件数据!!!"<<endl;
fprintf(Lei_90,"%s,%s\n","时间t","l雷克子波");
}
fprintf(fpLei,"%.4f,%.4f\n",T*(__int64)(i-50),Ricker_90[i]);
if(i==99)
cout<<"雷克子波数据全部读入成功!!!"<<endl;
}
fclose(Lei_90);
double Ricker_180[199];
double Zheji_180[100];
for (int i=0;i<100;i++)
{
Zheji_180[i]=2*(i-50)/PI;
}
convolution(Ricker,Zheji_180,Ricker_180,100,100);
Lei_180=fopen("雷克子波180.cvs","w");
if(!Lei_180)
{
cout<<"雷克子波180文件读入失败!!!"<<endl;
exit(0);
}
for(int i=0;i<100;++i)
{
if(0==i)
{
cout<<"开始读入雷克子波180文件!!!"<<endl;
fprintf(Lei_180,"%s,%s\n","时间t","l雷克子波");
}
fprintf(fpLei,"%.4f,%.4f\n",T*(__int64)(i-50),Ricker_180[i]);
if(i==99)
cout<<"雷克子波180数据读入完成!!!"<<endl;
}
fclose(Lei_180);
double Ricker_270[199];
double Zheji_270[100];
for(int i=0;i<100;++i)
{
Zheji_270[i]=3*(i-50)/PI;
}
convolution(Ricker,Zheji_270,Ricker_270,100,100);
Lei_270=fopen("雷克子波270.cvs","w");
if(!Lei_270)
{
cout<<"雷克子波270数据读入失败!!!"<<endl;
exit(0);
}
for(int i=0;i<100;++i)
{
if(0==i)
{
cout<<"开始读入雷克子波270文件数据!!!"<<endl;
fprintf(Lei_270,"%S,%S\n","时间t","l雷克子波");
}
fprintf(fpLei,"%.4f,%.4f\n",T*(__int64)(i-50),Ricker_270[i]);
if(i==99)
cout<<"雷克子波270文件数据读入完成!!!"<<endl;
}
fclose(Lei_270);
system("pause");
return 0;
}
/*褶积计算*/
void convolution(double a[],double b[],double c[],int M,int N)
{
int L=M+N-1;
for(int i=0;i<L;i++)
{
double tp=0.0;
for(int j=0;j<M;j++)
{
if((i-j)>=0&&(i-j)<N)
tp+=a[j]*b[i-j];
}
c[i]=tp;
tp=0.0;
}
}
在程序运行之后,不太熟悉的同学将会发现运行后的数据记录找不到了,要先试着读懂程序之后将数据保存在指定路径。
数据保存之后打开会是.cvs文件,这种类型的文件可以用excel打开,之后框选所得数据,生成可视化表格便可以完成本次课程设计,最后便是撰写课程报告了。
以上程序运行结果仅供参考。