C
C-------------------------------------------------------------------
C MAIN PROGRAM: MAIN LINE PROGRAM TO COMPUTE THE REAL CEPSTRUM
C               AND THE MINIMUM-PHASE RECONSTRUCTION OF AN
C               ARBITRARY-PHASE, REAL SEQUENCE X(N)
C AUTHORS:      THOMAS F. QUATIERI
C               M.I.T., CAMBRIDGE MASS. 02139
C               JOSE M. TRIBOLET
C               INSTITUTO SUPERIOR TECNICO, LISBON, PORTUGAL
C-------------------------------------------------------------------
C
      DIMENSION X(512), CX(1026)
      REAL ICX(1026)
C
C DIMENSION REQUIREMENTS:
C
C (SEE SUBROUTINES RCEPS AND ICEPS)
C
C
C DESCRIPTION OF ARRAYS:
C
C   X = ARRAY CONTAINING THE SEQUENCE X(N)
C  CX = ARRAY CONTAINING THE ZERO-PHASE(REAL)
C       OR MINIMUM-PHASE CEPSTRUM OF X(N)
C ICX = ARRAY CONTAINING THE INVERSE CEPSTRUM
C       (MINIMUM-PHASE RECONSTRUCTION OF X(N))
C
C
C DESCRIPTION OF VARIABLES:
C
C    NX = LENGTH OF THE SEQUENCE X(N)
C NFFTF = LENGTH OF THE FFTS IN FORWARD CEPSTRUM
C NFFTI = LENGTH OF THE FFTS IN INVERSE CEPSTRUM
C         (INVERSE LENGTH NEED NOT EQUAL FORWARD LENGTH)
C
C
C SUBROUTINES CALLED:
C
C RCEPS = SUBROUTINE TO COMPUTE THE REAL CEPSTRUM OF A
C         REAL SEQUENCE
C
C ICCEPS = SUBROUTINE TO COMPUTE THE INVERSE COMPLEX CEPSTRUM
C          (NOTE THAT "COMPLEX" REFERS TO AN ARBITRARY,BUT
C          REAL-VALUED SEQUENCE)
C
C
C SET FFT LENGTHS AND DEFINE OUTPUT DEVICE IODEV
C
      NFFTF = 1024
      NFFTI = 1024
      IODEV = I1MACH(2)
C
C GENERATE THE TEST SIGNAL WITH Z-TRANSFORM:
C
C X(Z)=1.0+(1.0/.99)*Z**(-1) ---ONE ZERO AT 1.01010101
C
      NX = 2
      X(1) = 1.0
      X(2) = 1.01010101
C
C PAD X(N) WITH ZEROS:STORE IN CX
C
      DO 10 I=1,NX
        CX(I) = X(I)
  10  CONTINUE
      INITL = NX + 1
      IEND = NFFTF + 2
      DO 20 I=INITL,IEND
        CX(I) = 0.0
  20  CONTINUE
C
C COMPUTE THE REAL CEPSTRUM
C
      CALL RCEPS(NFFTF, CX)
C
C WRITE LAST AND FIRST 32 VALUES OF THE REAL CEPSTRUM
C
      WRITE (IODEV,9999)
9999  FORMAT (1X, 10X, 35H*** RCPCIP- TEST PROGRAM OUTPUT ***//)
      INITL = NFFTF - 31
      WRITE (IODEV,9998)
9998  FORMAT (1X, 15H REAL CEPSTRUM //)
      WRITE (IODEV,9997) (I,CX(I),I=INITL,NFFTF)
      WRITE (IODEV,9997) (I,CX(I),I=1,32)
9997  FORMAT (1X/4(5H  CX(, I4, 2H)=, F7.3))
C
C FIND THE MINIMUM-PHASE COUNTERPART OF THE ZERO-PHASE(REAL)
C CEPSTRUM AND PAD WITH ZEROS:STORE IN ICX (NOTE THAT THE FORWARD
C FFT LENGTH EQUALS THE INVERSE FFT LENGTH IN THIS EXAMPLE.
C IN GENERAL,THE FOLLOWING PROCEDURE SHOULD BE TAILORED TO THE
C PARTICULAR FFT LENGTHS).
C
      ICX(1) = CX(1)
      IEND = NFFTF/2 + 1
      DO 30 I=2,IEND
        ICX(I) = 2.0*CX(I)
  30  CONTINUE
      INITL = IEND + 1
      IEND = NFFTI + 2
      DO 40 I=INITL,IEND
        ICX(I) = 0.0
  40  CONTINUE
C
C COMPUTE THE INVERSE CEPSTRUM
C
      CALL ICCEPS(NFFTI, ICX)
C
C WRITE LAST AND FIRST 32 VALUES OF THE MINIMUM-PHASE RECONSTRUCTION
C
      INITL = NFFTI - 31
      WRITE (IODEV,9996)
9996  FORMAT (///29H MINIMUM-PHASE RECONSTRUCTION//)
      WRITE (IODEV,9995) (I,ICX(I),I=INITL,NFFTI)
      WRITE (IODEV,9995) (I,ICX(I),I=1,32)
9995  FORMAT (1X/4(5H ICX(, I4, 2H)=, F7.3))
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: RCEPS
C COMPUTE THE REAL CEPSTRUM OF A REAL SEQUENCE X(N)
C-----------------------------------------------------------------------
C
      SUBROUTINE RCEPS(NFFT, CX)
C
C DESCRIPTION OF ARGUMENTS:
C
C NFFT = LENGTH OF THE FFT
C   CX = ARRAY CONTAINING THE REAL CEPSTRUM
C        (ALSO USED AS AN AUXILIARY ARRAY)
C
      DIMENSION CX(1)
C
C DIMENSION REQUIREMENTS:
C
C DIM(CX) .GE. NFFT+2
C
C
C SUBROUTINES CALLED:
C
C FFA,FFS-SUBROUTINES TO COMPUTE THE FFT AND IFFT (RADIX 2)
C
C FUNCTIONS CALLED:
C
C AMODSQ = FUNCTION TO COMPUTE THE MODULUS SQUARED
C          OF A COMPLEX NUMBER
C
C
C
C
C TRANSFORM X(N):FFT
C
      CALL FFA(CX, NFFT)
C
C COMPUTE THE LOGMAGNITUDE OF THE SPECTRUM:STORE IN THE
C ODD-INDEXED VALUES OF CX;CLEAR THE EVEN-INDEXED VALUES OF CX
C
      IEND = NFFT/2 + 1
      IO = -1
      DO 10 I=1,IEND
        IO = IO + 2
        IE = IO + 1
        AMAGSQ = AMODSQ(CX(IO),CX(IE))
        CX(IO) = .5*ALOG(AMAGSQ)
        CX(IE) = 0.0
  10  CONTINUE
C
C COMPUTE THE REAL CEPSTRUM:IFFT
C
      CALL FFS(CX, NFFT)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ICCEPS
C COMPUTE THE INVERSE COMPLEX CEPSTRUM (NOTE THAT "COMPLEX" REFERS
C TO AN ARBITRARY, BUT REAL-VALUED CEPSTRAL SEQUENCE)
C-----------------------------------------------------------------------
C
      SUBROUTINE ICCEPS(NFFT, ICX)
C
C DESCRIPTION OF ARGUMENTS:
C
C NFFT = LENGTH OF THE FFT
C  ICX = ARRAY CONTAINING THE INVERSE CEPSTRUM
C        (ALSO USED AS AN AUXILIARY ARRAY)
C
      REAL ICX(1)
C
C DIMENSION REQUIREMENTS:
C
C DIM(ICX) .GE. NFFT+2
C
C
C SUBROUTINES CALLED:
C
C FFA,FFS-SUBROUTINES TO COMPUTE THE FFT AND IFFT (RADIX 2)
C
C
C
C TRANSFORM CX(N):FFT
C
      CALL FFA(ICX, NFFT)
C
C EXPONENTIATE:STORE THE REAL AND IMAGINARY VALUES IN THE
C ODD-INDEXED AND EVEN-INDEXED LOCATIONS OF ICX, RESPECTIVELY
C
      IEND = NFFT/2 + 1
      IO = -1
      DO 10 I=1,IEND
        IO = IO + 2
        IE = IO + 1
        RMAG = EXP(ICX(IO))
        ICX(IO) = RMAG*COS(ICX(IE))
        ICX(IE) = RMAG*SIN(ICX(IE))
  10  CONTINUE
C
C INVERSE TRANSFORM:IFFT
C
      CALL FFS(ICX, NFFT)
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: AMODSQ
C COMPUTE THE SQUARE OF THE MODULUS OF A COMPLEX NUMBER
C-----------------------------------------------------------------------
C
      FUNCTION AMODSQ(ZR, ZI)
C
C DESCRIPTION OF ARGUMENTS:
C
C ZR = REAL PART OF THE COMPLEX NUMBER
C ZI = IMAGINARY PART OF THE COMPLEX NUMBER
C
      AMODSQ = SNGL(DBLE(ZR)*DBLE(ZR)+DBLE(ZI)*DBLE(ZI))
      RETURN
      END