C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++                                           C
       SUBROUTINE PERDGR ( X, N, DELT, NF, WORK, FREQ, PERDG, AMP,
C
     1                     SUMP, PAVE, PMAX, FPMAX, FISHER, PMSIG,
C                                 
     2                     CUMPER, NOUT, FLB, CLB, FUB, CUB, WHN, INFO )
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  Subroutine to calculate amplitude sepctrum, periodogram, and cumulative 
C             periodogram of a real time-series, using P.N. Swartztrauber's
C             fast Fourier transform subroutine EZFFTF.
C
C  Written by:  Bert W. Rust
C               Computing and Applied Mathematics Laboratory
C               National Institue of Standards and Technology
C               June 1995
C-------------------------------------------------------------------------------
C
C  Input:
C  ------
C  X       Real array (length = N)
C          Time-series to be analyzed
C
C  N       Integer (N .GT. 0)
C          No. terms in the time-series X
C
C  DELT    Real
C          Sample spacing for the time-series
C
C  NF      Integer 
C          No. frequencies desired in spectrum estimate (including freq. = 0)
c           (      .GE. N/2 + 1 , if N is even )
C           ( NF <                             )
C           (      .GE. (N+1)/2 , if N is odd  )
C
C  WORK    Real array (length = 10*NF + 9)
C          Working vector
C
C Output:
C -------
C  FREQ    Real array (length = NF)
C          Equally spaced spectral frequencies with: 
C          FREQ(1) = 0
C                       1 / (2*DELT) ,         if 2*(NF - 1) .GE. N
C          FREQ(NF) = < 
C                       (1 - 1/N) / (2*DELT) , if 2*(NF - 1) .LT. N,
C                                              i.e. if N odd and NF = (N+1)/2
C
C  PERDG   Real array (length = NF)
C          Periodogram estimate
C
C  AMP     Real array (length = NF)
C          Amplitude spectrum estimate
C
C  SUMP    Real
C          Sum of the periodogram ordinates
C
C  PAVE    Real
C          Average value of the periodogram ordinates
C
C  PMAX    Real
C          The maximum periodogram ordinate
C
C  FPMAX   Real
C          Frequency at which the maximum periodogram ordinate occurs
C
C  FISHER  Fisher statistic for testing (at the 95% level) the significance 
C          of the maximum periodogram ordinate
C
C  PMSIG   Logical variable
C          Indicates whether the maximum periodogram ordinate is significantly
C          greater (at the 95% level) than the average value.
C          .TRUE.  if  ( PMAX / PAVE ) .GT. FISHER
C          .FALSE. if  ( PMAX / PAVE ) .LE. FISHER
C
C  CUMPER  Real array (length = NF)
C          Cumulative periodogram estimate
C
C  NOUT    Integer ( 0 .LE. NOUT .LE. NF )
C          Number of points at which the cumululative periodogram ordinate
C          falls outside the 95% confidence band for white noise
C
C  FLB     Real array (length = 2)
C          Endpoint abcissas for the lower bound line of the 95% confidence 
C          band for a cumulative periodogram of white noise
C
C  CLB     Real array (length = 2)
C          Endpoint ordinates for the lower bound line of the 95% confidence 
C          band for a cumulative periodogram of white noise
C
C  FUB     Real array (length = 2)
C          Endpoint abcissas for the upper bound line of the 95% confidence 
C          band for a cumulative periodogram of white noise
C
C  CUB     Real array (length = 2)
C          Endpoint ordinates for the upper bound line of the 95% confidence 
C          band for a cumulative periodogram of white noise
C
C  WHN     Logical variable
C          Indicates whether the cumulative periodogram is consistent, at the
C          95% level, with a white noise spectrum:
C          .TRUE. if  ( NOUT / NF ) .LE. 0.05
C          .FALSE. if ( NOUT / NF ) .GT. 0.05
C
C  INFO    Integer
C          Error status indicator:
C                                 INFO  OK       Description
C                                 ---------------------------------
C                                   0   Yes  Successfull execution
C                                   1   No   N .LE. 0
C                                   2   No   N is even and NF .LT. N/2 + 1
C                                   3   No   N is odd  and NF .LT. (N+1)/2
C
C Subroutines called:  EZFFTI, EZFFTF
c -------------------
C
C-------------------------------------------------------------------------------
C
      EXTERNAL  EZFFTI, EZFFTF
C
      LOGICAL  PMSIG, WHN, NPEVEN
C
      REAL  X(N), WORK(10*NF+9)
      REAL  FREQ(NF), PERDG(NF), AMP(NF), CUMPER(NF)
      REAL  FLB(2), CLB(2), FUB(2), CUB(2)  
C
C-------------------------------------------------------------------------------
c
      IF ( N .LE. 0 ) THEN
                      INFO = 1
C                    ++++++++
                      RETURN
C                    ++++++++
      END IF
C
      IF ( 2*(N/2) .EQ. N ) THEN
                                IF ( 2*(NF - 1) .LT. N ) THEN
                                                             INFO = 2
C                                                           ++++++++
                                                             RETURN
C                                                           ++++++++
                                END IF
                   ELSE
                       IF ( (2*NF - 1) .LT. N ) THEN
                                                    INFO = 3
C                                                  --------
                                                    RETURN
C                                                  --------
                       END IF
      END IF
C
      INFO = 0
C
C  Zero-pad the data vector (if necessary) and compute the Fourier transform
C
      DO 20 I=1,N
         WORK(I) = X(I)
   20 CONTINUE
C
      NP = 2*(NF - 1)
      IF ( NP .LT. N ) NP = N
C
      IF ( 2*(NP/2) .EQ. NP ) THEN
                                  NPEVEN = .TRUE.
                  ELSE
                      NPEVEN = .FALSE.
      END IF
C
      IF ( NP .GT. N ) THEN
         DO 30 I = N+1, NP
            WORK(I) = 0.0
   30    CONTINUE
      END IF
C
C    ++++++++++++++++++++++++++++++++++
      CALL EZFFTI ( NP, WORK(2*NF) )
C    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      CALL EZFFTF ( NP, WORK(1), A0, WORK(8*NF+12),  
     1                               WORK(9*NF+11), WORK(2*NF) )
C    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  Compute the periodogram and amplitude spectrum
C
      RCPNPT = 1.0 / ( NP*DELT )
C
      DO 50 K = 1, NF
         FREQ(K) = (K - 1) * RCPNPT
   50 CONTINUE
C 
      CFAC = FLOAT( NP ) / FLOAT( N ) 
C        
      AMP(1) = CFAC*ABS( A0 )
      PERDG(1) = N * AMP(1)**2
C
      IF ( NPEVEN ) THEN
                        KMAX = NF - 2
              ELSE
                  KMAX = NF - 1
      END IF
C
      DO 70 K = 2, KMAX+1
         AMP(K) = CFAC * SQRT(   WORK(8*NF+K+10)**2
     1                         + WORK(9*NF+K+ 9)**2  )
         PERDG(K) = N * AMP(K)**2 / 4.0
   70 CONTINUE
C
      IF ( NPEVEN ) THEN
                        AMP(NF) = CFAC * ABS( WORK(9*NF+10) )
                        PERDG(NF) =  N * AMP(NF)**2
      END IF
C
C  Find the maximum peridogram ordinate and compute Fisher statistic
C
      PSUM = 0.0
      PMAX = 0.0
      IMAX = 1
      DO 80 I=1,NF
         PSUM = PSUM + PERDG(I)
         IF ( PERDG(I) .GT. PMAX )  THEN
                                    PMAX = PERDG(I)
                                    IMAX = I
         END IF
   80 CONTINUE
      FPMAX = FREQ(IMAX)
      FNF = FLOAT( NF )
      PAVE = PSUM / FNF
      RATIO = PMAX / PAVE
      EXPON = 1.0 / (FNF - 1.0)
      FISHER = FNF*( 1.0 - ( 0.05/FNF )**EXPON )
      IF ( RATIO .GT. FISHER ) THEN
                                   PMSIG = .TRUE.
              ELSE
                  PMSIG = .FALSE.
      END IF
C
C  Compute cumulative periodogram and endpoint coordinates of upper
C          and lower 95% bound lines for white noise confidence band
C
      C0 = 1.36 / SQRT( FLOAT(N) / 2.0 )
C
      SUMP = PSUM
      PSUM = 0.0
      NOUT = 0
      CUMPER(1) = 0.0
      SLOPE = 2.0 * DELT
      DO 100 J=1,NF
         PSUM = PSUM + PERDG(J)
         CUMPER(J) = PSUM / SUMP
         IF ( ( CUMPER(J) .LT. -C0 + SLOPE*FREQ(J) ) .OR.
     1        ( CUMPER(J) .GT.  C0 + SLOPE*FREQ(J) ) )  NOUT = NOUT + 1
  100 CONTINUE
C
      IF ( FLOAT( NOUT ) / FLOAT( NF ) .LE. 0.05 ) THEN
                                                   WHN = .TRUE.
                     ELSE
                     WHN = .FALSE.
      END IF
C
      FUB(1) = 0.0
      CUB(1) = C0
      FUB(2) = (1.0 - C0) / (2.0*DELT)
      CUB(2) = 1.0
      FLB(1) = C0 / (2.0*DELT)
      CLB(1) = 0.0
      FLB(2) = 1.0 / (2.0*DELT)
      CLB(2) = 1.0 - C0
C
      RETURN
C
      END

