1、地震资料处理实习报告 地震资料数据处理课程设计总结报告专业班级： 地球物理学1302班 姓 名： 学 号： 成 绩： 2016年12月31日一 、 设计内容（1）褶积滤波（2）快变滤波（3）褶积滤波与快变滤波的比较（4）设计高通滤波因子（5）频谱分析（6）分析补零对振幅谱的影响（7）线性褶积与循环褶积（8）最小平方反滤波（9）零相位转换（11）静校正二、程序一、设计内容 1、褶积滤波 理想低通滤波因子 理想带通滤波因子原始信号低通滤波结果带通滤波结果2、快变滤波原始数据快变低通滤波结果快变带通滤波结果3褶积滤波与递归滤波的比较原始数据零相位褶积滤波结果非零相位褶积滤波结果零相位递归滤波非零相位
2、递归滤波4设计高通滤波因子原始数据高通滤波因子频率域高通滤波因子5频谱分析正弦函数频谱尖脉冲频谱地震波频谱6分析补零对振幅谱的影响正弦函数n=60正弦函数n=128非周期波形n=60非周期波形n=64非周期波形n=1287线性褶积与圆周褶积线性褶积模型圆周褶积结果长度与圆周褶积相等的线性褶积8最小平方反滤波原始反射系数序列求出的反射系数序列9零相位转换非零相位子波零相位子波11、地震记录原始数据选择第一道为参考到静校正结果二、程序1. 褶积滤波# include stdio.h# include math.h# include conv.c# define pi 3.1415926# defi
3、ne N 100# define dt 0.002 main()float x100, h101, h1101, y_low200, y_band200; float df;int i,m=100,n=101,l=m+n-1; float f=70.0;float f1=10.0; float f2=80.0;FILE *fp1,*fp2,*fp3,*fp4,*fp5; fp1=fopen(INPUT1.DAT,r);for(i=0;i=N;i+)fscanf(fp1,%f,&xi);fp4=fopen(h_low.dat,w);/低通滤波设计 for(i=-50;i=50;i+)if(i=0
4、)hi+50=2*pi*f/pi;elsehi+50=sin(2*pi*f*i*dt)/(pi*i*dt);fprintf(fp4,%f ,hi+50); /output lowpass filterfp2=fopen(synthesisdata_lowpass.DAT,w); conv(x,m,h,n,y_low,l);for(i=50;il-50;i+)fprintf(fp2,%fn,y_lowi); /带通滤波器 fp5=fopen(h_band.dat,w); for(i=-50;i=50;i+)if(i=0)h1i+50=140;elseh1i+50=sin(2*pi*f2*i*dt
5、)/(pi*i*dt)-sin(2*pi*f1*i*dt)/(pi*i*dt);fprintf(fp5,%fn,h1i+50); / output bandpass filterfclose(fp5); conv(x,m,h1,n,y_band,l);fp3=fopen(synthesisdata_bandpass.DAT,w); for(i=50;ipow(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,si
6、zeof(double); H=(float*)calloc(nfft,sizeof(float);/ read x(n) fp1=fopen(INPUT1.DAT,r);for(i=0;i100;i+)fscanf(fp1,%f,&z); xri=z;fclose(fp1);/补零至 128 位for(i=100;infft;i+)xri=0.0;for(i=0;infft;i+)xii=0.0;/transfer to frequency doamin fft(xr,xi,k,1);/generate lowpass filter(zero-phase) for(i=0;infft/2;i
7、+)f=i*df;if(f=fc1)Hi=1.0; else Hi=0.0;/滤波器对称 for(i=nfft/2;infft;i+)f=i*df; Hi=Hnfft-i;/filtering in frequency domain for(i=0;infft;i+)xri=xri*Hi;xii=xii*Hi;fft(xr,xi,k,-1); fp2=fopen(lowpass2.dat,w); for(i=0;infft;i+)fprintf(fp2,%fn,xri);fclose(fp2);/获取高通截至频率 fpar=fopen(bandpass.par,r); fscanf(fpar,
8、%f%f,&fc1,&fc2); fp1=fopen(INPUT1.DAT,r);for(i=0;i100;i+)fscanf(fp1,%f,&z); xri=z;for(i=100;infft;i+)xri=0.0;for(i=0;infft;i+)xii=0.0;/transfer to frequency doamin fft(xr,xi,k,1);/generate lowpass filter(zero-phase) for(i=0;i=nfft/2;i+)f=i*df;if(f=fc1)Hi=1.0; elseHi=0.0;for(i=nfft/2+1;infft;i+)f=i*d
9、f; Hi=Hnfft-i;/filtering in frequency domain for(i=0;infft;i+)xri=xri*Hi;xii=xii*Hi;fft(xr,xi,k,-1); fp3=fopen(bandpass2.dat,w); for(i=0;infft;i+)fprintf(fp3,%fn,xri);3. 褶积滤波与递归滤波褶积滤波# include stdio.h# include math.h# include stdlib.h# include conv.c# include fft.c# define PI 3.1415926 main()void co
10、nv();float x50,h20,y69,hreverse20,hzero39,yreverse69; float dt=0.002;int i,m,n,l,p,q;FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6; m=50;n=20;l=m+n-1;/read x(n) fp1=fopen(INPUT3.DAT,r);for(i=0;i50;i+)fscanf(fp1,%f,&xi);fclose(fp1);/read filterfactor h(n) fp2=fopen(hn.dat,r); fp5=fopen(h_reverse.dat,w); for(i=0
11、;i20;i+)fscanf(fp2,%f,&hi); hreversei=h19-i; fprintf(fp5,%fn,hreversei);fclose(fp2);fclose(fp5);conv(x,m,h,n,y,l);/非零相位褶积滤波 fp3=fopen(con_filter.dat,w); for(i=0;il;i+)fprintf(fp3,%fn,yi);fclose(fp3);p=n+n-1; q=m+p-1;/构造零相位滤波因子conv(h,n,hreverse,n,hzero,p); fp6=fopen(zerophasefilterfactor.dat,w); for(
12、i=0;ip;i+)fprintf(fp6,%fn,hzeroi);fclose(fp6);/零相位滤波 conv(x,m,hzero,p,yreverse,q); fp4=fopen(convfilterreverse.dat,w); for(i=0;iq;i+)fprintf(fp4,%fn,yreversei);fclose(fp4);递归滤波#include#include#include#define np 50 void main()float *x,*a,*b,*fil1,*fil2; int i;void recur1(); void recur2();FILE *fp1,*f
13、p2,*fp3,*fp4,*fp5; 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(INPUT3.DAT,r);for(i=0;inp;i+)fscanf(fp1,%f,x+i); fclose(fp1);/输入 a 数组值 fp2=fopen(a(n
14、).txt,r); for(i=0;i5;i+)fscanf(fp2,%f,a+i);fclose(fp2); for(i=5;inp;i+)ai=0.0;/输入 b 数组值 fp3=fopen(b(n).txt,r); for(i=0;i5;i+)fscanf(fp3,%f,b+i); fclose(fp3); for(i=5;inp;i+)bi=0.0;/正向递归滤波recur1(x,a,b,fil1); fp4=fopen(正向递归结果.DAT,wb); for(i=0;inp;i+)fprintf(fp4,%12.4fn,fil1i);fclose(fp4); for(i=0;inp;
15、i+)printf(%12.4fn,fil1i);printf(n);/反向递归滤波recur2(fil1,a,b,fil2); fp5=fopen(反向递归结果.DAT,wb); for(i=0;inp;i+)fprintf(fp5,%12.4fn,fil2i);fclose(fp5); for(i=0;inp;i+)printf(%12.4fn,fil2i);void recur1(float x,float a,float b,float fil1)int i,j,k;float y1np,y2np; for(i=0;inp;i+) y1i=0.0;y2i=0.0;for(j=0;j=i
16、;j+)y1i=y1i+aj*xi-j;if(i=0)y2i=0.0;else for(k=1;k=0;i-)elsey3i=0.0;y4i=0.0;for(j=0;j=np-1-i;j+)y3i=y3i+aj*fil1i+j; if(i=np-1)y4i=0.0;for(k=1;kpow(2,k)k=k+1; nfft=pow(2,k);h=(double*)malloc(nfft*sizeof(double); hi=(double*)malloc(nfft*sizeof(double); for(i=-50;i=50;i+)if(i=0)hi+50=1.0/dt-60;elsehi+50
17、=-sin(2*PI*30.0*i*dt)/(PI*i*dt);for(i=100;i128;i+)hi=0.0;for(i=0;i128;i+)hii=0.0;fp1=fopen(timedomain.dat,w); for(i=0;i128;i+)fprintf(fp1,%fn,hi);fclose(fp1); fft(h,hi,k,1);fp2=fopen(frequencydoamin.dat,w); for(i=0;i128;i+)fprintf(fp2,%fn,sqrt(hi*hi+hii*hii);fclose(fp2); printf(it is ok!n);5. 分析补零对振
18、幅谱的影响1、 不补零# include stdio.h# include math.h# include dft.c# define N 60# define PI 3.1415926# define dt 0.004 main()float xN,xrN,xiN,wN,wrN,wiN,z; int i,k;FILE *fp,*fp1,*fp2,*fp3;fp3=fopen(sin.dat,w); for(i=0;iN;i+)xi=sin(2.0*PI*(i+1)/30);fprintf(fp3,%fn,xi);fclose(fp3);fp=fopen(WAVE.DAT,r); for(i=
19、0;iN;i+)fscanf(fp,%f,&z); wi=z;fclose(fp);printf(it is ok !n); dft(x,xr,xi,N);dft(w,wr,wi,N); fp1=fopen(frequencydomain1_60.dat,w); fp2=fopen(frequencydomain2_60.dat,w); for(i=0;ipow(2,k)k=k+1; nfft=pow(2,k);x1=(double*)malloc(nfft*sizeof(double); xi1=(double*)malloc(nfft*sizeof(double); x2=(double*
20、)malloc(nfft*sizeof(double)xi2=(double*)malloc(nfft*sizeof(double); for(i=0;iN;i+)x1i=sin(2*PI*(i+1)/30);fp=fopen(WAVE.DAT,r); for(i=0;iN;i+)fscanf(fp,%f,&z); x2i=z;for(i=0;iN;i+)xi1i=0.0;xi2i=0.0;/补到 64 位for(i=N;i64;i+)x1i=0.0;x2i=0.0;xi1i=0.0;xi2i=0.0;fft(x1,xi1,6,1);fft(x2,xi2,6,1); fp1=fopen(fre
21、quencydomain1_64.dat,w); fp2=fopen(frequencydomain2_64.dat,w); for(i=0;infft;i+)fprintf(fp1,%fn,sqrt(x1i*x1i+xi1i*xi1i)*dt);fprintf(fp2,%fn,sqrt(x2i*x2i+xi2i*xi2i)*dt);3、 补到 128 位# include stdio.h# include math.h# include stdlib.h# include fft.c# define N 60# define PI 3.1415926# define dt 0.004 mai
22、n()void fft();double *x1,*xi1,*x2,*xi2; float z;int i,nfft;FILE *fp,*fp1,*fp2; nfft=128;x1=(double*)malloc(nfft*sizeof(double); xi1=(double*)malloc(nfft*sizeof(double); x2=(double*)malloc(nfft*sizeof(double); xi2=(double*)malloc(nfft*sizeof(double);for(i=0;iN;i+)x1i=sin(2*PI*(i+1)/30);fp=fopen(WAVE.DAT,r); for(i=0;iN;i+)fscanf(fp,%f,&z); x2i=z;for(i=0;iN;i+)xi1i=0.0;xi2i=0.0;/补到 128 位for(i=N;i128;i+)x1i=0.0;x2i=0.0;xi1i=0.0;xi2i=0.0;fft(x1,xi1,7,1);
