C
      SUBROUTINE CORRELOGRAMPSD ( N, LAG, T, X, Y, PSD, NPSD, W,
C
     1                            R, WINDOW, MODE, TAPER )
C
C   Note that Rxy[-k] is computed as Conjugate(Ryx[k]).
C
C   Input Parameters:
C
C     N      - Number of data samples
C     LAG    - Highest lag index to compute
C     T      - Sample interval in seconds
C     X      - Array of complex data samples X(1) to X(N)
C     Y      - Array of complex data samples Y(1) to Y(N)
C     NPSD   - NO. FREQUENCIES IN THE CORRELOGRAM (MUST BE POWER OF TWO)
C     W      - COMPLEX WORKING VECTOR (LENGTH = NPSD)
C     R      - COMPLEX WORKING VECTOR (LENGTH = LAG)
C     WINDOW - REAL WORKING VECTOR (LENGTH = LAG)
C     MODE   - INTEGER ( MODE .EQ. 0 FOR UNBIASED CORRELATION ESTIMATE
C                        MODE .NE. 0 FOR BIASED CORRELATION ESTIMATE )
C     TAPER  - LOGICAL VARIABLE (TAPER .EQ. .TRUE. FOR TAPERING WINDOW)
C
C   Output Parameters:
C
C     PSD - Array of complex power spectral density estimate values
C
C   Notes:
C
C     External arrays X and Y must be dimensioned .GE. N and array PSD
C     must be dimensioned .GE. NPSD in the calling program.   Internal
C     array W must be dimensioned .GE. NPSD,  arrays R and WINDOW must
C     be dimensioned .GE. LAG. External subroutines used:  CORRELATION
C     (Appendix 5.A), PREFFT and FFT (Appendix 2.C).
C
      LOGICAL  TAPER
      COMPLEX  X(1), Y(1), PSD(1), W(NPSD), R(LAG), R0
      REAL  WINDOW(LAG)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      PI = ACOS( -1.0 )
C
C  Can use any window here
C
      DO 10 K=1,LAG
        IF ( TAPER ) THEN
                     WINDOW(K)=0.538+0.462*COS(PI*FLOAT(K)/FLOAT(LAG))
              ELSE
                  WINDOW(K) = 1.0
        END IF
   10 CONTINUE
C
c    ------------------------------------------------
      CALL CORRELATION ( N, LAG, MODE, X, Y, R0, R )
C    ------------------------------------------------
C
      PSD(1)=R0
      DO 20 K=1,LAG
20      PSD(K+1)=R(K)*WINDOW(K)
C
C    ------------------------------------------------
      CALL CORRELATION ( N, LAG, MODE, Y, X, R0, R )
C    ------------------------------------------------
      DO 30 K=1,LAG
30      PSD(NPSD+1-K)=CONJG(R(K))*WINDOW(K)
      DO 40 K=LAG+2,NPSD-LAG
40      PSD(K)=(0.,0.)
C
C    ---------------------------------------
      CALL PREFFT ( NPSD, 0, NEXP, W )
C
      CALL FFT ( NPSD, 0, T, NEXP, W, PSD )
C    ---------------------------------------
C
      RETURN
      END


