《地震数据处理课程设计-毕业论文.doc》由会员分享，可在线阅读，更多相关《地震数据处理课程设计-毕业论文.doc（37页珍藏版）》请在金锄头文库上搜索。
1、地震资料数据处理方法课程设计报告 地震资料数据处理课程设计总结报告专业班级： 姓 名： 学 号： 设计时间： 指导老师： 目 录一 、 设计内容（1）褶积滤波（2）快变滤波（3）褶积滤波与快变滤波的比较（4）设计高通滤波因子（5）频谱分析（6）分析补零对振幅谱的影响（7）线性褶积与循环褶积（8）最小平方反滤波（9）零相位转换（10）最小相位转换（11）静校正二、附录（1）附录1：相关程序 （2）附录2：相关图件【附录1：有关程序】1. 褶积滤波CCCCCCCCCCCCCCCCC 褶积滤波 CCCCCCCCCCCCCCCCC PROGRAM MAINDIMENSION X(100),H1(-50
2、:50),H2(-50:50),Y_LOW(200),Y_BAND(200)PARAMETER (PI=3.141592654)CCCCCCCC H1是低通滤波因子，H2为带通滤波因子CCCCCCREAL X,H1,H2,Y_LOW,Y_BANDREAL dt,F,F1,F2INTEGER Idt=0.002F=70.0F1=10.0F2=80.0OPEN(1,FILE=INPUT1.DAT,FORM=FORMATTED,STATUS=UNKNOWN)READ(1,*)(X(I),I=1,100)CCCCCCCCCCCCCCCCCC低通滤波器CCCCCCCCCCCCCCCCCDO 10 I=-
3、50,50 IF (I.EQ.0)THEN H1(I)=2*F*PI/PI ELSE H1(I)=SIN(2*PI*F*I*dt)/(PI*I*dt) END IF10CONTINUECCCCCCCCCCCCCCCC输出低通滤波因子CCCCCCCCCCCCCCCCOPEN(2,FILE=H1_LOW.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(2,*)(H1(I),I=-50,50)CLOSE(2)CALL CON(X,H1,Y_LOW,100,101,200)CCCCCCCCCCCCCCCC输出滤波后的数据CCCCCCCCCCCCCCCCOPEN(3,FI
4、LE=Y_LOW.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(3,*)(Y_LOW(I),I=51,150)CLOSE(3)CCCCCCCCCCCCCCCCCC带通滤波器CCCCCCCCCCCCCCCCCCCC DO 20 I=-50,50 IF(I.EQ.0)THEN H2(I)=140 ELSE H2(I)=SIN(2*PI*F2*I*dt)/(PI*I*dt)-SIN(2*PI*F1*I*dt)/(PI*I*dt) END IF 20CONTINUECCCCCCCCCCCCCCC输出带通滤波因子CCCCCCCCCCCCCCCCCOPEN(4,FILE=
5、H2_BAND.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(4,*)(H2(I),I=-50,50) CLOSE(4)CALL CON(X,H2,Y_BAND,100,101,200)CCCCCCCCCCCCCCCC输出滤波后的数据CCCCCCCCCCCCCCCCCOPEN(5,FILE=Y_BAND.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(5,*)(Y_BAND(I), I=51,150)CLOSE(5) ENDCCCCCCCCCCCCCCCCCCCCC褶积函数CCCCCCCCCCCCCCCCCCCCSUBROUTI
6、NE CON(A,B,C,I,J,K)DIMENSION A(I),B(J),C(K)DO 1 K1=1,K 1C(K1)=0.0DO 2 I1=1,IDO 2 I2=1,JII=I1+I2-1 2C(II)=C(II)+A(I1)*B(I2)*0.002RETURNEND2. 快变滤波CCCCCCCCCCCCCCC频率滤波CCCCCCCCCCCCCCCCCCCC PROGRAM MAINPARAMETER (PI=3.141592654)REAL H_LOW(1:200),H_BAND(1:200),Y_LOW(1:200),Y_BAND(1:200) REAL X(1:200)INTEGE
7、R I,NFFT,K,NPCOMPLEX C(1:200),TEMP(1:200)REAL DT,DF,F,F1,F2F=70.0F1=10.0F2=80.0DT=0.002OPEN(1,FILE=INPUT1.DAT,FORM=FORMATTED,STATUS=UNKNOWN)READ(1,*)(X(I), I=1,100)NP=100K=LOG(100.0)/LOG(2.0)IF(2*K.LT.100)THENK=K+1 NFFT=2*KEND IF DF=1.0/(NFFT*DT)DO 10 I=101,128 X(I)=0.010CONTINUE CCCCCCCCCCC 数据变换成复数
8、形式 CCCCCCCCCCC DO 20 I=1,NFFT20C(I)=CMPLX(X(I),0.0)CCCCCCCCCCC 向频率域转换 CCCCCCCCCCCCCCCCCCALL FFT(NFFT,C,1.0)CCCCCCCCCCC 低通滤波因子的设计 CCCCCCCCCCC DO 30 I=1,NFFT/2IF(I*DF.LE.F)THENH_LOW(I)=1.0 ELSEH_LOW(I)=0.0END IF30CONTINUECCCCCCCCC 使滤波器理想化（对称）CCCCCCCCCC DO 40 I=NFFT/2+1,NFFT F=I*DFH_LOW(I)=H_LOW(NFFT-I
9、+1)40CONTINUECCCCCCCCCCCCCCC 实现滤波 CCCCCCCCCCCCCCCCC DO 50 I=1 , NFFT50TEMP(I)=C(I)*H_LOW(I)CCCCCCCCCCCCCCC 向时间域变换 CCCCCCCCCCCCC CALL FFT(NFFT,TEMP,-1.0)DO 60 I=1 , NFFT60 Y_LOW(I)=REAL(TEMP(I)OPEN(2,FILE=LOWPASS.DAT,FORM=FORMATTED,STATUS=UNKNOWN) WRITE(2,*)(Y_LOW(I),I=1,NFFT)CLOSE(2)CCCCCCCCCCC 带通滤波
10、因子的设计 CCCCCCCCCCC DO 70 I=1,NFFT/2IF(I*DF.GE.F1).AND.(I*DF.LE.F2)THENH_BAND(I)=1.0 ELSEH_BAND(I)=0.0END IF70CONTINUECCCCCCCCC 使滤波器理想化（对称）CCCCCCCCCC DO 80 I=NFFT/2+1,NFFT F=I*DFH_LOW(I)=H_BAND(NFFT-I+1)80CONTINUECCCCCCCCCCCCCCC 实现滤波 CCCCCCCCCCCCCCCCC DO 90 I=1 , NFFT90TEMP(I)=C(I)*H_BAND(I)CCCCCCCCCC
11、CCCCC 向时间域变换 CCCCCCCCCCCCC CALL FFT(NFFT,TEMP,-1.0)DO 100 I=1 , NFFT100 Y_BAND(I)=REAL(TEMP(I)OPEN(3,FILE=BANDPASS.DAT,FORM=FORMATTED,STATUS=UNKNOWN) WRITE(3,*)(Y_BAND(I),I=1,NFFT)CLOSE(3)CLOSE(1)ENDCCCCCCCCCCCCCCCCC子程序CCCCCCCCCCCCCCCCCCCCC LX 为输入数据的点数 CCCCCCCCCCCCCCCCCCCCCC CX 为复型数组 CCCCCCCCCCCCCCC
12、CCCCCCCCCCCC SIGNI 为转换标志 CCCCCCCCCCCCCCCCCCCC SUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1.0)SC=1.0SIG=-SIGNIDO 300 I=1,LXIF(I.GT.J)GO TO 100CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP 100 M=LX/2 200IF(J.LE.M)GO TO 300J=J-MM=M/2IF(M.GE.1)GO TO 200 300J=J+ML=1 400
13、ISTEP=2*LDO 500 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1)/LCW=CEXP(CARG)DO 500 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP 500CX(I)=CX(I)+CTEMPL=ISTEPIF(L.LT.LX)GO TO 400RETURNEND3. 褶积滤波与快变滤波的比较CCCCCCCCCCCCC 褶积滤波与递归滤波的比较 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 褶积滤波的设计 CCCCCCCCCCCCCCCCCC PROGRAM MAIN PARAME
