C PROGRAM 3.1 PERIOD C C ... This program computes periodogram ... 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 MJ: adjustable dimension of COV C NF: adjustable dimension of PE and SPE C IDEV: input device C JDEV: output device C IWINDW: window type (0: box-car, 1: Hanning, 2: Hamming) C @TEST.PN31: 5/16/89, 12/5/90, 12/21/90, 8/5/91 C PARAMETER( NMAX=200,MJ=NMAX,NF=512,IDEV=1,JDEV=6,IWINDW=1 ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), PE(0:NF), SPE(0:NF), COV(0:MJ) DATA OUTMIN/-1.0D30/, OUTMAX/ 1.0D30/ C CALL READTS( IDEV,Y,N ) C IF( IWINDW.EQ.0 ) LAG = MIN0( N-1,MJ ) IF( IWINDW.GT.0 ) LAG = 2*DSQRT( DBLE(N) ) IF( IWINDW.EQ.0 ) NP = (LAG+1)/2 IF( IWINDW.GT.0 ) NP = LAG C CALL PERIOD( Y,N,LAG,OUTMIN,OUTMAX,NP,0,COV,PE ) C CALL WINDOW( PE,NP,IWINDW,SPE ) C CALL PRPER( JEV,PE,SPE,N,NP,IWINDW ) C STOP E N D SUBROUTINE FOURIE( X,N,M,FC,FS ) C C ... Discrete Fourier transformation by Goertzel method ... C C Inputs: C X(I): data (I=1,N) C N: data length C M: number of Fourier components C FC(J): Fourier cosine transform (J=1,M) C FS(J): Fourier sine transform (J=1,M) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N), FC(M), FS(M) DATA PI/3.14159265358979D0/ C W = PI/(M-1) DO 20 I=1,M CI = DCOS(W*(I-1)) SI = DSIN(W*(I-1)) T1 = 0.0 T2 = 0.0 DO 10 J=N,2,-1 T0 = 2*CI*T1 - T2 + X(J) T2 = T1 10 T1 = T0 FC(I) = CI*T1 - T2 + X(1) 20 FS(I) = SI*T1 C RETURN E N D SUBROUTINE PERIOD( Y,N,LAG,OUTMIN,OUTMAX,NN,ISW,COV,P ) C C ... Periodogram computation by Fourier transf. of autocovariance ... C C Inputs: C Y(I): time series C N: data length C LAG: maximum lag of autocovariance C NN: number of frequencies C OUTMIN: lower bound for detecting outliers C OUTMAX: upper bound for detecting outliers C ISW = 0: autocovariance is not given C = 1: given C COV: autocovariance (if ISW=1) C Outputs: C COV(I): autocovariance (if ISW=0) C P(I): periodogram (raw spectrum) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N), COV(0:1024), P(0:512), FC(513), FS(513) C IF( LAG.GE.1023 ) LAG = 1023 IF( ISW.EQ.0 ) CALL AUTCOV( Y,N,LAG,OUTMIN,OUTMAX,COV ) C C ... Fourier transformation by Goertzel method ... C CALL FOURIE( COV(0),LAG+1,NN+1,FC,FS ) C DO 10 I=0,NN 10 P(I) = 2.0D0*FC(I+1) - COV(0) 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: output device specification 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: ',A72 ) 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 SUBROUTINE AUTCOV( Y,N,MAXLAG,OUTMIN,OUTMAX,C ) C C ... THIS SUBROUTINE COMPUTES SAMPLE AUTOCOVARIANCE ... C C INPUTS: C Y(I): TIME SERIES C N: DATA LENGTH C MAXLAG: MAXIMUM LAG OF AUTOCOVARIANCE C OUTMIN: BOUND FOR OUTLIERS IN LOW SIDE C OUTMAX: BOUND FOR OUTLIERS IN HIGH SIDE C OUTPUT: C C(I): AUTOCOVARIANCE FUNCTION C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION Y(N), C(0:MAXLAG ) C C ... SAMPLE MEAN ... C CALL MEAN( Y,N,OUTMIN,OUTMAX,NSUM,YMEAN ) C C ... SAMPLE AUTOCOVARIANCE ... C DO 20 LAG = 0,MAXLAG SUM = 0.0D0 DO 10 I=LAG+1,N IF( Y(I).GT.OUTMIN .AND. Y(I).LT.OUTMAX ) THEN IF( Y(I-LAG).GT.OUTMIN .AND. Y(I-LAG).LT.OUTMAX ) THEN SUM = SUM + ( Y(I)-YMEAN)*( Y(I-LAG) - YMEAN ) END IF END IF 10 CONTINUE 20 C(LAG) = SUM/NSUM C RETURN E N D SUBROUTINE MEAN( Y,N,OUTMIN,OUTMAX,NSUM,YMEAN ) C C ... This subroutine computes sample mean ... C C Inputs: C Y(I): time series C N: data length C OUTMIN: bound for outliers in low side C OUTMAX: bound for outliers in high side C Outputs: C NSUM: number of non-outlier observaions C YMEAN: sample mean C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N) C NSUM = 0 YSUM = 0.0D0 DO 10 I=1,N IF( Y(I) .GT.OUTMIN .AND. Y(I).LT.OUTMAX ) THEN NSUM = NSUM + 1 YSUM = YSUM + Y(I) END IF 10 CONTINUE YMEAN = YSUM/NSUM C RETURN E N D