c
c
C++++ Subroutine for complex SPECTRAL ESTIMATION for complex data  ++++
C 
c    --------------------------------------------------------------------
      subroutine CMPLXPER ( T, X, NP, DELT, NPSD, FREQ, RPSD, AMP,
c
     1                      WINDOW, dbPSD, YLABEL, TS,
c
     1                      W, CWRK1, CWRK2, PSD,
c
     2                      R, ARA, BMA, CWRK3, CWRK4, IOU,
c
     3                      TSIDENT )
c    --------------------------------------------------------------------
c
c         T = array of times (length = NP)
c         X = complex time-series array (length = NP)
c        NP = no. terms in the time series
c      DELT = time spacing for time-series terms
c      NPSD = no. frequency points desired in the periodogram
c      FREQ = array of desired frequencies (output)
c      RPSD = calculated periodogram estimate (output)
c       AMP = calculated amplitude spectrum (output)
c    WINDOW = real working vector for windowing (length = NPSD) 
c     dbPSD = real working vector for calculating periodogram in dB units 
c    YLABEL = character string label for the time-series values
c        TS = logical variable indicating the data are a time series
c         W = \
c     CWRK1 =  \ complex working vectors
c     CWRK2 =  /    (of length NPSD)
c       PSD = /
c         R = \
c       ARA =  \
c       BMA =   > complex working vectors
c     CWRK3 =  /     (of length NP)
c     CWRK4 = /
c       IOU = hardcopy output unit
c   TSIDENT = character string identifying the input time series (9 characters)
C
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      LOGICAL  TAPER, TS
C
      CHARACTER*3  IANS3
      CHARACTER*(*)  TSIDENT
      CHARACTER*(*) YLABEL
C
      REAL  T(NP), FREQ(NPSD), RPSD(NPSD), AMP(NPSD), WINDOW(NPSD),
     1      dbPSD(NPSD)
C
      COMPLEX  X(NP), W(NPSD), CWRK1(NPSD), CWRK2(NPSD),
     1         PSD(NPSD)
C
      COMPLEX  R0, R(NP), ARA(NP), BMA(NP), CWRK3(NP), CWRK4(NP)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C
      WRITE (IOU,10)
   10 FORMAT ( ' '/' ' )
C
C
      TAPER = .FALSE.
C
C--------------
  999 CONTINUE
C--------------     THE COMMAND POINT
C
      PRINT *, ' TYPE  ACS  FOR AUTOCORRELATION ESTIMATES'
      PRINT *, '       PER  FOR CLASSICAL PERIODOGRAM'
      PRINT *, '       COR  FOR CLASSICAL CORRELOGRAM'
      PRINT *, '       YUW  FOR YULE-WALKER EST. OF AR PARAMS. AND PSD'
      PRINT *, '       BRG  FOR BURG EST. OF AR PARAMS. AND PSD'
      PRINT *, '       ARM  FOR ARMA PSD  '
      PRINT *, '       COV  FOR COVARIANCE EST. AND PSD'
      PRINT *, '       MCV  FOR MODIFIED COVARIANCE EST. AND PSD'
      PRINT *, '       RTN  TO RETURN TO CALLING PROGRAM  '
C
      READ *, IANS3
C                                                       --------
      IF ( (IANS3 .EQ. 'RTN') .OR. (IANS3 .EQ. 'rtn') )  RETURN
C                                                       --------
C
      IF ( (IANS3 .EQ. 'PER') .OR. (IANS3 .EQ. 'per') ) GO TO 50
C
      IF ( (IANS3 .EQ. 'ACS') .OR. (IANS3 .EQ. 'acs') ) GO TO 100
C
      IF ( (IANS3 .EQ. 'COR') .OR. (IANS3 .EQ. 'cor') ) GO TO 150
C
      IF ( (IANS3 .EQ. 'YUW') .OR. (IANS3 .EQ. 'yuw') ) GO TO 200
C
      IF ( (IANS3. EQ. 'BRG') .OR. (IANS3 .EQ. 'brg') ) GO TO 250
C
      IF ( (IANS3 .EQ. 'ARM') .OR. (IANS3 .EQ. 'arm') ) GO TO 300
C
      IF ( (IANS3 .EQ. 'COV') .OR. (IANS3 .EQ. 'cov') ) GO TO 350
C
      IF ( (IANS3 .EQ. 'MCV') .OR. (IANS3 .EQ. 'mcv') ) GO TO 400
C
      GO TO 999
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
   50 CONTINUE
C                 COMPUTE PERIODOGRAM
C
      NSHIFT =  1
      NSAMP  = NP
C
C   ---------------------------------------------------------------
      CALL PERIODOGRAM ( NP, NSHIFT, NSAMP, DELT, X, X, NSEG, PSD,
C
     1                   NPSD, W, CWRK1, CWRK2, WINDOW, TAPER )
C   ---------------------------------------------------------------
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
   99 CONTINUE
c
c--- Complete calculation of Periodogram and Plot it
C
C
      NPSDO2 = NPSD / 2
      DELF = 1.0 / ( NPSD*DELT )
C
      DO 51 K=1,NPSD
         FREQ(K) = (K - 1 - NPSDO2) * DELF
   51 CONTINUE
C
      DO 53 K=1,NPSDO2
         RPSD(K) = CABS( PSD(NPSDO2+K) )
   53 CONTINUE
C
      DO 55 K=(NPSDO2+1),NPSD
         RPSD(K) = CABS( PSD(K-NPSDO2) )
   55 CONTINUE
C
c----- Compute amplitude spectrum
c
      DO K = 1, NPSD
         AMP(K) = 2.0 * SQRT( RPSD(K) / NP )
      END DO
C
C    -----------------------------------------------------------
      CALL SCPLT1 ( FREQ, RPSD, NPSD, 0, 1, 
c
     2             'Complex '//TSIDENT//' PSD',
C
     1             'Frequency [kHz]', 'Variance  ( Power )' )
C   ------------------------------------------------------------
C
C---  PLOT ESTIMATED PSD IN DECIBELS ---
C
      PMAX = RPSD(1)
      DO 90 K=1,NPSD
         IF ( RPSD(K) .LT. 0.0 ) RPSD(K) = -RPSD(K)
         IF ( RPSD(K) .GT. PMAX ) PMAX = RPSD(K)
   90 CONTINUE
      WRITE (IOU,10)
      WRITE (IOU,*) ' PMAX = ', PMAX
      DO 95 K=1,NPSD
         dbPSD(K) = 10.0*ALOG10( RPSD(K) / PMAX )
   95 CONTINUE
C
C    -------------------------------------------------------------
      CALL SCPLT1 ( FREQ, dbPSD, NPSD, 0, 1, 
c
     1              'Complex '//TSIDENT//' PSD',
C
     1              'Frequency [kHz]', 'Rel. Variance [dB]' )
C    -----------------------------------------------------------
c
      DO K = 1, NPSD
         WINDOW(K) = 0.0
      END DO
c
c   -------------------------------------------------------------------
      CALL TWOPLT ( FREQ, RPSD, NPSD, FREQ, WINDOW, NPSD, 0, 1,
C
     1              'Complex '//TSIDENT//' PSD',
C
     2              'Frequency [kHz]', 'Power (Variance)',
C
     3              FREQ, dbPSD, NPSD, FREQ, WINDOW, NPSD, 0, 1,
C
     4              'Complex '//TSIDENT//' PSD',
c
     5              'Frequency [kHz]', 'Rel. Variance [dB]', TS )
c   ---------------------------------------------------------------------
c
c    --------------------------------------------------------------------
      CALL TWOPLT ( FREQ, dbPSD, NPSD, FREQ, WINDOW, NPSD, 0, 1,
C
     1              'Complex '//TSIDENT//' PSD',
C
     2              'Frequency [kHz]', 'Rel. Power [dB]',
C
     3              FREQ, AMP, NPSD, FREQ, WINDOW, NPSD, 0, 1,
C
     4              'Complex '//TSIDENT//' Amplitude Spectrum',
c
     5              'Frequency [kHz]', YLABEL, TS )
C    --------------------------------------------------------------------
c
C
      GO TO 999
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
  100 CONTINUE
C               COMPUTE AUTOCORRELATION SERIES
C
      MLAG = NP / 2
      MODE = 0
C
C    --------------------------------------------------
      CALL CORRELATION ( NP, MLAG, MODE, X, X, R0, R )
C    --------------------------------------------------
C
      WRITE (IOU,*) ' R[ 0] =', R0
      DO 120 I=1,MLAG
         WRITE (IOU,*) ' R[', I, '] =', R(I)
  120 CONTINUE
C
      GO TO 999
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
  150 CONTINUE
C              COMPUTE CORRELOGRAM
C
      MLAG = NP / 2
C   ---------------------------------------------------------------
      CALL CORRELOGRAMPSD ( NP, MLAG, DELT, X, X, PSD, NPSD, W,
C
     1                      R, WINDOW, 0, TAPER )
C   ---------------------------------------------------------------
C
      GO TO 99
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
  200 CONTINUE
C              COMPUTE YULE-WALKER ESTIMATE OF AUTOREGRESSIVE PARAMETERS
C
      IP = NP / 2
      LBIAS = 1
C
C    ----------------------------------------------------------
      CALL YULEWALKER ( NP, IP, LBIAS, X, RHO, ARA, ISTAT, R )
C    ----------------------------------------------------------
C
      WRITE (IOU,10)
      WRITE (IOU,10)
      WRITE (IOU,*) ' YULEWALKER EST. OF AUTOREGRESSIVE PARAMETERS'
      WRITE (IOU,*)
      DO 205 I=1,IP
          WRITE(IOU,*) ' ARA(', I, ') = ', ARA(I)
  205 CONTINUE
C
      IF ( ISTAT .NE. 0 ) THEN
                          WRITE (IOU,210)
  210    FORMAT (' '/' '/' CORRELATION MATRIX WAS SINGULAR  ')
                          PRINT 210
              GO TO 999
      END IF
C
      IQ = 0
C
      GO TO 310
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
  250 CONTINUE
C
      IP = NP / 2
C
C    --------------------------------------------------------
      CALL BURG ( NP, IP, X, RHO, ARA, ISTAT, CWRK3, CWRK4 )
C    --------------------------------------------------------
C
      WRITE (IOU,10)
      WRITE (IOU,10)
      WRITE (IOU,*) ' BURG EST. OF AUTOREGRESSIVE PARAMETERS'
      WRITE (IOU,*)
      DO 255 I=1,IP
          WRITE(IOU,*) ' ARA(', I, ') = ', ARA(I)
  255 CONTINUE
C
      IF ( ISTAT .NE. 0 ) THEN
                          WRITE (IOU,260)
  260    FORMAT (' '/' '/' ILL-CONDITIONED CORRELATION MATRIX',
     1              /' '/' DRIVING NOISE VARIANCE RHO IS NEGATIVE ')
                          PRINT 260
              GO TO 999
      END IF
C
      IQ = 0
C
      GO TO 310
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
  300 CONTINUE
C              COMPUTE  ARMA PSD SPECTRUM
C
      IP = 1
      IQ = 1
      ARA(1) = CMPLX( 0.8 , 0.9 )
      BMA(1) = CMPLX( 0.1 , - 0.3 )
      RHO = 2.0
C
  310 CONTINUE
C
C   -----------------------------------------------------------
      CALL ARMAPSD ( IP, IQ, RHO, DELT, ARA, BMA, RPSD, ISTAT,
C
     1               NPSD, W, CWRK1, CWRK2 )
C   -----------------------------------------------------------
C
      DO 320 I=1,NPSD
          PSD(I) = CMPLX( RPSD(I), 0.0 )
  320 CONTINUE
C
      GO TO 99
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
  350 CONTINUE
C             SOLVE COVARIANCE EQUATIONS FOR LINEAR PREDICTION COEFFICIENTS
C             AND COMPUTE PSD ESTIMATE FOR FORWARD PREDICTION COEFFICIENTS.
C
      IP = NP / 2
C
C   -------------------------------------------------------------
      CALL COVAR ( NP, IP, X, RHOF, ARA, RHOB, BMA, ISTAT,
C
     1                                         CWRK3, CWRK4, R )
C   -------------------------------------------------------------
C
      WRITE (IOU,10)
      WRITE (IOU,*) ' COVARIANCE ESTIMATES OF FORWARD PREDICTION FILTER'
      WRITE (IOU,*) ' '
      WRITE (IOU,*) ' RHOF =', RHOF
      WRITE (IOU,*) ' '
      DO 365 I=1,IP
         WRITE (IOU,*) ' AF(', I, ') =', ARA(I)
  365 CONTINUE
      IF ( ISTAT .EQ. 0 ) THEN
                          IQ = 0
                          RHO = RHOF
                          GO TO 310
              ELSE
              WRITE (IOU,*)
              WRITE (IOU,*) ' *********** ISTAT =', ISTAT
              GO TO 999
      END IF
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
  400 CONTINUE
C             SOLVE MODIFIED COVARIANCE EQNS FOR LIN. PREDICTION COEFFS.
C             AND COMPUTE CORRESPONDING PSD ESTIMATE
C
      IP = NP / 2
C
C    ---------------------------------------------------------------
      CALL MODCOVAR ( NP, IP, X, RHO, ARA, ISTAT, CWRK3, CWRK4, R )
C    ---------------------------------------------------------------
C
      WRITE (IOU,10)
      WRITE (IOU,*) ' MODIFIED COVARIANCE EST. OF PREDICTION FILTER'
      WRITE (IOU,*) ' '
      WRITE (IOU,*) ' RHO =', RHO
      WRITE (IOU,*) ' '
      DO 415 I=1,IP
         WRITE (IOU,*) ' A(', I, ') =', ARA(I)
  415 CONTINUE
      IF ( ISTAT .EQ. 0 ) THEN
                          IQ = 0
                          GO TO 310
              ELSE
              WRITE (IOU,*)
              WRITE (IOU,*) ' *********** ISTAT =', ISTAT
              GO TO 999
      END IF
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      END



