C SPECTRAL.FOR - COMPUTATION OF THE DISCRETE FOURIER TRANSFORM C AND THE PERIODOGRAM OF A TIME SERIES. THE PERIODOGRAM IS C COMPUTED AS THE SQUARED AMPLITUDE R(W)**2 RATHER THAN THE C ACTUAL PERIODOGRAM (N/(8*PI))R(W)**2 FOR COMPUTATIONAL USE C IN SUBSEQUENT PROCESSING. C DIMENSION X(1024),Y(1024) CHARACTER * 40 DF WRITE(1,2050) 2050 FORMAT(/" SPECTRAL - Fourier Transform And Periodogram") 2001 WRITE(1) "Enter UPPERCASE Input File Name: " READ(1) DF IF(IOREAD(5,2,0,DF)) GOTO 2001 M=5 DO 10 I=1,1024 X(I)=0.0 10 Y(I)=0.0 CALL DATIN(X,N,START,STEP,M) WRITE(1,7001) 7001 FORMAT(' Order of Trend to Remove (0=Mean, 1=Linear): '$ READ(1,7002) K 7002 FORMAT(I0) CALL DETRND(X,N,K) WRITE(1,7003) 7003 FORMAT(' Proportion of Data to Taper: '$ READ(1,7004) P 7004 FORMAT(F0.0) CALL TAPER(X,N,P) NP2=1 20 NP2=NP2*2 IF(NP2.LT.N) GOTO 20 INV=0 CALL FT01A(X,Y,NP2,INV) IF(INV.EQ.0) GOTO 30 WRITE(1) "*** ERROR IN SUBROUTINE FT01A ***" GOTO 99 30 NPGM=NP2/2+1 2052 WRITE(1) "Enter UPPERCASE File for Transform Output: " READ(1) DF IF(IOWRIT(8,2,0,DF)) GOTO 2052 DF="REAL PART OF TRANSFORM" CALL DATOUT(X,NPGM,0.0,1.0,8,DF) DF="IMAGINARY PART OF TRANSFORM" CALL DATOUT(Y,NPGM,0.0,1.0,8,DF) INV=IOCLOS(8) CON=(2.0*FLOAT(NP2)/FLOAT(N))**2 DO 40 I=1,NPGM 40 X(I)=(X(I)**2+Y(I)**2)*CON 2053 WRITE(1) "Enter UPPERCASE File for Periodogram Output: " READ(1) DF IF(IOWRIT(8,2,0,DF)) GOTO 2053 DF="SQUARED AMPLITUDE" CALL DATOUT(X,NPGM,0.0,1.0,8,DF) INV=IOCLOS(8) 99 STOP END C SUBROUTINE FT01A(XR,XI,N,INV) C C THIS SUBROUTINE IMPLEMENTS THE SANDE-TUKEY RADIX-2 FAST FOURIER C TRANSFORM. EITHER THE DIRECT OR INVERSE TRANSFORM MAY BE COMP- C UTED. PARAMETERS ARE: C C XR - REAL ARRAY - REAL PART OF THE - RETURNS REAL PART C SERIES OF TRANSFORM C XI - REAL ARRAY - IMAGINARY PART OF - RETURNS IMAGINARY C THE SERIES PART OF TRANSFORM C N - INTEGER - LENGTH OF SERIES - (SAME) C INV - INTEGER - 0 FOR DIRECT - RETURNS ERROR CODE C 1 FOR INVERSE (-1 IF ERROR) C C DIRECT IS: C N C (1/N) SUM X(T)*EXP(-2*PI*I*(T-1)*(J-1)/N), J=1,...,N C T=1 C C INVERSE IS: C N C SUM X(T)*EXP( 2*PI*I*(T-1)*(J-1)/N), J=1,...,N C T=1 C DIMENSION XR(N),XI(N),UR(15),UI(15) LOGICAL FIRST DATA FIRST /.TRUE./ IF(.NOT.FIRST) GOTO 120 UR(1)=0.0 UI(1)=1.0 DO 110 I=2,15 UR(I)=SQRT((1.0+UR(I-1))/2.0) 110 UI(I)=UI(I-1)/(2.0*UR(I)) FIRST=.FALSE. 120 IF((N.GT.0).AND.(FLOAT(N).LE.2.0**16)) GOTO 130 INV=-1 RETURN 130 N0=1 II=0 140 N0=N0+N0 II=II+1 IF(N0.LT.N) GOTO 140 I1=N0/2 I3=1 I0=II WRITE(1,7001) 7001 FORMAT(' FFTing data pass - '$ DO 260 I4=1,II WRITE(1,7002) I4 7002 FORMAT(' ',I2$ DO 250 K=1,I1 WR=1.0 WI=0.0 KK=K-1 DO 230 I=1,I0 IF(KK.EQ.0) GOTO 240 IF(MOD(KK,2).EQ.0) GOTO 230 J0=I0-I WS=WR*UR(J0)-WI*UI(J0) WI=WR*UI(J0)+WI*UR(J0) WR=WS 230 KK=KK/2 240 IF(INV.EQ.0) WI=-WI L=K DO 250 J=1,I3 L1=L+I1 ZR=XR(L)+XR(L1) ZI=XI(L)+XI(L1) Z=WR*(XR(L)-XR(L1))-WI*(XI(L)-XI(L1)) XI(L1)=WR*(XI(L)-XI(L1))+WI*(XR(L)-XR(L1)) XR(L1)=Z XR(L)=ZR XI(L)=ZI 250 L=L1+I1 I0=I0-1 I3=I3+I3 260 I1=I1/2 WRITE(1,7003) 7003 FORMAT(' *') UM=1.0 IF(INV.EQ.0) UM=1.0/FLOAT(N0) DO 310 J=1,N0 K=0 J1=J-1 DO 320 I=1,II K=2*K+MOD(J1,2) 320 J1=J1/2 K=K+1 IF(K.LT.J) GOTO 310 ZR=XR(J) ZI=XI(J) XR(J)=XR(K)*UM XI(J)=XI(K)*UM XR(K)=ZR*UM XI(K)=ZI*UM 310 CONTINUE RETURN END C SUBROUTINE DETRND(X,N,K) C C COMPUTES AND SUBTRACTS EITHER THE SERIES MEAN OR A LINEAR C LEAST SQUARES TREND FROM THE SERIES DATA, E.G. POLYMONIAL C OF DEGREE K C C X - REAL ARRAY - THE INPUT SERIES - RETURNS DETRENDED C N - INTEGER - SERIES LENGTH - (UNCHANGED) C K - INTEGER - DEGREE OF POLYN - (UNCHANGED) C DIMENSION X(N) WRITE(1,7001) 7001 FORMAT(' Detrending data...') SUMX=0.0 DO 20 I=1,N 20 SUMX=SUMX+X(I) XBAR=SUMX/FLOAT(N) WRITE(1,7002) XBAR 7002 FORMAT(5X,'Mean Removed = ',F13.5) DO 30 I=1,N 30 X(I)=X(I)-XBAR IF(K.LE.0) RETURN TBAR=FLOAT(N+1)/2.0 SUMTT=(FLOAT(N)*(FLOAT(N)**2-1))/12.0 SUMTX=0.0 DO 40 I=1,N 40 SUMTX=SUMTX+X(I)*(FLOAT(I)-TBAR) BETA=SUMTX/SUMTT WRITE(1,7003) BETA 7003 FORMAT(5X,'Slope Removed = ',F13.5) DO 50 I=1,N 50 X(I)=X(I)-BETA*(FLOAT(I)-TBAR) RETURN END C SUBROUTINE TAPER(X,N,P) C C APPLIES SPLIT-COSINE-BELL TAPERING TO THE TIME SERIES. THE C ARGUMENT P IS THE TOTAL PROPORTION OF THE SERIES TO TAPER. C TUKEY SUGGESTS .1 TO .2, .25 IS THE DEFAULT IN MANY PACKAGES. C C X - REAL ARRAY - THE INPUT SERIES - RETURNS TAPERED C N - INTEGER - SERIES LENGTH - (UNCHANGED) C P - REAL - PROP. TO TAPER - (UNCHANGED) C DATA PI /3.141593/ DIMENSION X(N) IF((P.LE.0.0).OR.(P.GE.1.0)) RETURN WRITE(1,7001) 7001 FORMAT(' Tapering data...') M=INT(P*FLOAT(N)+0.5)/2 DO 10 I=1,M WEIGHT=0.5-0.5*COS(PI*(FLOAT(I)-0.5)/FLOAT(M)) X(I)=X(I)*WEIGHT 10 X(N+1-I)=X(N+1-I)*WEIGHT RETURN END C SUBROUTINE DATIN(X,N,START,STEP,M) C C INPUT A SERIES OF VALUES IN RUNTIME FORMAT. THE FIRST FOUR C CARDS MUST BE (FREE FORMAT): C 1) A TITLE CARD (<72 COLS.) C 2) NUMBER OF CASES (I5) C 3) DATA FORMAT (<72 COLS.) C 4) START TIME VALUE AND STEP (2F10.5) C C X - REAL ARRAY - RETURN THE SERIES C N - INTEGER - RETURN SERIES LENGTH C START - REAL - RETURN SERIES TIME VALUE AT 1ST POINT C STEP - REAL - RETURN TIME INCREMENT BETWEEN POINTS C M - INTEGER - UNIT NUMBER (UNCHANGED) C DIMENSION X(1024) CHARACTER * 72 HEAD,FMT WRITE(1,7001) 7001 FORMAT(' Reading data...') READ(M,2) HEAD,N,FMT,START,STEP 2 FORMAT(A0/I5/A0/2F10.5) WRITE(1,3) HEAD,N,FMT,START,STEP 3 FORMAT(5X,'The Input Data Header Reads:'/5X,A0/ * 5X,'Series Length Input is = ',I5/5X,'Format Of Input Data:'/ * 5X,A0/5X,'Time Origin = ',F11.5,' Time Increment = ',F11.5) READ(M,FMT) (X(I),I=1,N) RETURN END C SUBROUTINE DATOUT(X,N,START,STEP,M,HEAD) C C OUTPUT A SERIES OF VALUES IN FORMAT REQUIRED FOR DATIN. C C X - REAL ARRAY - THE SERIES C N - INTEGER - SERIES LENGTH C START - REAL - SERIES TIME VALUE AT 1ST POINT C STEP - REAL - TIME INCREMENT BETWEEN POINTS C M - INTEGER - UNIT NUMBER (UNCHANGED) C HEAD - CHAR*40 - HEADER TITLE C DIMENSION X(N) CHARACTER * 40 HEAD WRITE(1,7001) 7001 FORMAT(' Writing data...') WRITE(M,1) HEAD,N,START,STEP 1 FORMAT(1X,A0/1X,I5/' (5E15.8)'/1X,2F10.5) WRITE(M,2) (X(I),I=1,N) 2 FORMAT(1X,5E15.8) RETURN END