C
      SUBROUTINE PERIODOGRAM ( N, NSHIFT, NSAMP, T, X, Y, NSEG, PSD,
C
     1                         NPSD, W, XSEG, YSEG, WINDOW, TAPER )
C
C   This subroutine computes the averaged periodogram based on Eq.(5.43).
C   It computes the nearest integer number of segments, NSEG.  A Hamming
C   window is used, although other windows may be substituted from Table
C   5.1 .  NPSD power spectral density estimate values over the frequency
C   range -1/2T to 1/2T, where T is the sample interval in seconds, are
C   returned.  Some additional computational efficiency may be gained in
C   the autospectral estimate by eliminating the lines of code that deal
C   with data vector Y.  The PSD is complex, in general, for cross PSD.
C   The PSD is real for auto PSD.
C
C   Input Parameters:
C
C     N      - Number of data samples
C     NSHIFT - Number of samples shift between segments
C     NSAMP  - Number of samples per segment (must be even)
C     T      - Sample interval in seconds
C     X      - Array of complex samples X(1) to X(N)
C     Y      - Array of complex samples Y(1) to Y(N)
C     NPSD   - NO. FREQUENCIES IN THE PERIODOGRAM
C              ( MUST BE A POWER OF 2 )  --  CHANGE BY BWR 6/23/91
C     W      - WORKING VECTOR ( LENGTH = NPSD )
C     XSEG   - WORKING VECTOR ( LENGTH = NPSD )
C     YSEG   - WORKING VECTOR ( LENGTH = NPSD )
C     WINDOW - WORKING VECTOR ( LENGTH .GE. NSAMP )
C     TAPER  - LOGICAL VARIABLE ( .TRUE. TO USE TAPERING )
C
C   Output Parameters:
C
C     NSEG   - Number of segments averaged
C     PSD    - Real array of NPSD power spectral density estimate values
C              ( LENGTH = NPSD )
C
C   Notes:
C
C     External arrays X,Y must be dimensioned .GE. N and array PSD must be
C     dimensioned .GE. NPSD in the calling program.    Internal arrays W,
C     XSEG,YSEG must be dimensioned .GE. NPSD and array WINDOW must be
C     dimensioned .GE. NSAMP .  External subroutines used:  PREFFT and FFT
C     in Appendix 2.C .
C
C---------PARAMETER  ( NPSD = 4096 )     DISABLED BY BWR 6/23/91
C         ----------!  Must be a power of 2 -- 
C
       LOGICAL  TAPER
C
       COMPLEX  X(1),Y(1),PSD(1),W(NPSD),XSEG(NPSD),YSEG(NPSD)
C
       REAL  WINDOW(NSAMP)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
       PI2 = 8.0*ATAN(1.0)
C
       CALL PREFFT (NPSD,0,NEXP,W)
C
C   If TAPER .EQ. .TRUE. data is multiplied by a tapering window.
C   A Hamming window is used.  The two line calculation of WINDOW(K) in the
C   loop below may be replaced by any other desired window function.
C
      SUMW = 0.0
C                Accumulate the sum of the window ordinates for computing
C                          normalization factor used in DO 60 loop below.
      DO 10 K=1,NSAMP
       IF ( TAPER ) THEN
                    WINDOW(K)=0.538
     1                  +0.462*COS(PI2*(-.5+FLOAT(K-1)/FLOAT(NSAMP-1)))
              ELSE
                  WINDOW(K) = 1.0
       END IF
       SUMW = SUMW + WINDOW(K)
   10 CONTINUE
C
      UNORM = ( SUMW / FLOAT( NSAMP ) )**2
C
C   Compute the number of whole segments NSEG such that NSEG .LE. N
C
      NSEG=(N-NSAMP)/NSHIFT + 1
      DO 50 K=1,NSEG
        DO 20 J=1,NSAMP
          INDEX=J+(K-1)*NSHIFT
          XSEG(J)=X(INDEX)*WINDOW(J)
20        YSEG(J)=Y(INDEX)*WINDOW(J)
        DO 30 J=NSAMP+1,NPSD
          XSEG(J)=(0.,0.)
30        YSEG(J)=(0.,0.)
        CALL FFT (NPSD,0,T,NEXP,W,XSEG)
        CALL FFT (NPSD,0,T,NEXP,W,YSEG)
        DO 40 J=1,NPSD
          IF (K .EQ. 1)  PSD(J)=XSEG(J)*CONJG(YSEG(J))
40        IF (K .NE. 1)  PSD(J)=PSD(J)+XSEG(J)*CONJG(YSEG(J))
50      CONTINUE
C
C   Needed to add window energy correction factor
C
      DO 60 K=1,NPSD
         PSD(K) = PSD(K) / FLOAT( NSEG*NSAMP )
         PSD(K) = PSD(K) / UNORM
   60 CONTINUE
C----------------------60    PSD(K)=PSD(K)/FLOAT(NSEG*NSAMP)
C
      RETURN
      END
