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 ) CALL PTPER( SPE,N,LAG,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 PTPER( PE,N,LAG,NP,IWINDW ) C C ... This subroutine draws log-periodogram ... C C Inputs: C PE(I): log-periodogam C N: data length C LAG: maximum lag of autocovariance C NP: number of frequencies C IWINDW: window type (0: box-car, 1: Hanning, 2: Hamming) C IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 VNAME(5) DIMENSION PE(0:NP), VALUE(5) WY = 12.0 WX = 16.0 VNAME(1) = 'N = ' VNAME(2) = 'LAG = ' VNAME(3) = 'IWINDW =' VALUE(1) = N VALUE(2) = LAG VALUE(3) = IWINDW C CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 3.2: PERIODOGRAM BY FFT',33,3,VNAME,VALUE CALL MAXMIN( PE(1),NP,YMIN,YMAX,DY) CALL DRAWY2( 'PERIODOGRAM (IN LOG SCALE)',26,3.0D0,2.0D0, * PE(1),NP,WX,WY,1,2,0.0D0,0.5D0,0.1D0,YMIN,YMAX,DY ) CALL PLOTE C call plot( 0.0,0.0,999 ) 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 AXISXY(X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC,IBOX,IX,IY) C C ... This subroutine draws X and Y axes. C C Inputs: C X,Y: location of the left bottom of the figure C WX,WY: width and height of the figure C X0,X1: lower and upper bounds of the X axis C Y0,Y1: lower and upper bounds of the Y axis C DX,DY: step width in X and Y axes C DWC: size of characters C IBOX: = 0 draw X and Y axes only C = 1 draw window doundary C IX,IY: number of additional click in each step C Modified: 8/31/90, 11/21/90 C REAL*8 X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC C FX = X FY = Y FWX = WX FWY = WY FWC = DWC CALL NEWPEN( 2 ) CALL PLOT( FX,FY,-3 ) CALL PLOT( 0.0,FWY,3 ) CALL PLOT( 0.0,0.0,2 ) CALL PLOT( FWX,0.0,2 ) IF ( IBOX.EQ.1 ) THEN CALL PLOT( FWX,FWY,2 ) CALL PLOT( 0.0,FWY,2 ) END IF C C ... draw X axis ... C WC = DWC NX = (X1-X0)/DX + 0.001 NY = (Y1-Y0)/DY + 0.001 MX = -DLOG10( DX ) + 0.001 IF(DX.LT.1.0D0) MX = -DLOG10( DX ) + 0.999 IF( MX.LE.0 ) MX = -1 D = WX*DX/(X1-X0) DO 10 I=1,NX+1 XN = X0 + DX*(I-1) XX = D*(I-1) CALL PLOT( XX,-FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF LX = 0 IF(XN.GT.0.0D0) LX = ALOG10( XN ) + 0.001 IF( MX.GT.0 ) XXX = WC*(LX+MX)*0.5 + WC*4/7.0 IF( MX.LT.0 ) XXX = WC*LX*0.5 + WC*4/7.0 CALL NUMBER( XX-XXX,-(1.5*FWC+0.1),FWC,XN,0.0,MX ) 10 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 20 I=1,NX+1 DO 20 J=1,IX XX = D*(I-1) + (D*J)/IX IF( XX.LE.WX ) THEN CALL PLOT( XX,FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF END IF 20 CONTINUE CALL NEWPEN( 2 ) END IF C C ... draw Y axis ... C D = WY*DY/(Y1-Y0) MY = -DLOG10( DY ) + 0.99999 IF( MY.LE.0 ) MY = -1 DO 30 I=1,NY+1 YN = Y0 + DY*(I-1) YY = D*(I-1) CALL PLOT( -FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF LY = 1 IF(YN.LT. 0.0) LY = 2 IF(YN.GE. 1.0) LY = ALOG10( YN )+1.00001 IF(YN.LE.-1.0) LY = ALOG10(-YN )+2.00001 IF( MY.GT.0 ) YYY = WC*(LY+MY+1) + 0.20 IF( MY.LT.0 ) YYY = WC*LY + 0.20 CALL NUMBER( -YYY,YY-FWC/2,FWC,YN,0.0,MY ) 30 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 40 I=1,NY+1 DO 40 J=1,IY YY = D*(I-1) + (D*J)/IY IF( YY.LE.WY ) THEN CALL PLOT( FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF END IF 40 CONTINUE CALL NEWPEN( 2 ) END IF C RETURN E N D SUBROUTINE HEADER( PNAME,NP,M,VNAME,VALUE ) C C ... This subroutine draw title of data, program name, names C and values of parameters ... C C Inputs: C PNAME: program name C NP: number of characters of program name C M: number of parameters C VNAME(I): name of the I-th parameter C VALUE(I): value of the I-th parameter C Input via common: C TITLE: title of the data set C CHARACTER PNAME*80, TITLE*72 CHARACTER*8 VNAME(M) REAL*8 VALUE DIMENSION VALUE(M), DAY(2), TIME(3) COMMON /CMDATA/ TITLE C CALL NEWPEN( 2 ) CALL DATE( DAY ) CALL CLOCK( TIME,1 ) C CALL SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 ) CALL SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP ) CALL NEWPEN( 1 ) CALL SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 ) CALL SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 ) YY = 17.9 DO 10 I=1,M IF( MOD(I,4).EQ.1 ) THEN YY = YY - 0.40 X1 = 2.0 X2 = 3.5 ELSE X1 = X1 + 5.0 X2 = X2 + 5.0 END IF CALL SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 ) 10 CALL NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 ) C RETURN E N D SUBROUTINE DRAWY2( FNAME,NC,XO,YO,Y,N,WX,WY,IPOS,ISW,XMIN, * XMAX,DX,YMIN,YMAX,DY ) C C ... This subroutine draws time series data ... C C Inputs: C FNAME: title of the data set C NC: number of charcters of the name C XO,YO: origin of the figure C Y(I): time series C N: data length C WX,WY: width and height of the figure C IPOS: starting position of time series C ISW: = 1 connect data by line C = 2 connect data by dashed line C XMIN,XMAX: minimum and maximum of the X axis C DX: step width of the X axis C YMIN,YMAX: minimum and maximum of the Y axis C DY: step width of the Y axis C IMPLICIT REAL*8(A-H,O-Z) CHARACTER FNAME*80 DIMENSION Y(N) C CALL AXISXY(XO,YO,WX,WY,XMIN,XMAX,YMIN,YMAX,DX,DY,0.25D0,1,10,1) IF( ISW.LE.1 ) CALL PLOTY ( Y,N,YMIN,YMAX,WX,WY,IPOS,1 ) IF( ISW.EQ.1 ) CALL PLOTY2( Y,N,YMIN,YMAX,WX,WY,IPOS,1,0.0D0 ) IF( ISW.EQ.2 ) CALL PLOTY2( Y,N,YMIN,YMAX,WX,WY,IPOS,1,YMIN ) CALL NEWPEN( 1 ) CALL SYMBOL( 0.0,SNGL(WY)+0.2,0.3,FNAME,0.0,NC ) RETURN E N D SUBROUTINE MAXMIN( X,N,XMIN0,XMAX0,DXL ) C C ... This subroutine determines the minimum, the maximum and C the step width ... C C Inputs: C X(I): data C N: data length C Outputs: C XMIN0: the minimum bound for the figure C XMAX0: the maximum bound for the figure C DXL: step width C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N) C XMIN = 1.0D30 XMAX =-1.0D30 C DO 10 I=1,N IF( X(I) .LT. XMIN ) XMIN = X(I) 10 IF( X(I) .GT. XMAX ) XMAX = X(I) DX = XMAX-XMIN IF( DLOG10(DX) .GE. 0.0D0 ) DXL = INT( DLOG10(DX) ) IF( DLOG10(DX) .LT. 0.0D0 ) DXL = INT( DLOG10(DX) )-1.0 DXL = 10.0**DXL IF( DX/DXL.GT.6.0D0 ) DXL = DXL*2.0 DIF = INT( DX/DXL ) XMIN0 = INT( XMIN/DXL )*DXL XMAX0 = XMIN0 + DIF*DXL IF( XMIN0 .GT. XMIN ) XMIN0 = XMIN0 - DXL 30 IF( XMAX0 .GE. XMAX ) GO TO 40 XMAX0 = XMAX0 + DXL GO TO 30 40 CONTINUE C RETURN E N D SUBROUTINE PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C C ... This subroutine draws time vs data plot of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 connect data by straight line C = 2 connect data by dashed line C = 3 express data by character (character code=ISW) C Modified: 8/31/90 C REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IPOS) YY = (Y(I) - YMIN)*DY IF( I.EQ.1 .AND. ISW.LT.3 ) THEN CALL PLOT( XX,YY,3 ) ELSE IF(ISW.EQ.1) CALL PLOT( XX,YY,2 ) IF(ISW.EQ.2) CALL DASHP( XX,YY,0.2 ) IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) END IF 100 CONTINUE C RETURN E N D SUBROUTINE PLOTY2( Y,N,YMIN,YMAX,WX,WY,IOFF,ISW,YLEVEL ) C C ... This subroutine draws box car graph of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 draw box with thick line C = 2 draw box with thin line C YLEVEL: ground level for the box car C REAL*8 Y(N), YMIN, YMAX, WX, WY, YLEVEL C DX = WX/(N-1+IOFF) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 3 ) IF( N.GT.75 ) CALL NEWPEN( 2 ) IF( N.GT.150) CALL NEWPEN( 1 ) IF( ISW.EQ.2 ) CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IOFF) YY = (Y(I) - YMIN)*DY Y0 = (YLEVEL-YMIN)*DY CALL PLOT( XX,Y0,3 ) CALL PLOT( XX,YY,2 ) 100 CONTINUE CALL NEWPEN( 2 ) C RETURN E N D SUBROUTINE READTS( IDEV,Y,N ) REAL*8 Y(N) CHARACTER*72 TITLE COMMON /CMDATA/ TITLE C C OPEN( IDEV,FILE='temp.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N READ( IDEV,* ) (Y(I),I=1,N) C CLOSE( IDEV ) RETURN 1 FORMAT( A72 ) E N D