C....*...1.........2.........3.........4.........5.........6.........7.*.......8 C DPER 10/2/74 C C PURPOSE C COMPUTE THE PERIODOGRAM OF A TIME SERIES C C USAGE C CALL DPER(X,N,P,NP,IER) C C ARGUMENTS C X - INPUT VECTOR OF LENGTH N C REAL*8 C N - LENGTH OF X C INTEGER*4 C P - COMPUTED PERIODOGRAM ORDINATES. P MUST BE DIMENSIONED C LARGER THAN N+2 TO ALLOW WORKSPACE FOR THE FFT ALGORITHM C REAL*8 C NP - COMPUTED NUMBER OF PERIODOGRAM ORDINATES. C INTEGER*4 C IER - ERROR PARAMETER C IER=0 FFT ALGORITHM EMPLOYED C IER=1 FFT ALGORITHM NOT EMPLOYED C C SUBROUTINE DPER(X,N,P,NP,IER) IMPLICIT REAL*8 (A-H,O-Z) save REAL*8 X(1),P(1) IF(N.EQ.1) GO TO 60 MODN2=MOD(N,2) IF(MODN2.EQ.1)GO TO 40 M=N/2 NP=M+1 DO 10 I=1,M P(I)=X(2*I-1) 10 P(NP+I)=X(2*I) ISN=1 CALL FFT(P(1),P(NP+1),M,M,M,ISN) IF(ISN.EQ.0) GO TO 40 IER=0 CALL REALTR(P(1),P(NP+1),M,ISN) FN=DFLOAT(N) C=1.D0/(2.D0*FN) DO 30 I=1,NP 30 P(I)=C*(P(I)**2+P(NP+I)**2) RETURN 40 IF(MODN2.EQ.0) NP=N/2+1 IF(MODN2.EQ.1)NP=(N-1)/2+1 IER=1 FN=DFLOAT(N) C=2.D0/FN TWOPI=8.D0*DATAN(1.D0) DO 50 K=1,NP WK=TWOPI*DFLOAT(K-1)/FN CALL DCXSX(X,N,WK,CX,SX) 50 P(K)=C*(CX*CX+SX*SX) RETURN 60 P(1)=2.D0*X(1)**2 IER=1 RETURN END