毕业设计-地震数据处理

《毕业设计-地震数据处理》由会员分享，可在线阅读，更多相关《毕业设计-地震数据处理（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 PARA

14、METER(DT=0.002)C H_NZ为非零相位滤波因子,H_Z为零相位滤波因子 CC Y_CNZ为非零相位滤波结果,Y_CZ为零相位滤波结果 CREAL Y_CNZ(1:100),Y_CZ(1:200) REAL H_Z(1:20),H_NZ(1:20) REAL X(1:50) INTEGER I,J,K,NREAL A(0:100),B(0:100)REAL Y_DNZ(1:100),Y_DZ(1:100)REAL X3(1:100),X4(1:100),X1(1:100),X2(1:100)REAL DFDATA H_NZ/1.0,5.254,11.6,13.7,8.47,0.77

15、5,-3.325,-2.764,-0.364, $ 1.099,0.97,0.15,-0.37,-0.344,-0.06,-0.126,0.122,0.025,-0.042, $ -0.043/ OPEN(1,FILE=INPUT3.DAT,FORM=FORMATTED,STATUS=UNKNOWN)READ(1,*)(X(I),I=1,50)CLOSE(1) CALL CON(X,H_NZ,Y_CNZ,50,20,69) OPEN(2,FILE=Y_CNZ.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(2,*)(Y_CNZ(i),I=1,69)CLOSE(2

16、) CALL AUTCOR(20,20,H_NZ,H_HZ) CALL con(X,H_CZ,Y_CZ,50,20,69) OPEN(3,FILE=Y_CZ.DAT,FORM=FORMATTED,STATUS=UNKNOWN)write(3,*)(Y_CZ(i),I=1,69) CLOSE(3)CLOSE(2)CCCCCCCCCCCCCCC 递归滤波的设计 CCCCCCCCCCCCCCCC X 存放INPUT3里数据的数组 CC Y_DNZ 为非零相位递归滤波，Y_DZ 为零相位递归滤波 CC X1,2,X,3,X4 为运算过程中的过度变量 CA(0)=1.0A(1)=4.0A(2)=6.0A

17、(3)=4.0A(4)=1.0b(0)=0.0B(1)=-1.254B(2)=0.987B(3)=-0.341B(4)=0.0524CCCCCCCCCCCCCCC 对 A 补零 CCCCCCCCCCCCCCCCCCDO 40 i=5,4940A(I)=0.0CCCCCCCCCCCCCCC 对 B 补零 CCCCCCCCCCCCCCCCCC DO 50 i=5,4950B(I)=0.0CCCCCCC 从外部数据向 X 导入数据 CCCCCCCCCCCCOPEN(1,FILE=INPUT3.DAT,FORM=FORMATTED,STATUS=UNKNOWN)READ(1,*)(X(I),I=1,5

18、0)CCCCCCCCCCC 非零相位递归滤波的设计 CCCCCCCCC DO 10 I=0,49 X1(I)=0.0X2(I)=0.0 DO 20 J=1,I20 X1(I)=X1(I)+A(J)*X(I-J)IF(I.EQ.0) THENX2(I)=0.0ELSEDO 30 K=1,I30X2(i)=X2(i)+B(k)*Y_DNZ(I-K) END IF Y_DNZ(I)=X1(I)-X2(I)10 CONTINUEOPEN(2,FILE=Y_DNZ.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(2,*)(Y_DNZ(I),I=1,49)CLOSE(2)C

19、CCCCCCCCC 零相位递归滤波的设计 CCCCCCCCCCCCDO 60 I=49,0,-1X3(I)=0.0X4(I)=0.0DO 70 J=0,49-I70X3(I)=X3(I)+A(J)*Y_DNZ(I+J) IF(I.EQ.49) THENX4(I)=0.0ELSEDO 80 K=2,49-I80 X4(I)=X4(I)+B(K)*Y_DZ(I+K) END IF Y_DZ(I)=X3(I)-X4(I)60 CONTINUE OPEN(3,FILE=Y_DZ.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(3,*)(Y_DZ(I),I=1,49)CL

20、OSE(3)CLOSE(1)ENDCCCCCCCCCCCCCCCCCCC 褶积子程序 CCCCCCCCCCCCCCCCCCC SUBROUTINE 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.004RETURNENDCCCCCCCCCCCCCCCCCC 自相关子程序 CCCCCCCCCCCCCCCCCC SUBROUTINE AUTCOR(J,K,C,A)DIMENSION C(J),A(K)DO

21、 7 N=1,KA(N)=0.0I=NDO 6 IN=I,JIL=IN-N+1 6A(N)=A(N)+C(IN)*C(IL) 7CONTINUERETURNEND4. 设计高通滤波因子CCCCCCCCCCC 高通滤波器的设计 CCCCCCCCCCC PROGRAM MAIN PARAMETER (N=101,DT=0.004,PI=3.141592654,F=30.0)REAL H(150),H2_R(128),H2_I(128),H_S(128)INTEGER K,NFFTCOMPLEX H2(128)OPEN(1,FILE=H1_HIGHPASS.DAT,FORM=FORMATTED,ST

22、ATUS=UNKNOWN) OPEN(2,FILE=H2_HIGHPASS.DAT,FORM=FORMATTED,STATUS=UNKNOWN)K=LOG(101.0)/LOG(2.0) IF(2*K.LT.101)THENK=K+1 ENDIF NFFT=2*KDO 10,I=-50,50IF(I.EQ.0)THENH(I+51)=1.0/DT-2*FELSEH(I+51)=-SIN(2*PI*30.0*I*DT)/(PI*I*DT)END IF10 CONTINUE DO 20,I=102,12820H(I)=0.0 DO 30,I=1,12830H2(I)=CMPLX(H(I),0.0)

23、CALL FFT(128,H2,1.0) DO 40,I=1,128H2_R(I)=REAL(H2(I)H2_I(I)=AIMAG(H2(I)H_S(I)=(H2_R(I)*2+H2_I(I)*2)40CONTINUEWRITE(2,*)(H_S(I),I=1,128)WRITE(1,*)(H(I),I=1,128)CLOSE(1) CLOSE(2)ENDCCCCCCCCCCCCC FFT 子程序 CCCCCCCCCCCCCCSUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1

24、.0)SC=1.0SIG=-SIGNIDO 30 I=1,LXIF(I.GT.J)GO TO 10CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP 10 M=LX/2 20IF(J.LE.M)GO TO 30J=J-MM=M/2IF(M.GE.1)GO TO 20 30J=J+ML=1 40ISTEP=2*LDO 50 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1)/LCW=CEXP(CARG)DO 50 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP 50CX(I)=CX(I)+

25、CTEMPL=ISTEPIF(L.LT.LX)GO TO 40RETURNEND5. 频谱分析CCCCCCCCCCCCC 频谱分析 CCCCCCCCCCCCCC PROGRAM MAIN PARAMETER (PI=3.141592654,DT=0.004)REAL XSIN(200),X(200),WAVE(200)REAL X_SIN_R(200),X_SIN_I(200) REAL X_R(200),X_I(200) REAL X_WAVE_R(200),X_WAVE_I(200)REAL SINSPRT(200),XSPRT(200),WAVESPRT(200)REAL DFCOMPL

26、EX XSIN_C(200),X_C(200),X_WAVE_C(200) CCCCCCCCCCCC 对 SIN 的分析 CCCCCCCCCCCCDO 100 I=1,128XSIN(I)=SIN(2*I*PI/64.0)100CONTINUECCCCCCCCC 是数据转换成复数形式 CCCCCCCCCC DO 200 I=1,128200XSIN_C(I)=CMPLX(XSIN(I),0.0)CALL FFT(128,XSIN_C,1.0) DO 300 I=1,128 300 X_SIN_R(I)=REAL(XSIN_C(I)DO 400 I=1,128400X_SIN_I(I)=AIMA

27、G(XSIN_C(I) OPEN(1,FILE=XSINSPRT.DAT,FORM=FORMATTED,STATUS=UNKNOWN) DO 1 I=1,1281SINSPRT(I)=SQRT(X_SIN_R(I)*2+X_SIN_I(I)*2)WRITE(1,*)(SINSPRT(I), I=1,128)CLOSE(1)CCCCCCCCCCCCCC 对脉冲信号的分析 CCCCCCCCCCCCC X(1)=1.0DO 800 I=2,128800X(I)=0.0 DO 900 I=1,128900X_C(I)=CMPLX(X(I),0.0)CALL FFT(128,X_C,1.0) DO 10

28、00 I=1,1281000X_R(I)=REAL(X_C(I)DO 1100 I=1,1281100X_I(I)=AIMAG(X_C(I) OPEN(5,FILE=XSPRT.DAT,FORM=FORMATTED,STATUS=UNKNOWN)DO 2 I=1,1282XSPRT(I)=SQRT(X_R(I)*2+X_I(I)*2)WRITE(5,*)(XSPRT(I), I=1,128)CLOSE(5)CCCCCCCCCCCC 对地震波 WAVE 的分析 CCCCCCCCCCC OPEN(3,FILE=WAVE.DAT,FORM=FORMATTED,STATUS=UNKNOWN)READ(

29、3,*)(WAVE(I), I=1,128) DO 500 I=1,128500X_WAVE_C(I)=CMPLX(WAVE(I),0.0) CALL FFT(128,X_WAVE_C,1.0)DO 600 I=1,128600X_WAVE_R(I)=REAL(X_WAVE_C(I)DO 700 I=1,128700X_WAVE_I(I)=AIMAG(X_WAVE_C(I) OPEN(4,FILE=WAVESPRT.DAT,FORM=FORMATTED,STATUS=UNKNOWN) DO 3 I=1,1283 WAVESPRT(I)=SQRT(X_WAVE_R(I)*2+X_WAVE_I(I

30、)*2)WRITE(4,*)(WAVESPRT(I), I=1,128)CLOSE(4)CLOSE(3)ENDCCCCCCCCCCCCCCC FFT 子程序 CCCCCCCCCCCCCCCCCSUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1.0)SC=1.0SIG=-SIGNIDO 30 I=1,LXIF(I.GT.J)GO TO 10CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP 10 M=LX/2 20IF(J.LE.M)GO TO 30

31、J=J-MM=M/2IF(M.GE.1)GO TO 20 30J=J+ML=1 40ISTEP=2*LDO 50 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1)/LCW=CEXP(CARG)DO 50 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP 50CX(I)=CX(I)+CTEMPL=ISTEPIF(L.LT.LX)GO TO 40RETURNEND6. 分析补零对振幅谱的影响CCCCCCCCCCC 分析补零对振幅谱的影响 CCCCCCCCCCC PROGRAM MAINCCCCCCCCCCC 不补零时的

32、振幅普 CCCCCCCCCCCCCCCCCCCCCCCCCCCC 对函数 SIN 的分析 CCCCCCCCCCCCCCCC PARAMETER (PI=3.1415) REAL XSIN(150),XWAVE(150)REAL XWAVE_R(150),XWAVE_I(150)REAL XSIN_Z(150),XWAVE_Z(150)REAL XSIN_R(150),XSIN_I(150) COMPLEX XSIN_60C(150),XWAVE_60C(150) COMPLEX XSIN_64C(150),XWAVE_64C(150) COMPLEX XSIN_128C(150),XWAVE_

33、128C(150)REAL XSIN_60S(150),XWAVE_60ZS(150)REAL XSIN_64S(150),XWAVE_64ZS(150)REAL XSIN_128S(150),XWAVE_128ZS(150)OPEN(1,FILE=XSIN_60S.DAT,FORM=FORMATTED,STATUS=UNKNOWN)DO 1 I=1,60 1XSIN(I)=SIN(2*PI*I/30.0) DO 11 I=1,60 11XSIN_60C(I)=CMPLX(XSIN(I),0.0)CALL DFT(60,XSIN_60C,1.0)DO 12 I=1,6012XSIN_R(I)=

34、REAL(XSIN_60C(I) DO 13 I=1,6013XSIN_I(I)=AIMAG(XSIN_60C(I) DO 14 I=1,6014XSIN_60S(I)=SQRT(XSIN_R(I)*2+XSIN_I(I)*2)WRITE(1,*)(XSIN_60S(I),I=1,60) CLOSE(1)CCCCCCCCCCCCC 对地震数据的分析 CCCCCCCCCCCCCC OPEN(2,FILE=WAVE.DAT,FORM=FORMATTED,STATUS=UNKNOWN)READ(2,*)(XWAVE(I),I=1,60) CLOSE(2)DO 21 I=1,6021XWAVE_60C

35、(I)=CMPLX(XWAVE(I),0.0) CALL DFT(60,XWAVE_60C,1.0)DO 22 I=1,6022XWAVE_R(I)=REAL(XWAVE_60C(I) DO 23 I=1,6023XWAVE_I(I)=AIMAG(XWAVE_60C(I) DO 24 I=1,6024XWAVE_60ZS(I)=SQRT(XWAVE_R(I)*2+XWAVE_I(I)*2) OPEN(3,FILE=XWAVE_60ZS.DAT,FORM=FORMATTED,STATUS=UNKNOWN) WRITE(3,*)(XWAVE_60ZS(I),I=1,60)CLOSE(3)CCCCC

36、CCCCCCCC 对补零数据的分析 CCCCCCCCCCCCCCCCCCCCCCCCCCCC 对函数 SIN 的分析 CCCCCCCCCCCCCCCC OPEN(4,FILE=XSIN_64S.DAT,FORM=FORMATTED,STATUS=UNKNOWN)DO 31 I=61,64 31XSIN(I)=0.0 DO 32 I=1,64 32XSIN_64C(I)=CMPLX(XSIN(I),0.0)CALL DFT(64,XSIN_64C,1.0)DO 33 I=1,6433XSIN_R(I)=REAL(XSIN_64C(I) DO 34 I=1,6434XSIN_I(I)=AIMAG(

37、XSIN_64C(I) DO 35 I=1,6435XSIN_64S(I)=SQRT(XSIN_R(I)*2+XSIN_I(I)*2)WRITE(4,*)(XSIN_64S(I),I=1,64)CLOSE(4)OPEN(5,FILE=XSIN_128S.DAT,FORM=FORMATTED,STATUS=UNKNOWN)DO 41 I=65,12841XSIN(I)=0.0 DO 42 I=1,128 42XSIN_128C(I)=CMPLX(XSIN(I),0.0)CALL DFT(128,XSIN_128C,1.0)DO 43 I=1,12843XSIN_R(I)=REAL(XSIN_12

38、8C(I) DO 44 I=1,12844XSIN_I(I)=AIMAG(XSIN_128C(I) DO 45 I=1,12845XSIN_128S(I)=SQRT(XSIN_R(I)*2+XSIN_I(I)*2)WRITE(5,*)(XSIN_128S(I),I=1,128)CLOSE(5)CCCCCCCCCCCCC 对地震数据的分析 CCCCCCCCCCCCCCOPEN(6,FILE=XWAVE_64ZS.DAT,FORM=FORMATTED,STATUS=UNKNOWN) DO 51 I=61,6451 XWAVE(I)=0.0DO 52 I=1,6452XWAVE_64C(I)=CMP

39、LX(XWAVE(I),0.0) CALL DFT(64,XWAVE_64C,1.0)DO 53 I=1,6453XWAVE_R(I)=REAL(XWAVE_64C(I) DO 54 I=1,6454XWAVE_I(I)=AIMAG(XWAVE_64C(I) DO 55 I=1,6455XWAVE_64ZS(I)=SQRT(XWAVE_R(I)*2+XWAVE_I(I)*2) WRITE(6,*)(XWAVE_64ZS(I),I=1,64)CLOSE(6)OPEN(7,FILE=XWAVE_128ZS.DAT,FORM=FORMATTED,STATUS=UNKNOWN) DO 61 I=65,

40、12861 XWAVE(I)=0.0DO 62 I=1,12862XWAVE_128C(I)=CMPLX(XWAVE(I),0.0) CALL DFT(128,XWAVE_128C,1.0)DO 63 I=1,12863XWAVE_R(I)=REAL(XWAVE_128C(I) DO 64 I=1,12864XWAVE_I(I)=AIMAG(XWAVE_128C(I) DO 65 I=1,12865XWAVE_128ZS(I)=SQRT(XWAVE_R(I)*2+XWAVE_I(I)*2) WRITE(7,*)(XWAVE_128ZS(I),I=1,128)CLOSE(7)ENDCCCCCCC

41、CCCCCCCC DFT 子程序 CCCCCCCCCCCCCCCCSUBROUTINE DFT(N,C,S)COMPLEX C(N),Y(512),WIF(S.EQ.1.0)Z=1.0IF(S.EQ.-1.0)Z=1.0/FLOAT(N)DO 20 I=1,NY(I)=(0.0,0.0)DO 20 J=1,NA=COS(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N)B=SIN(2.0*3.14159265*(I-1)*(J-1)/FLOAT(N)*(-S)W=CMPLX(A,B) 20Y(I)=Y(I)+C(J)*WDO 30 I=1,N 30C(I)=Y(I)*ZRE

42、TURNEND7. 线性褶积与循环褶积CCCCCCCCC 线性滤波与循环褶积比较 CCCCCCCCC PROGRAM MAIN CCCCCCCCCCCCC 线性滤波 CCCCCCCCCCCCCCCCCCCCCCCCCCC H是低通滤波因 CCCCCCCCCCCCPARAMETER (PI=3.141592654,DT=0.002,F=70.0)REAL X(100),X2(101),H(-50:50),Y_LINEAR(200),Y_CIRCLE(200)OPEN(1,FILE=INPUT1.DAT,FORM=FORMATTED,STATUS=UNKNOWN)READ(1,*)(X(I),I=

43、1,100)CCCCCCCCCCCCCCC低通滤波器CCCCCCCCCCCCC DO 10 I=-50,50 IF (I.EQ.0)THEN H(I)=2*F*PI/PI ELSE H(I)=SIN(2*PI*F*I*DT)/(PI*I*DT) END IF10CONTINUECALL CON(X,H,Y_LINEAR,100,101,200)CCCCCCCCCCCC输出滤波后的数据CCCCCCCCCCCCCCCOPEN(2,FILE=Y_LINEAR.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(2,*)(Y_LINEAR(I),I=1,200)CLOSE(

44、2)CCCCCCCCCCCCC 圆周滤波的设计 CCCCCCCCCCCC OPEN(3,FILE=Y_CIRCLE.DAT,FORM=FORMATTED,STATUS=UNKNOWN) DO 20 I=1,10020X2(I)=X(I)X2(100)=0.0CALL CIRCON(X2,H,Y_CIRCLE,101) WRITE(3,*)(Y_CIRCLE(I),I=1,101) CLOSE(3)CCCCCCC 线性卷积输出点等于圆周卷积 CCCCCC OPEN(4,FILE=Y_LINEARCUT.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(4,*)(Y_

45、LINEAR(I),I=1,101)CLOSE(4)ENDCCCCCCCCCCCCCC 线性卷积子程序 CCCCCCCCCCCCCCSUBROUTINE 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.004RETURNENDCCCCCCCCCCCCCCC 圆周卷积子程序 CCCCCCCCCCCCCCSUBROUTINE CIRCON(X,H,Y,N)REAL X(N),H(N),Y(N)DO I=1

46、,NY(I)=0.0ENDDODO I=1,N DO J=1,N IF(I-J.GT.0) Y(I)=Y(I)+X(J)*H(I-J) ENDDO ENDDOEND8. 最小平方反滤波CCCCCCCCCC 两种反射系数序列的对比 CCCCCCCCCC PROGRAM MAIN REAL B(12),R(200),X_CONV(211),RXX(211) REAL A(211),A_CONV(270)DATA B/17.5,12.5,7.0,0.0,-3.5,-7.0,-3.0,-1.0,0.0,1.4,3.5,2.0/ OPEN(1,FILE=INPUT8.DAT,FORM=FORMATTED

47、,STATUS=UNKNOWN)READ(1,*)(R(I),I=1,200)CCCCCCCCCCCCCC 求合成地震记录 CCCCCCCCCCCCCCCALL CON(R,B,X_CONV,200,12,211)CCCCCCCCCCCC 求地震子波的自相关 CCCCCCCCCCCCCALL AUTCOR(211,211,X_CONV,RXX)CCCCCCCCCCCCCCC 求反滤波因子 CCCCCCCCCCCCCCCCALL PEO(211,RXX,A)OPEN(2,FILE=A.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(2,*)(A(I),I=1,60

48、)CCCCCCC 通过反滤波因子求反射系数序列 CCCCCCCCALL CON(A,X_CONV,A_CONV,60,211,270)OPEN(3,FILE=A_CONV.DAT,FORM=FORMATTED,STATUS=UNKNOWN)WRITE(3,*)(A_CONV(I),I=1,270)CLOSE(3)CLOSE(2)CLOSE(1)ENDCCCCCCCCCCCCCCC 褶积子程序 CCCCCCCCCCCCCCCCSUBROUTINE 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

49、 2 I2=1,JII=I1+I2-1 2C(II)=C(II)+A(I1)*B(I2)*0.004RETURNENDCCCCCCCCCCCCCCC 自相关子程序 CCCCCCCCCCCCCCCSUBROUTINE AUTCOR(J,K,C,A)DIMENSION C(J),A(K)DO 7 N=1,KA(N)=0.0I=NDO 6 IN=I,JIL=IN-N+1 6A(N)=A(N)+C(IN)*C(IL) 7CONTINUERETURNENDCCCCCCCCCCCCCC 求反滤波因子 CCCCCCCCCCCCCCCCSUBROUTINE PEO(LR,R,A)CPEO IS THE AUX

50、ILIARY TOEPLITZ RECURSIONCIT GIVES THE PREDICTION-ERROR OPERATOR A(I)DIMENSION R(LR),A(LR)V=R(1)A(1)=1.0IF(LR.EQ.1)RETURNDO 6 L=2,LRD=0.0L3=L-1DO 1 J=1,L3K=L-J+1 1D=D+A(J)*R(K)C=D/VIF(L.EQ.2)GOTO 4L1=(L-2)/2L2=L1+1IF(L2.LT.2)GOTO 3DO 2 J=2,L2HOLD=A(J)K=L-J+1A(J)=A(J)-C*A(K) 2A(K)=A(K)-C*HOLDIF(2*L1.

51、EQ.L-2)GOTO 4 3LT3=L2+1A(LT3)=A(LT3)-C*A(LT3) 4A(L)=-C 6V=V-C*DRETURNEND9. 零相位转换CCCCCCCCCCC 零相位的转换 CCCCCCCCCCCCC PROGRAM MAINREAL B(32),B1(32),BC_S(32),B_CONV(32)REAL BC_R(32),BC_I(32)DATA B/0,1,3,5,7,9,10,8,6,4,2,0,-1,0,1,2,3,1.5,0,-1,0,2,1,0,0/ COMPLEX B_C(32),BC_SZ(32) DO 10,I=1,32IF(I.GE.26)THEN

52、 B_C(I)=(0.0,0.0)ELSE B_C(I)=CMPLX(B(I),0.0) END IF10CONTINUECALL FFT(32,B_C,1.0) DO 20,I=1,32 BC_R(I)=REAL(B_C(I)BC_I(I)=AIMAG(B_C(I)BC_S(I)=SQRT(BC_R(I)*2+BC_I(I)*2)20BC_SZ(I)=CMPLX(BC_S(I),0.0) CALL FFT(32,BC_SZ,-1.0)DO 30,I=1,3230B_CONV(I)=REAL(BC_SZ(I) OPEN (1,FILE=B_CONV.DAT,FORM=FORMATTED,STA

53、TUS=UNKNOWN)WRITE(1,*)(B_CONV(I),I=1,32)CLOSE(1)ENDCCCCCCCCCCCCCC FFT子程序 CCCCCCCCCCCCCCSUBROUTINE FFT(LX,CX,SIGNI)COMPLEX CX(LX),CARG,CEXP,CW,CTEMPJ=1SC=1.0/LXIF(SIGNI.EQ.1.0)SC=1.0SIG=-SIGNIDO 30 I=1,LXIF(I.GT.J)GO TO 10CTEMP=CX(J)*SCCX(J)=CX(I)*SCCX(I)=CTEMP 10 M=LX/2 20IF(J.LE.M)GO TO 30J=J-MM=M/

54、2IF(M.GE.1)GO TO 20 30J=J+ML=1 40ISTEP=2*LDO 50 M=1,LCARG=(0.0,1.0)*(3.14159265*SIG*(M-1)/LCW=CEXP(CARG)DO 50 I=M,LX,ISTEPCTEMP=CW*CX(I+L)CX(I+L)=CX(I)-CTEMP 50CX(I)=CX(I)+CTEMPL=ISTEPIF(L.LT.LX)GO TO 40RETURNEND10. 最小相位转换CCCCCCCCCCCC 最小相位转换 CCCCCCCCCCCC PROGRAM MAIM REAL B(25),B0(75)DATA B/0,1,2,3,1.5,0,-1,0,1,3,5,7,9,10,8,6,4,2,0,-1,0,2,3,2,0/ REAL RBB(25),Y0(25),A0(25)REAL SUMX,C(49),G(49)OPEN(1,FILE=B0.DAT,FORM=FORMATTED,STATUS=UNKNOWN) OPEN(2,FILE=B.D