地震资料处理实习报告Word文件下载.docx
- 文档编号:8205754
- 上传时间:2023-05-10
- 格式:DOCX
- 页数:43
- 大小:298.92KB
地震资料处理实习报告Word文件下载.docx
《地震资料处理实习报告Word文件下载.docx》由会员分享,可在线阅读,更多相关《地震资料处理实习报告Word文件下载.docx(43页珍藏版)》请在冰点文库上搜索。
高通滤波因子
频率域高通滤波因子
5.频谱分析
正弦函数频谱
尖脉冲频谱
地震波频谱
6.分析补零对振幅谱的影响
正弦函数n=60
正弦函数n=128
\
非周期波形n=60
非周期波形n=64
非周期波形n=128
7.线性褶积与圆周褶积
线性褶积模型
圆周褶积结果
长度与圆周褶积相等的线性褶积
8.最小平方反滤波
原始反射系数序列
求出的反射系数序列
9.零相位转换
非零相位子波
零相位子波
11、地震记录
选择第一道为参考到静校正结果
二、程序
1.褶积滤波
#include"
stdio.h"
math.h"
conv.c"
#definepi3.1415926
#defineN100
#definedt0.002main()
{
floatx[100],h[101],h1[101],y_low[200],y_band[200];
floatdf;
inti,m=100,n=101,l=m+n-1;
floatf=70.0;
floatf1=10.0;
floatf2=80.0;
FILE*fp1,*fp2,*fp3,*fp4,*fp5;
fp1=fopen("
INPUT1.DAT"
"
r"
);
for(i=0;
i<
=N;
i++)
fscanf(fp1,"
%f"
&
x[i]);
}
fp4=fopen("
h_low.dat"
w"
//低通滤波设计for(i=-50;
=50;
if(i==0)
h[i+50]=2*pi*f/pi;
else
h[i+50]=sin(2*pi*f*i*dt)/(pi*i*dt);
fprintf(fp4,"
%f"
h[i+50]);
//outputlowpassfilter
fp2=fopen("
synthesisdata_lowpass.DAT"
conv(x,m,h,n,y_low,l);
for(i=50;
l-50;
fprintf(fp2,"
%f\n"
y_low[i]);
}
//带通滤波器fp5=fopen("
h_band.dat"
for(i=-50;
h1[i+50]=140;
h1[i+50]=sin(2*pi*f2*i*dt)/(pi*i*dt)-sin(2*pi*f1*i*dt)/(pi*i*dt);
fprintf(fp5,"
h1[i+50]);
//outputbandpassfilter
fclose(fp5);
conv(x,m,h1,n,y_band,l);
fp3=fopen("
synthesisdata_bandpass.DAT"
for(i=50;
fprintf(fp3,"
y_band[i]);
fclose(fp1);
fclose(fp2);
fclose(fp3);
2.快变滤波
stdlib.h"
fft.c"
main()
double*xr,*xi;
float*H;
inti,np,nfft,k;
floatt,dt,df,f,z,fc1,fc2;
FILE*fpar,*fp1,*fp2,*fp3;
//从参数文件中获得截至频率
fpar=fopen("
lowpassfilter.par"
fscanf(fpar,"
%f%f"
fc1,&
fc2);
np=100;
k=log(np)/log
(2);
if(np>
pow(2,k))k=k+1;
nfft=pow(2,k);
dt=0.002;
df=1.0/(nfft*dt);
xr=(double*)calloc(nfft,sizeof(double));
xi=(double*)calloc(nfft,sizeof(double));
H=(float*)calloc(nfft,sizeof(float));
//readx(n)fp1=fopen("
100;
z);
xr[i]=z;
//补零至128位
for(i=100;
nfft;
xr[i]=0.0;
xi[i]=0.0;
//transfertofrequencydoaminfft(xr,xi,k,1);
//generatelowpassfilter(zero-phase)for(i=0;
nfft/2;
f=i*df;
if(f<
=fc2&
&
f>
=fc1)H[i]=1.0;
elseH[i]=0.0;
//滤波器对称for(i=nfft/2;
H[i]=H[nfft-i];
//filteringinfrequencydomainfor(i=0;
xr[i]=xr[i]*H[i];
xi[i]=xi[i]*H[i];
fft(xr,xi,k,-1);
fp2=fopen("
lowpass2.dat"
for(i=0;
xr[i]);
//获取高通截至频率fpar=fopen("
bandpass.par"
=nfft/2;
else
H[i]=0.0;
for(i=nfft/2+1;
fp3=fopen("
bandpass2.dat"
3.褶积滤波与递归滤波
褶积滤波
#definePI3.1415926main()
voidconv();
floatx[50],h[20],y[69],hreverse[20],hzero[39],yreverse[69];
floatdt=0.002;
inti,m,n,l,p,q;
FILE*fp1,*fp2,*fp3,*fp4,*fp5,*fp6;
m=50;
n=20;
l=m+n-1;
//readx(n)fp1=fopen("
INPUT3.DAT"
50;
//readfilterfactorh(n)fp2=fopen("
hn.dat"
fp5=fopen("
h_reverse.dat"
20;
fscanf(fp2,"
h[i]);
hreverse[i]=h[19-i];
fprintf(fp5,"
hreverse[i]);
conv(x,m,h,n,y,l);
//非零相位褶积滤波fp3=fopen("
con_filter.dat"
l;
y[i]);
p=n+n-1;
q=m+p-1;
//构造零相位滤波因子
conv(h,n,hreverse,n,hzero,p);
fp6=fopen("
zerophasefilterfactor.dat"
p;
fprintf(fp6,"
hzero[i]);
fclose(fp6);
//零相位滤波conv(x,m,hzero,p,yreverse,q);
fp4=fopen("
convfilterreverse.dat"
q;
yreverse[i]);
fclose(fp4);
递归滤波
#include<
stdio.h>
math.h>
stdlib.h>
#definenp50voidmain()
float*x,*a,*b,*fil1,*fil2;
inti;
voidrecur1();
voidrecur2();
x=(float*)malloc(np*sizeof(float));
a=(float*)malloc(np*sizeof(float));
b=(float*)malloc(np*sizeof(float));
fil1=(float*)malloc(np*sizeof(float));
fil2=(float*)malloc(np*sizeof(float));
//输入地震数据
fp1=fopen("
np;
i++)fscanf(fp1,"
x+i);
fclose(fp1);
//输入a数组值fp2=fopen("
a(n).txt"
5;
i++)fscanf(fp2,"
a+i);
for(i=5;
i++)a[i]=0.0;
//输入b数组值fp3=fopen("
b(n).txt"
i++)fscanf(fp3,"
b+i);
fclose(fp3);
i++)b[i]=0.0;
//正向递归滤波
recur1(x,a,b,fil1);
正向递归结果.DAT"
wb"
%12.4f\n"
fil1[i]);
printf("
\n"
//反向递归滤波
recur2(fil1,a,b,fil2);
反向递归结果.DAT"
fil2[i]);
voidrecur1(floatx[],floata[],floatb[],floatfil1[])
inti,j,k;
floaty1[np],y2[np];
{y1[i]=0.0;
y2[i]=0.0;
for(j=0;
j<
=i;
j++)y1[i]=y1[i]+a[j]*x[i-j];
if(i==0)y2[i]=0.0;
else{for(k=1;
k<
k++)y2[i]=y2[i]+b[k]*fil1[i-k];
fil1[i]=y1[i]-y2[i];
voidrecur2(floatfil1[],floata[],floatb[],floatfil2[])
floaty3[np],y4[np];
for(i=np-1;
i>
=0;
i--)
else{
y3[i]=0.0;
y4[i]=0.0;
=np-1-i;
j++)y3[i]=y3[i]+a[j]*fil1[i+j];
if(i==np-1)y4[i]=0.0;
for(k=1;
k++)y4[i]=y4[i]+b[k]*fil2[i+k];
fil2[i]=y3[i]-y4[i];
4.设计高通滤波因子
#defineN101
#definedt0.004
#definePI3.1415926
double*h,*hi;
inti,k,nfft;
FILE*fp1,*fp2;
k=log(N)/log
(2);
if(N>
h=(double*)malloc(nfft*sizeof(double));
hi=(double*)malloc(nfft*sizeof(double));
h[i+50]=1.0/dt-60;
h[i+50]=-sin(2*PI*30.0*i*dt)/(PI*i*dt);
128;
h[i]=0.0;
hi[i]=0.0;
timedomain.dat"
fprintf(fp1,"
h[i]);
fft(h,hi,k,1);
frequencydoamin.dat"
sqrt(h[i]*h[i]+hi[i]*hi[i]));
printf("
itisok!
5.分析补零对振幅谱的影响
1、不补零
dft.c"
#defineN60
#definedt0.004main()
floatx[N],xr[N],xi[N],w[N],wr[N],wi[N],z;
inti,k;
FILE*fp,*fp1,*fp2,*fp3;
sin.dat"
N;
x[i]=sin(2.0*PI*(i+1)/30);
x[i]);
fp=fopen("
WAVE.DAT"
fscanf(fp,"
w[i]=z;
fclose(fp);
itisok!
dft(x,xr,xi,N);
dft(w,wr,wi,N);
frequencydomain1_60.dat"
frequencydomain2_60.dat"
wr[i]);
2、补到64位
voidfft();
double*x1,*xi1,*x2,*xi2;
floatz;
inti,k,nfft;
FILE*fp,*fp1,*fp2,*fp3,*fp4;
k=log(N)/log
(2);
if(N>
x1=(double*)malloc(nfft*sizeof(double));
xi1=(double*)malloc(nfft*sizeof(double));
x2=(double*)malloc(nfft*sizeof(double))
xi2=(double*)malloc(nfft*sizeof(double));
x1[i]=sin(2*PI*(i+1)/30);
x2[i]=z;
xi1[i]=0.0;
xi2[i]=0.0;
//补到64位
for(i=N;
64;
x1[i]=0.0;
x2[i]=0.0;
fft(x1,xi1,6,1);
fft(x2,xi2,6,1);
frequencydomain1_64.dat"
frequencydomain2_64.dat"
sqrt(x1[i]*x1[i]+xi1[i]*xi1[i])*dt);
sqrt(x2[i]*x2[i]+xi2[i]*xi2[i])*dt);
3、补到128位
inti,nfft;
FILE*fp,*fp1,*fp2;
nfft=128;
x2=(double*)malloc(nfft*sizeof(double));
xi2=(double*)malloc(nfft*sizeof(double));
//补到128位
fft(x1,xi1,7,1);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 地震 资料 处理 实习 报告