C Find autocorrelation to el Nino data using direct and FFT methods.
C
INTEGER N
PARAMETER(N=168)
REAL EL(0:2*N-1),WSAVE(4*(2*N)+15),ACOV(0:N-1),A(2*N),B(2*N),
* ACOVR(0:N-1)
COMPLEX CEL(0:2*N-1),CORR(0:2*N-1)
LOGICAL EX
C
C needs N locations for el, 2N for acov, cel, corr
C and 4*(2N)+15 for wsave in complex case. N has maximum of 168.
INQUIRE (FILE='ELNINO.DAT',EXIST=EX,ERR=1000)
IF(.NOT.EX)GOTO 1000
OPEN (UNIT=8,FILE='ELNINO.DAT',ERR=1000)
C
C Read data, find mean.
SUM = 0.0
DO 100 I = 0,N-1
READ (8,*) EL(I)
SUM = SUM + EL(I)
100 CONTINUE
DO 101 I = 0,N-1
EL(I) = EL(I) - SUM/N
C Subtract mean, and add N zeros for either complex or real FFT usage
CEL(I) = CMPLX(EL(I),0.0)
EL(I+N) = 0.0
CEL(I+N) = 0.0
101 CONTINUE
C
C----------------------------------------------------------------------
C
C Direct calculation. Only sum as far as there is data.
C Simple, but slow.
C
DO 110 J = 0,N-1
ACOV(J) = 0.0
DO 110 M = 0,N-1-J
ACOV(J) = ACOV(J)+ EL(M) * EL(M+J)
110 CONTINUE
C
C Write, scaled correlation
WRITE (*,*)
* 'EX 11.6: AUTOCORRELATION (DIRECT) OUTPUT SUPRESSED'
CCCCC WRITE (*,'(5E14.6)') ( (ACOV(I)/ACOV(0)), I = 0,N-1)
C----------------------------------------------------------------------
C
C FFT approach (complex).
C Compute FFT of data of length 2N.
C Compute square of magnitude of transform components and place
C in complex array as real parts.
C Compute inverse transform, throwing away second half and
C imaginary parts (which are zero), and multiply by length of
C sequence, 2N.
CALL CFFTI (2*N,WSAVE)
CALL CFFTF (2*N,CEL,WSAVE)
C CFFTF returns unscaled transforms. Actual transforms are output
C divided by (2N).
DO 120 I = 0,2*N-1
CORR(I) = ABS(CEL(I) / (2*N)) **2
120 CONTINUE
C Since we compute transform times its conjugate, must divide by
C (2N) for each, i.e., (2N)**2.
CALL CFFTB (2*N,CORR,WSAVE)
C
DO 121 I = 0,N-1
ACOV(I) = REAL(CORR(I))*(2*N)
121 CONTINUE
C Autocovariance is inverse transform times sequence length, 2N.
C Normally, all the scaling would be done only once
C by dividing by 2N. We've broken it up for exposition.
WRITE (*,*)
* 'EX 11.6: AUTOCORRELATION (COMPLEX FFT) OUTPUT SUPRESSED'
CCCCC WRITE (*,'(5E14.6)') ( (ACOV(I)/ACOV(0) ), I = 0,N-1)
C
C----------------------------------------------------------------------
C
C FFT approach (real).
C Compute FFT of data of length 2N.
C EZFFTF produces correctly scaled A's and B's so no extra scaling
C needed to get transform.
C Compute array of square of each frequency component and place
C in cosine array (A's) to be back transformed. Set B's to 0.
C There are N A's, and N B's.
C Note that care must be taken to compute magnitude correctly,
C 0.5*(A(I)**2+B(I)**2) for I < N, twice that for I=N.
C Compute back transform throwing away its second half.
C
CALL EZFFTI (2*N,WSAVE)
CALL EZFFTF (2*N,EL,AZERO,A,B,WSAVE)
AZERO = AZERO*AZERO
C
DO 150 I = 1,N
IF(I.NE.N) THEN
A(I) = (A(I)**2 + B(I)**2) / 2.0
ELSE
A(I) = (A(I)**2 + B(I)**2)
ENDIF
B(I) = 0.0
150 CONTINUE
CALL EZFFTB (2*N,ACOVR,AZERO,A,B,WSAVE)
WRITE (*,*)
* 'EX 11.6: AUTOCORRELATION (REAL FFT) OUTPUT REDUCED'
WRITE (*,'(5E14.6)') ( (ACOVR(I)/ACOVR(0) ), I = 0,19)
CCCCC WRITE (*,'(5E14.6)') ( (ACOVR(I)/ACOVR(0) ), I = 0,N-1)
WRITE (*,*)
WRITE (*,*)
* 'REFERENCE RESULTS (EX 11.6 PARTIAL-LAST 4 LINES) FROM IBM PC/AT'
WRITE(*,*)' 0.100000E+01 0.606740E+00 0.353857E+00 0.184303E+00
* -0.152914E-01'
WRITE(*,*)'-0.219872E+00 -0.296218E+00 -0.291459E+00 -0.155056E+00
* 0.368081E-01'
WRITE(*,*)' 0.173559E+00 0.266598E+00 0.304989E+00 0.201169E+00
* 0.178037E-01'
WRITE(*,*)'-0.210367E+00 -0.377387E+00 -0.437960E+00 -0.460333E+00
* -0.425590E+00'
C
STOP
1000 WRITE (*,*) 'CANNOT FIND THE DATA FILE: ELNINO.DAT '
END