【正文】 DO 32 I=1,64 32 XSIN_64C(I)=CMPLX(XSIN(I),) CALL DFT(64,XSIN_64C,) DO 33 I=1,6433 XSIN_R(I)=REAL(XSIN_64C(I)) DO 34 I=1,6434 XSIN_I(I)=AIMAG(XSIN_64C(I)) DO 35 I=1,6435 XSIN_64S(I)=SQRT(XSIN_R(I)**2+XSIN_I(I)**2) WRITE(4,*)(XSIN_64S(I),I=1,64) CLOSE(4) OPEN(5,FILE=39。,STATUS=39。,FORM=39。) WRITE(3,*)(XWAVE_60ZS(I),I=1,60) CLOSE(3)CCCCCCCCCCCCC 對補零數據的分析 CCCCCCCCCCCCCCCCCCCCCCCCCCCC 對函數 SIN 的分析 CCCCCCCCCCCCCCCC OPEN(4,FILE=39。,STATUS=39。,FORM=39。) READ(2,*)(XWAVE(I),I=1,60) CLOSE(2) DO 21 I=1,6021 XWAVE_60C(I)=CMPLX(XWAVE(I),) CALL DFT(60,XWAVE_60C,) DO 22 I=1,6022 XWAVE_R(I)=REAL(XWAVE_60C(I)) DO 23 I=1,6023 XWAVE_I(I)=AIMAG(XWAVE_60C(I)) DO 24 I=1,6024 XWAVE_60ZS(I)=SQRT(XWAVE_R(I)**2+XWAVE_I(I)**2) OPEN(3,FILE=39。,STATUS=39。,FORM=39。) DO 1 I=1,60 1 XSIN(I)=SIN(2*PI*I/) DO 11 I=1,60 11 XSIN_60C(I)=CMPLX(XSIN(I),) CALL DFT(60,XSIN_60C,) DO 12 I=1,6012 XSIN_R(I)=REAL(XSIN_60C(I)) DO 13 I=1,6013 XSIN_I(I)=AIMAG(XSIN_60C(I)) DO 14 I=1,6014 XSIN_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=39。,STATUS=39。,FORM=39。) DO 3 I=1,1283 WAVESPRT(I)=SQRT(X_WAVE_R(I)**2+X_WAVE_I(I)**2) WRITE(4,*)(WAVESPRT(I), I=1,128) CLOSE(4) CLOSE(3) ENDCCCCCCCCCCCCCCC FFT 子程序 CCCCCCCCCCCCCCCCC SUBROUTINE FFT(LX,CX,SIGNI) COMPLEX CX(LX),CARG,CEXP,CW,CTEMP J=1 SC= IF()SC= SIG=SIGNI DO 30 I=1,LX IF()GO TO 10 CTEMP=CX(J)*SC CX(J)=CX(I)*SC CX(I)=CTEMP 10 M=LX/2 20 IF()GO TO 30 J=JM M=M/2 IF()GO TO 20 30 J=J+M L=1 40 ISTEP=2*L DO 50 M=1,L CARG=(,)*(*SIG*(M1))/L CW=CEXP(CARG) DO 50 I=M,LX,ISTEP CTEMP=CW*CX(I+L) CX(I+L)=CX(I)CTEMP 50 CX(I)=CX(I)+CTEMP L=ISTEP IF()GO TO 40 RETURN END6. 分析補零對振幅譜的影響CCCCCCCCCCC 分析補零對振幅譜的影響 CCCCCCCCCCC PROGRAM MAINCCCCCCCCCCC 不補零時的振幅普 CCCCCCCCCCCCCCCCCCCCCCCCCCCC 對函數 SIN 的分析 CCCCCCCCCCCCCCCC PARAMETER (PI=) 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_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=39。,STATUS=39。,FORM=39。) READ(3,*)(WAVE(I), I=1,128) DO 500 I=1,128500 X_WAVE_C(I)=CMPLX(WAVE(I),) CALL FFT(128,X_WAVE_C,) DO 600 I=1,128600 X_WAVE_R(I)=REAL(X_WAVE_C(I)) DO 700 I=1,128700 X_WAVE_I(I)=AIMAG(X_WAVE_C(I)) OPEN(4,FILE=39。,STATUS=39。,FORM=39。) DO 2 I=1,1282 XSPRT(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=39。,STATUS=39。,FORM=39。) DO 1 I=1,1281 SINSPRT(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)= DO 800 I=2,128800 X(I)= DO 900 I=1,128900 X_C(I)=CMPLX(X(I),) CALL FFT(128,X_C,) DO 1000 I=1,1281000 X_R(I)=REAL(X_C(I)) DO 1100 I=1,1281100 X_I(I)=AIMAG(X_C(I)) OPEN(5,FILE=39。,STATUS=39。,FORM=39。) K=LOG()/LOG() IF(2**)THEN K=K+1 ENDIF NFFT=2**K DO 10,I=50,50 IF()THEN H(I+51)=*F ELSE H(I+51)=SIN(2*PI**I*DT)/(PI*I*DT) END IF10 CONTINUE DO 20,I=102,12820 H(I)= DO 30,I=1,12830 H2(I)=CMPLX(H(I),) CALL FFT(128,H2,) DO 40,I=1,128 H2_R(I)=REAL(H2(I)) H2_I(I)=AIMAG(H2(I)) H_S(I)=(H2_R(I)**2+H2_I(I)**2)40 CONTINUE WRITE(2,*)(H_S(I),I=1,128) WRITE(1,*)(H(I),I=1,128) CLOSE(1) CLOSE(2) END CCCCCCCCCCCCC FFT 子程序 CCCCCCCCCCCCCC SUBROUTINE FFT(LX,CX,SIGNI) COMPLEX CX(LX),CARG,CEXP,CW,CTEMP J=1 SC= IF()SC= SIG=SIGNI DO 30 I=1,LX IF()GO TO 10 CTEMP=CX(J)*SC CX(J)=CX(I)*SC CX(I)=CTEMP 10 M=LX/2 20 IF()GO TO 30 J=JM M=M/2 IF()GO TO 20 30 J=J+M L=1 40 ISTEP=2*L DO 50 M=1,L CARG=(,)*(*SIG*(M1))/L CW=CEXP(CARG) DO 50 I=M,LX,ISTEP CTEMP=CW*CX(I+L) CX(I+L)=CX(I)CTEMP 50 CX(I)=CX(I)+CTEMP L=ISTEP IF()GO TO 40 RETURN END5. 頻譜分析CCCCCCCCCCCCC 頻譜分析 CCCCCCCCCCCCCC PROGRAM MAIN PARAM