C PROGRAM 3.2 FFTPER C C ... This program computes periodogram via FFT ... C C The following inputs are required in READTS. C TITLE: title of the data C N: data length C Y(I): time series C Parameters: C NMAX: adjustable dimension of Y (NMAX.GE.N) C NF: adjustable dimension of PE and SPE C MJ: adjustable dimension of COV C @TEST.PN31: 5/16/89, 12/5/90, 12/21/90, 8/5/91, 7/19/92 C PARAMETER( NMAX=3000,NF=512,IDEV=1,JDEV=6,IWINDW=1 ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), PE(0:NF), SPE(0:NF) C CALL READTS( IDEV,Y,N ) C IF( IWINDW.EQ.0 ) LAG = N-1 IF( IWINDW.GT.0 ) LAG = 2*DSQRT( DFLOAT(N) ) IF( IWINDW.EQ.0 ) NP = (LAG+1)/2 IF( IWINDW.GT.0 ) NP = LAG C CALL FFTPER( Y,N,NN,PE,NP ) CALL WINDOW( PE,NP,IWINDW,SPE ) CALL PRPER( JDEV,PE,SPE,N,NP,IWINDW ) C STOP E N D SUBROUTINE FFTPER( Y,N,NN,P,N2 ) C C ... periodogram computation by FFT ... C C Inputs: C Y(I): data C N: data length C NN: basic span for computing peiodogram (if NN>0) C Outputs: C P(J): periodogram C NN: basic span for computing peiodogram C N2: number of frequencies C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N),P(0:1024),X(1024),FY(1024),WK(1024) C IF(NN.LE.0) THEN IF( N.LE.1024 ) THEN NP = ALOG( N-0.01 )/ALOG(2.0) + 1 NN = 2**NP NPOOL = 1 ELSE NN = 1024 NPOOL = (N-1)/NN + 1 END IF ELSE NP = ALOG(NN-0.01)/ALOG(2.0) +1 NN = 2**NP IF( NN.GT.1024) NN=1024 NPOOL = (N-1)/NN +1 END IF C N2 = NN/2 DO 10 I=0,N2 10 P(I) = 0.0D0 C DO 100 II=1,NPOOL ND = MIN(N,II*NN) - (II-1)*NN DO 20 I=1,ND 20 X(I) = Y((II-1)*NN+I) DO 30 I=ND+1,NN 30 X(I) = 0.0D0 C CALL FFTR2 ( X,NN,1,FY,WK ) C P(0) = P(0) + FY(1)**2 P(N2) = P(N2) + FY(N2+1)**2 DO 40 I=0,N2-1 40 P(I) = P(I) + FY(I+1)**2 + FY(N2+1)**2 100 CONTINUE C DO 50 I=0,N2 50 P(I) = P(I)/N C RETURN E N D SUBROUTINE FFTR2( X,N,ISW,FX,WRK ) C C ... Fast Fourier Transform ... C C Inputs: C X(I): data C N: data length C ISW: C WRK: working area C Output: C FX(J): Fourier transform C J=1,...,N/2+1 cosine transform C J=N/2+2,...,N sine transform C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N), FX(N), WRK(N/4) DATA PI2 /6.2831853071796D0/ C NP = DLOG( DFLOAT(N) )/DLOG( 2.0D0 ) + 1.0D-5 N2 = N/2 N4 = N/4 C IF( ISW.EQ.1 ) THEN DO 10 I=2,N4 10 WRK(I) = DSIN( (I-1)*PI2/N ) END IF C DO 20 I=1,N2 FX(I) = X(I) + X(I+N2) 20 FX(I+N2) = X(I) - X(I+N2) C IFG = 1 JFG = 1 K = 1 M = N4 C DO 30 I=1,NP-1 M2 = M*2 K2 = K*2 IF( M.GE.K ) THEN IF( IFG.LT.0 ) THEN CALL FFTSB1( X,WRK,K,M,M2,K2,FX ) ELSE CALL FFTSB1( FX,WRK,K,M,M2,K2,X ) END IF ELSE IF( JFG.EQ.1 ) THEN IF( IFG.LT.0 ) THEN CALL FFTSB2( X,K2,M2,FX ) ELSE CALL FFTSB2( FX,K2,M2,X ) END IF JFG = 2 END IF IF( IFG.LT.0 ) THEN CALL FFTSB3( FX,WRK,K,M,X ) ELSE CALL FFTSB3( X,WRK,K,M,FX ) END IF END IF K = K*2 M = M/2 30 IFG = -IFG C IF( IFG.GT.0 ) THEN DO 40 I=1,N 40 FX(I) = X(I) END IF C RETURN E N D SUBROUTINE FFTSB2( X,M,L,Y ) C C ... slave subroutine for FFTR2 ... C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(L,M), Y(M,L) C IF( M.GE.L ) THEN DO 10 J=1,L DO 10 I=1,M 10 Y(I,J) = X(J,I) ELSE DO 20 I=1,M DO 20 J=1,L 20 Y(I,J) = X(J,I) END IF C RETURN E N D SUBROUTINE FFTSB1( X,SINE,K,M,MJ1,MJ2,Y ) C C ... slave subroutine for FFTR2 ... C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,MJ2), Y(M,K,4), SINE(M,K) C DO 10 I=1,M Y(I,1,1) = X(I,1) + X(I+M,1) Y(I,1,3) = X(I,1) - X(I+M,1) Y(I,1,2) = X(I,K+1) 10 Y(I,1,4) = X(I+M,K+1) C DO 20 J=2,K DO 20 I=1,M SUM1 = SINE(1,K-J+2)*X(I+M,J) - SINE(1,J) *X(I+M,J+K) SUM2 = SINE(1,J) *X(I+M,J) + SINE(1,K-J+2)*X(I+M,J+K) Y(I,J,1) = X(I,J) + SUM1 Y(I,K-J+2,2) = X(I,J) - SUM1 Y(I,J,3) = SUM2 + X(I,J+K) 20 Y(I,K-J+2,4) = SUM2 - X(I,J+K) C RETURN E N D SUBROUTINE FFTSB3( X,SINE,K,M,Y ) C C ... slave subroutine for FFTR2 ... C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(K,2,M,2), SINE(M,K), Y(K,4,M) C DO 10 J=1,M Y(1,1,J) = X(1,1,J,1) + X(1,1,J,2) Y(1,3,J) = X(1,1,J,1) - X(1,1,J,2) Y(1,2,J) = X(1,2,J,1) Y(1,4,J) = X(1,2,J,2) DO 10 I=2,K SUM1 = SINE(1,K-I+2)*X(I,1,J,2) - SINE(1,I) *X(I,2,J,2) SUM2 = SINE(1,I) *X(I,1,J,2) + SINE(1,K-I+2)*X(I,2,J,2) Y(I,1,J) = X(I,1,J,1) + SUM1 Y(K-I+2,2,J) = X(I,1,J,1) - SUM1 Y(I,3,J) = SUM2 + X(I,2,J,1) 10 Y(K-I+2,4,J) = SUM2 - X(I,2,J,1) C RETURN E N D SUBROUTINE PRPER( JDEV,PE,SPE,N,MJ,IWINDW ) C C ... This subroutine prints out periodogram ... C C Inputs: C JDEV: C PE: periodogram (raw spectrum) C SPE: logarithm of smoothed periodogram C N: data length C MJ: number of frequencies C IWINDW: type of specral window C CHARACTER*72 TITLE REAL*8 PE(0:MJ), SPE(0:MJ) COMMON /CMDATA/ TITLE C WRITE(JDEV,600) WRITE(JDEV,610) TITLE WRITE(JDEV,620) N, MJ, IWINDW WRITE(JDEV,630) (PE(I),I=0,MJ) WRITE(JDEV,640) (SPE(I),I=0,MJ) C 600 FORMAT( 1H ,'PROGRAM 3.1: ' ) 610 FORMAT( 1H ,' DATA: ',10A8 ) 620 FORMAT( 1H ,' N =',I5,3X,'MJ =',I4,3X,'IWINDW =',I2 ) 630 FORMAT( 1H0,'** PERIODOGRAM: P(I);(I=0,...,MJ) **',/ * (1H ,5D13.5) ) 640 FORMAT( 1H0,'** LOG-PERIODOGRAM: LOG P(I);(I=0,...,MJ) **',/ * (1H ,10F8.3) ) RETURN E N D SUBROUTINE WINDOW( PE,NP,IWINDW,SPE ) C C ... Smoothing by spectral window and log-transformation ... C C Inputs: C PE(I): raw specrum C NP: number of frequencies C IWINDW: window type (0: box-car, 1: Hanning, 2: Hamming) C Outputs: C SPE(I): logarithm of smoothed peiodogram C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PE(0:NP), SPE(0:NP) DIMENSION W(0:1,2) DATA W/0.5D0, 0.25D0, 0.54D0, 0.23D0/ C IF( IWINDW.EQ.0 ) THEN PMIN = 1.0D30 DO 10 I=0,NP 10 IF( PE(I).GT.0.0D0 .AND. PE(I).LT.PMIN ) PMIN = PE(I) DO 20 I=0,NP 20 SPE(I) = DMAX1( PE(I),PMIN ) ELSE SPE(0) = W(0,IWINDW)*PE(0) + W(1,IWINDW)*PE(1)*2 SPE(NP)= W(0,IWINDW)*PE(NP)+ W(1,IWINDW)*PE(NP-1)*2 DO 30 I=1,NP-1 30 SPE(I) = W(0,IWINDW)*PE(I) + W(1,IWINDW)*(PE(I-1) + PE(I+1)) END IF DO 50 I=0,NP 50 SPE(I) = DLOG10( SPE(I) ) C RETURN E N D SUBROUTINE READTS( IDEV,Y,N ) REAL*8 Y(N) CHARACTER*72 TITLE COMMON /CMDATA/ TITLE C OPEN( IDEV,FILE='temp.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N READ( IDEV,* ) (Y(I),I=1,N) CLOSE( IDEV ) RETURN 1 FORMAT( A72 ) E N D