C
C
      SUBROUTINE ARMAPSD ( IP, IQ, RHO, T, A, B, PSD, ISTAT,
C
     1                     NPSD, W, DEN, NUM )
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C   This subroutine computes a set of NPSD power spectral density values
C   across the frequency band from -1/2T to 1/2T Hz, where T is the sam-
C   ple interval in sec. ,  given the complex ARMA parameters of an ARMA
C   model.  The   FFT  is used to efficiently evaluate the numerator and
C   denominator polynomials of the ARMA PSD function.    For the case of
C   real  ARMA  parameters,  the spectral values are symmetric about the
C   element indexed as NPSD/2 + 1 .
C
C   The PSD for just an AR model may be obtained as a special case  of
C   this subroutine by setting IQ to zero. Similarly, the PSD for just
C   a MA model may be obtained as a special case by setting IP to zero.
C
C   Input parameters:
C
C          IP    - AR model order (integer)
C          IQ    - MA model order (integer)
C          RHO   - White noise variance of model input (real)
C          T     - Sample interval in seconds (real)
C          A     - Complex array of autoregressive parameters A(1) to A(IP)
C          B     - Complex array of moving average parameters B(1) to B(IQ)
C          NPSD  - No. Frequencies in the PSD Estimate
C          W     - Complex working vector ( length .ge. NPSD )
C          DEN   - Complex working vector ( length .ge. NPSD )
C          NUM   - Complex working vector ( length .ge. NPSD )
C
C   Output parameters:
C
C          PSD   - Real array of power spectral density values
C          ISTAT - Integer status indicator at time of exit:
C                    0 for normal exit
C                    1 if IP or IQ < 0
C
C   Notes:
C
C     External array A must be dimensioned .GE. IP, array B must be
C     dimensioned .GE. IQ, and array PSD must be dimensioned .GE.
C     NPSD in the calling program.  Subroutines PREFFT and FFT (in
C     Appendix 2.A) are required.
C
C-----      PARAMETER NPSD=4096         ! Must be a power of 2
C-----PARAMETER  ( NPSD = 4096 )
C                        ! Must be a power of 2
C
      COMPLEX A(1),B(1),W(NPSD),DEN(NPSD),NUM(NPSD)
      REAL PSD(1),RHO,T
      ISTAT=0
C   Determine if error condition exists with order values
      IF (IP .GE. 0 .AND. IQ .GE. 0)  GO TO 5
      ISTAT=1
      RETURN
C   Set up exponential table for the FFT (no zero padding)
5     CALL PREFFT (NPSD,0,NEXP,W)
      IF (IP .EQ. 0)  GO TO 30
      DEN(1)=(1.,0.)
      DO 10 K=1,IP
10      DEN(K+1)=A(K)
      DO 20 K=IP+2,NPSD
20      DEN(K)=(0.,0.)
      CALL FFT (NPSD,0,1.,NEXP,W,DEN)
30    IF (IQ .EQ. 0)  GO TO 60
      NUM(1)=(1.,0.)
      DO 40 K=1,IQ
40      NUM(K+1)=B(K)
      DO 50 K=IQ+2,NPSD
50      NUM(K)=(0.,0.)
      CALL FFT (NPSD,0,1.,NEXP,W,NUM)
60    IF (IP .EQ. 0)  GO TO 80
      IF (IQ .EQ. 0)  GO TO 100
C  Equation (6.8)
      DO 70 K=1,NPSD
70      PSD(K)=RHO*T*(REAL(NUM(K))**2+AIMAG(NUM(K))**2)/
     *         (REAL(DEN(K))**2+AIMAG(DEN(K))**2)
      RETURN
C  Equation (6.12)
80    DO 90 K=1,NPSD
90      PSD(K)=RHO*T*(REAL(NUM(K))**2+AIMAG(NUM(K))**2)
      RETURN
C  Equation (6.14)
100   DO 110 K=1,NPSD
110     PSD(K)=RHO*T/(REAL(DEN(K))**2+AIMAG(DEN(K))**2)
      RETURN
      END
