Make your own free website on Tripod.com
Computation of the Complex Cepstrum

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: MAIN LINE PROGRAM TO COMPUTE THE COMPLEX CEPSTRUM
C               OF A REAL SEQUENCE X(N)
C AUTHORS:      JOSE M. TRIBOLET
C               INSTITUTO SUPERIOR TECNICO, LISBON, PORTUGAL
C               THOMAS F. QUATIERI
C               M.I.T., CAMBRIDGE, MASS. 02139
C-----------------------------------------------------------------------
C
      DIMENSION X(1024), CX(1026), AUX(1026), SPECM(12), SPECP(12)
      COMMON PI, TWOPI, THLINC, THLCON, NFFT, NPTS, N, L, H, H1, DVTMN2
      LOGICAL ISSUC
C
C DIMENSION REQUIREMENTS:
C
C (SEE SUBROUTINE CCEPS)
C
C DESCRIPTION OF ARRAYS:
C
C     X = ARRAY CONTAINING THE SEQUENCE X(N)
C    CX = ARRAY CONTAINING THE COMPLEX CEPSTRUM OF X(N)
C   AUX = AUXILIARY ARRAY
C SPECM = MAGNITUDES OF SPECTRAL ZEROS (SEE COEFF)
C SPECP = PHASE OF SPECTRAL ZEROS (SEE COEFF)
C
C DESCRIPTION OF VARIABLES:
C
C  IOUTD = OUTPUT DEVICE NUMBER
C     NX = LENGTH OF THE SEQUENCE X(N)
C   NFFT = LENGTH OF THE FFT
C THLINC = PHASE INCREMENT THRESHOLD(USED TO OBTAIN MORE
C          CONFIDENT PHASE ESTIMATE NEAR SHARP ZEROS)
C THLCON = PHASE CONSISTENCY THRESHOLD
C  ISSUC = .TRUE. IF COMPLEX CEPSTRUM SUCCESSFULLY COMPUTED;
C          .FALSE. OTHERWISE(PHASE ESTIMATION INCOMPLETE)
C
C SUBROUTINES CALLED:
C
C COEFF = SUBROUTINE TO COMPUTE A SEQUENCE X(N),BY
C         PRESCRIBING ITS SPECTRAL ZEROS
C CCEPS = SUBROUTINE TO COMPUTE THE COMPLEX CEPSTRUM
C         OF A REAL SEQUENCE X(N)
C
C
C
C SET LENGTH OF FFT,THRESHOLD LEVELS,AND OTHER CONSTANTS
C
      IOUTD = I1MACH(2)
      NFFT = 1024
      THLINC = 1.5
      THLCON = .5
      PI = 4.*ATAN(1.)
      TWOPI = 2.*PI
C
C GENERATE THE TEST SIGNAL WITH ZEROS AT:
C
C |0.9|EXP(+/-JPI/4)
C |1.1|EXP(+/-J(PI/4+PI/8192))
C |0.9|EXP(+/-J(PI/4+2*PI/8192))
C
C ***************************************
C *                                     *
C *       P L E A S E    N O T E        *
C *                                     *
C *    USER IS STRONGLY ENCOURAGED TO   *
C *    RUN TEST PROGRAM WITH A VARIETY  *
C *    OF TEST SIGNALS TO BECOME FULLY  *
C *    ACQUAINTED WITH THIS PROGRAM'S   *
C *    CAPABILITIES.PHASE UNWRAPPING    *
C *    MAY OFTEN BE NOT POSSIBLE DUE    *
C *    EITHER TO THEORETICAL REASONS    *
C *    (ZEROS ON THE UNIT CIRCLE )      *
C *    OR TO COMPUTATIONAL REASONS      *
C *    (ZEROS TOO CLOSE TO UNIT CIRCLE) *
C *                                     *
C ***************************************
C
      NX = 6
      SPECM(1) = .9
      SPECP(1) = (PI/4.)
      SPECM(2) = 1.1
      SPECP(2) = (PI/4.+PI/8192.)
      SPECM(3) = .9
      SPECP(3) = (PI/4.+TWOPI/8192.)
      CALL COEFF(NX, X, SPECM, SPECP)
C
C WRITE INPUT DATA SPECIFICATIONS
C
      WRITE (IOUTD,9999)
9999  FORMAT (1X, 10X, 35H*** CCMAIN- TEST PROGRAM OUTPUT ***//)
      WRITE (IOUTD,9998) NX
9998  FORMAT (1X, 34HINPUT SIGNAL-LENGTH (IN SAMPLES) =, I4/)
      I = NX/4
      J = 4*I
      IF (NX.EQ.J) GO TO 20
      J = J + 4
      K = NX + 1
      DO 10 I=K,J
        X(I) = 0.
  10  CONTINUE
  20  WRITE (IOUTD,9997) (I,X(I),I=1,J)
9997  FORMAT (1X/4(4H  X(, I4, 2H)=, F8.4))
C
C COMPUTE THE COMPLEX CEPSTRUM
C
      CALL CCEPS(NX, X, ISNX, ISFX, ISSUC, CX, AUX)
C
C CHECK WHETHER PHASE UNWRAPPING SUCCESSFUL; IF SO
C WRITE OUT DATA;OTHERWISE STOP|
C
      IF (ISSUC) GO TO 30
      WRITE (IOUTD,9996)
9996  FORMAT (1X///1X, 23HPHASE ESTIMATION FAILED)
      STOP
C
C WRITE SIGN,LINEAR PHASE,AND LAST AND FIRST 32 VALUES OF
C THE COMPLEX CEPSTRUM
C
  30  WRITE (IOUTD,9995) ISNX, ISFX
9995  FORMAT (1X/1X, 5HSIGN=, I2//1X, 13HLINEAR PHASE=, I4//)
      INITL = NFFT - 31
      WRITE (IOUTD,9994)
9994  FORMAT (1X, 17HCOMPLEX CEPSTRUM //)
      WRITE (IOUTD,9993) (I,CX(I),I=INITL,NFFT)
      WRITE (IOUTD,9993) (I,CX(I),I=1,32)
9993  FORMAT (1X/4(4H CX(, I4, 2H)=, F8.4))
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COEFF
C COMPUTES A SEQUENCE C(N) WITH NC SAMPLES BY PRESCRIBING
C THE MAGNITUDES AND PHASES OF ITS SPECTRAL ZEROS.
C-----------------------------------------------------------------------
C
      SUBROUTINE COEFF(NC, C, SM, SP)
      DOUBLE PRECISION S(2,31), Z1, Z2, Z3, Z4
      REAL M(12), P(12), C(12), SM(1), SP(1)
C
C DESCRIPTION OF ARGUMENTS
C
C NC   -ON INPUT EQUALS THE TOTAL NUMBER OF ZEROS OF C(Z)
C      -ON OUTPUT EQUALS THE SEQUENCE LENGTH
C C    -DSEQUENCE WITH PRESCRIBED Z-TRANSFORM
C M    -ARRAY OF SPECTRAL MAGNITUDES.ONLY ONE ENTRY
C       PER EACH COMPLEX CONJUGATE ZERO PAIR MUST BE SPECIFIED
C P    -ARRAY OF SPECTRAL PHASES IN RADIANS.
C ONLY ONE ENTRY PER EACH COMPLEX CONJUGATE ZERO PAIR MUST
C BE SPECIFIED.
C
      N = NC
      I = 1
      IND = 1
  10  CONTINUE
      M(I) = SM(IND)
      P(I) = SP(IND)
      IND = IND + 1
      IF (P(I).EQ.0.) GO TO 20
      M(I+1) = M(I)
      P(I+1) = -P(I)
      I = I + 1
  20  I = I + 1
      IF (I.LT.N) GO TO 10
      Y = P(1)
      S(1,1) = -DBLE(M(1)*COS(Y))
      S(2,1) = -DBLE(M(1)*SIN(Y))
      S(1,2) = 1.
      S(2,2) = 0.
      IF (N.EQ.1) GO TO 50
      DO 40 J=2,N
        S(1,J+1) = S(1,J)
        S(2,J+1) = S(2,J)
        Y = P(J)
        Z3 = -DBLE(M(J)*COS(Y))
        Z4 = -DBLE(M(J)*SIN(Y))
        DO 30 K=1,J
          M0 = J - K
          M1 = M0 + 1
          Z1 = S(1,M1)
          Z2 = S(2,M1)
          IF (M0.EQ.0) GO TO 30
          S(1,M1) = S(1,M0) + Z1*Z3 - Z2*Z4
          S(2,M1) = S(2,M0) + Z1*Z4 + Z2*Z3
  30    CONTINUE
        S(1,1) = Z1*Z3 - Z2*Z4
        S(2,1) = Z1*Z4 + Z2*Z3
  40  CONTINUE
  50  CONTINUE
      NC = N + 1
      DO 60 I=1,NC
        J = NC - I + 1
        C(I) = S(1,J)
  60  CONTINUE
      C(1) = 1.
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CCEPS
C SUBROUTINE TO COMPUTE THE COMPLEX CEPSTRUM OF A SEQUENCE X(N)
C-----------------------------------------------------------------------
C
      SUBROUTINE CCEPS(NX, X, ISNX, ISFX, ISSUC, CX, AUX)
C
C DESCRIPTION OF ARGUMENTS:
C
C    NX = LENGTH OF THE SEQUENCE X(N)
C     X = ARRAY CONTAINING THE SEQUENCE X(N)
C  ISNX = INTEGER VALUE,EITHER +1 OR  -1 DEPENDING
C         ON THE SIGN REVERSAL OF X(N)
C  ISFX = INTEGER VALUE INDICATING THE AMOUNT OF SHIFT
C         ON X(N) DUE TO LINEAR PHASE REMOVAL
C ISSUC = .TRUE. IF PHASE ESTIMATION COMPLETE;
C         .FALSE. OTHERWISE
C    CX = ARRAY CONTAINING THE COMPLEX CEPSTRUM
C   AUX = AUXILIARY ARRAY
C
      DIMENSION X(1), CX(1), AUX(1)
C
C DIMENSION REQUIREMENTS:
C
C NX       .LE. NFFT
C DIM(X)   .LE. NFFT
C DIM(CX)  .GE. NFFT+2
C DIM(AUX) .GE. NFFT+2
C
      COMMON PI, TWOPI, THLINC, THLCON, NFFT, NPTS, N, L, H, H1, DVTMN2
      LOGICAL ISSUC
C
C DESCRIPTION OF VARIABLES:
C
C THLINC = PHASE INCREMENT THRESHOLD(USED TO OBTAIN MORE
C          CONFIDENT ESTIMATE NEAR SHARP ZEROS)
C THLCON = PHASE CONSISTENCY THRESHOLD
C   NFFT = LENGTH OF FFT
C   NPTS = HALF THE LENGTH OF THE FFT
C DVTMN2 = TWICE THE MEAN OF THE PHASE DERIVATIVE
C
C SUBROUTINES CALLED:
C
C FFA,FFS-SUBROUTINES TO COMPUTE FFT AND IFFT (RADIX 2)
C
C FUNCTIONS CALLED:
C
C AMODSQ = FUNCTION TO COMPUTE THE MODULUS
C          SQUARED OF A COMPLEX NUMBER
C PHADVT = FUNCTION TO COMPUTE THE PHASE DERIVATIVE OF
C          A SPECTRAL VALUE
C PPVPHA = FUNCTION TO COMPUTE THE PRINCIPLE VALUE OF THE
C          PHASE OF A SPECTRAL VALUE
C PHAUNW = FUNCTION TO COMPUTE THE UNWRAPPED PHASE OF A
C          SPECTRAL VALUE
C
C INITIALIZATION
C
      NPTS = NFFT/2
      N = 12
      L = 2**N
      H = FLOAT(L)*FLOAT(NFFT)
      H1 = PI/H
      ISSUC = .TRUE.
      ISNX = 1
C
C TRANSFORM X(N) AND NX(N):FFT
C
      DO 10 I=1,NX
        CX(I) = X(I)
        AUX(I) = FLOAT(I-1)*X(I)
  10  CONTINUE
      INITL = NX + 1
      IEND = NFFT + 2
      DO 20 I=INITL,IEND
        CX(I) = 0.0
        AUX(I) = 0.0
  20  CONTINUE
C
C USE RADIX 2 FFT
C
      CALL FFA(CX, NFFT)
      CALL FFA(AUX, NFFT)
C
C CHECK IF SIGN REVERSAL IS REQUIRED
C
      IF (CX(1).LT.0.0) ISNX = -1
C
C COMPUTE MAGNITUDE OF SPECTRUM:STORE IN ODD-INDEXED
C VALUES OF AUX
C COMPUTE PHASE DERIVATIVE OF SPECTRUM:STORE IN EVEN-INDEXED
C VALUES OF AUX
C COMPUTE LINEAR PHASE ESTIMATE(MEAN OF THE PHASE DERIVATIVE):
C STORE TWICE THE ESTIMATE IN DVTMN2
C
      IO = -1
      DVTMN2 = 0.0
      IEND = NPTS + 1
      DO 30 I=1,IEND
        IO = IO + 2
        IE = IO + 1
        AMAGSQ = AMODSQ(CX(IO),CX(IE))
        PDVT = PHADVT(CX(IO),CX(IE),AUX(IO),AUX(IE),AMAGSQ)
        AUX(IO) = AMAGSQ
        AUX(IE) = PDVT
        DVTMN2 = DVTMN2 + PDVT
  30  CONTINUE
      DVTMN2 = (2.*DVTMN2-AUX(2)-PDVT)/FLOAT(NPTS)
C
C COMPUTE LOGMAGNITUDE:STORE IN ODD-INDEXED
C VALUES OF CX
C COMPUTE UNWRAPPED PHASE:STORE IN EVEN-INDEXED
C VALUES OF CX
C
      PPDVT = AUX(2)
      PPHASE = 0.0
      PPV = PPVPHA(CX(1),CX(2),ISNX)
      CX(1) = .5*ALOG(AUX(1))
      CX(2) = 0.0
      IO = 1
      DO 50 I=2,IEND
        IO = IO + 2
        IE = IO + 1
        PDVT = AUX(IE)
        PPV = PPVPHA(CX(IO),CX(IE),ISNX)
        PHASE = PHAUNW(X,NX,ISNX,I,PPHASE,PPDVT,PPV,PDVT,ISSUC)
C
C IF PHASE ESTIMATION SUCCESSFUL,CONTINUE;OTHERWISE RETURN
C
        IF (ISSUC) GO TO 40
        ISSUC = .FALSE.
        RETURN
  40    PPDVT = PDVT
        PPHASE = PHASE
        CX(IO) = .5*ALOG(AUX(IO))
        CX(IE) = PHASE
  50  CONTINUE
C
C REMOVE LINEAR PHASE COMPONENT
C
      ISFX = (ABS(PHASE/PI)+.1)
      IF (PHASE.LT.0.0) ISFX = -ISFX
      H = PHASE/FLOAT(NPTS)
      IE = 0
      DO 60 I=1,IEND
        IE = IE + 2
        CX(IE) = CX(IE) - H*FLOAT(I-1)
  60  CONTINUE
C
C COMPUTE THE COMPLEX CEPSTRUM:IFFT
C
      CALL FFS(CX, NFFT)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SPCVAL
C SUBROUTINE TO COMPUTE A SPECTRAL VALUE AT FREQUENCY
C FREQ(RADIANS) FOR SEQUENCES X(N) AND N*X(N)
C-----------------------------------------------------------------------
C
      SUBROUTINE SPCVAL(NX, X, FREQ, XR, XI, YR, YI)
      DIMENSION X(1)
      DOUBLE PRECISION U0, U1, U2, W0, W1, W2, A, B, C, D, A1, A2, SA0,
     *    CA0
C
C DESCRIPTION OF ARGUMENTS:
C
C   NX = # VALUES IN SEQUENCE X(N)
C    X = ARRAY CONTAINING INPUT SEQUENCE X(N)
C FREQ = FREQUENCY(RADIANS)
C   XR = REAL PART OF THE SPECTRAL VALUE OF X(N)
C   XI = IMAGINARY PART OF THE SPECTRAL VALUE OF X(N)
C   YR = REAL PART OF THE SPECTRAL VALUE OF NX(N)
C   YI = IMAGINARY PART OF THE SPECTRAL VALUE OF NX(N)
C
C METHOD:
C--MODIFIED GOERTZEL ALGORITHM AS PROPOSED BY BONZANIGO
C
C (IEEE TRAS. ASSP,VOL. 26,NO. 1,FEB 78)
C
C INITIALIZATION
C
      CA0 = DBLE(COS(FREQ))
      SA0 = DBLE(SIN(FREQ))
      A1 = 2.D+0*CA0
      U1 = 0.D+0
      U2 = U1
      W1 = U1
      W2 = U1
C
C MAIN LOOP (GOERTZEL ALGORITHM)
C
      DO 10 J=1,NX
        XJ = DBLE(X(J))
        U0 = XJ + A1*U1 - U2
        W0 = DBLE(FLOAT(J-1))*XJ + A1*W1 - W2
        U2 = U1
        U1 = U0
        W2 = W1
        W1 = W0
  10  CONTINUE
C
C BONZANIGO'S PHASE CORRECTION
C
      A = U1 - U2*CA0
      B = U2*SA0
      C = W1 - W2*CA0
      D = W2*SA0
      A2 = DBLE(FREQ*FLOAT(NX-1))
      U1 = DCOS(A2)
      U2 = -DSIN(A2)
      XR = SNGL(U1*A-U2*B)
      XI = SNGL(U2*A+U1*B)
      YR = SNGL(U1*C-U2*D)
      YI = SNGL(U2*C+U1*D)
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: PHAUNW
C PHASE UNWRAPPING BASED ON TRIBOLET'S ADAPTIVE INTEGRATION SCHEME.
C THE UNWRAPPED PHASE ESTIMATE IS RETURNED IN PHAUNW.
C-----------------------------------------------------------------------
C
      FUNCTION PHAUNW(X, NX, ISNX, I, PPHASE, PPDVT, PPV, PDVT, ISCONS)
C
C DESCRIPTION OF ARGUMENTS:
C
C      X = ARRAY CONTAINING SEQUENCE X(N)
C     NX = # POINTS IN SEQUENCE X(N)
C   ISNX = INTEGER VALUE,EITHER +1 OR -1 DEPENDING
C          ON THE SIGN REVERSAL OF X(N)
C      I = INDEX OF PHASE ESTIMATE ON EQUALLY SPACED FFT
C          FREQUENCY GRID
C PPHASE = PHASE ESTIMATE AT INDEX I-1
C  PPDVT = PHASE DERIVATIVE AT INDEX I-1
C    PPV = PHASE PRINCIPLE VALUE AT INDEX I
C   PDVT = PHASE DERIVATIVE AT INDEX I
C ISCONS = .FALSE. IF PHASE ESTIMATION UNSUCCESSFUL
C          (STACK DIMENSION EXCEEDED BEFORE CONSISTENT ESTIMATE
C          FOUND);.TRUE. OTHERWISE
C
C
C SUBROUTINES CALLED:
C
C SPCVAL = SUBROUTINE TO COMPUTE SPECTRAL VALUE
C PHCHCK = SUBROUTINE TO CHECK PHASE CONSISTENCY
C
C FUNCTIONS CALLED:
C
C PPVPHA = FUNCTION TO COMPUTE PRINCIPLE VALUE OF PHASE
C PHADVT = FUNCTION TO COMPUTE PHASE DERIVATIVE
C AMODSQ = FUNCTION TO COMPUTE MODULUS SQUARED OF A
C          COMPLEX NUMBER
C
C
      DIMENSION SDVT(17), SPPV(17), X(1)
      INTEGER SINDEX(17), PINDEX, SP
      LOGICAL ISCONS, FIRST
      COMMON PI, TWOPI, THLINC, THLCON, NFFT, NPTS, N, L, H, H1, DVTMN2
C
C DESCRIPTION OF ARRAYS AND VARIABLES:
C
C SINDEX = INDEX STACK
C   SDVT = PHASE DERIVATIVE STACK
C   SPPV = PHASE PRINCIPLE VALUE STACK
C     SP = STACK POINTER
C
C
C
C
C
C INITIALIZATION
C
      FIRST = .TRUE.
      PINDEX = 1
      SP = 1
      SPPV(SP) = PPV
      SDVT(SP) = PDVT
      SINDEX(SP) = L + 1
C
C ENTER MAJOR LOOP
C
      GO TO 40
C
C UPDATE PREVIOUS ESTIMATE
C
  10  PINDEX = SINDEX(SP)
      PPHASE = PHASE
      PPDVT = SDVT(SP)
      SP = SP - 1
      GO TO 40
C
C IF SOFTWARE STACK DIMENSION DOES NOT ALLOW FURTHER STEP
C REDUCTION, RETURN.
C
  20  IF ((SINDEX(SP)-PINDEX).GT.1) GO TO 30
      ISCONS = .FALSE.
      PHAUNW = 0.
      RETURN
C
C DEFINE INTERMEDIATE FREQUENCY(I.F.):
C W = (TWOPI/NFFT)*(I-2+(K-1)/L)
C
  30  K = (SINDEX(SP)+PINDEX)/2
C
C CALCULATE I.F.
C
      FREQ = TWOPI*(FLOAT(I-2)*FLOAT(L)+FLOAT(K-1))/H
      CALL SPCVAL(NX, X, FREQ, XR, XI, YR, YI)
C
C COMPUTE PHASE DERIVATIVE AND PRINCIPLE VALUE OF THE PHASE
C AT I.F.;UPDATE STACK
C
      SP = SP + 1
      SINDEX(SP) = K
      SPPV(SP) = PPVPHA(XR,XI,ISNX)
      XMAG = AMODSQ(XR,XI)
      SDVT(SP) = PHADVT(XR,XI,YR,YI,XMAG)
C
C EVALUATE THE PHASE INCREMENT ACROSS SPECTRAL INTERVAL
C
  40  DELTA = H1*FLOAT(SINDEX(SP)-PINDEX)
      PHAINC = DELTA*(PPDVT+SDVT(SP))
C
C IF PHASE INCREMENT,REDUCED BY EXPECTED LINEAR PHASE INCREMENT,
C IS GREATER THAN SPECIFIED THRESHOLD,ADAPT STEP SIZE
C
      IF (ABS(PHAINC-DELTA*DVTMN2).GT.THLINC) GO TO 20
C
C FORM PHASE ESTIMATE;CHECK CONSISTENCY
C
      PHASE = PPHASE + PHAINC
      CALL PHCHCK(PHASE, SPPV(SP), ISCONS)
      IF (.NOT.ISCONS) GO TO 20
C
C IF RESULTING PHASE INCREMENT IS GREATER THAN PI,ADAPT STEP SIZE
C FOR MORE CONFIDENT ESTIMATE;OTHERWISE UPDATE PREVIOUS ESTIMATE
C IF STACK IS NOT EMPTY
C
      IF (ABS(PHASE-PPHASE).GT.PI) GO TO 20
C
C WHEN STACK IS EMPTY,THE UNWRAPPED PHASE AT
C W = TWOPI*(I-1)/NFFT IS HELD IN PHASE
C
      IF (SP.NE.1) GO TO 10
      PHAUNW = PHASE
      RETURN
      END
C
C---------------------------------------------------------------------
C FUNCTION: PPVPHA
C COMPUTE THE PRINCIPLE VALUE OF THE PHASE OF A SPECTRAL VALUE
C---------------------------------------------------------------------
C
      FUNCTION PPVPHA(XR, XI, ISNX)
C
C DESCRIPTION OF ARGUMENTS:
C
C   XR = REAL PART OF THE SPECTRAL VALUE
C   XI = IMAGINARY PART OF THE SPECTRAL VALUE
C ISNX = SIGN OF THE SPECTRAL VALUE AT ZERO FREQUENCY
C
      IF (ISNX.EQ.1) PPVPHA = (ATAN2((XI),(XR)))
      IF (ISNX.EQ.(-1)) PPVPHA = (ATAN2(-(XI),-(XR)))
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: PHADVT
C COMPUTE THE PHASE DERIVATIVE OF A SPECTRAL VALUE OF A SEQUENCE X(N)
C A SPECTRAL VALUE OF A SEQUENCE X(N)
C-----------------------------------------------------------------------
C
      FUNCTION PHADVT(XR, XI, YR, YI, XMAG)
C
C DESCRIPTION OF ARGUMENTS:
C
C   XR = REAL PART OF THE SPECTRAL VALUE OF X(N)
C   XI = IMAGINARY PART OF THE SPECTRAL VALUE OF X(N)
C   YR = REAL PART OF THE SPECTRAL VALUE OF NX(N)
C   YI = IMAGINARY PART OF THE SPECTRAL VALUE OF NX(N)
C XMAG = MAGNITUDE SQUARED OF THE SPECTRAL VALUE OF X(N)
C
      PHADVT = -SNGL((DBLE(XR)*DBLE(YR)+DBLE(XI)*DBLE(YI))/DBLE(XMAG))
      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
C
C-----------------------------------------------------------------------
C SUBROUTINE: PHCHCK
C SUBROUTINE TO CHECK CONSISTENCY OF A PHASE ESTIMATE
C-----------------------------------------------------------------------
C
      SUBROUTINE PHCHCK(PH, PV, ISCONS)
C
C DESCRIPTION OF ARGUMENTS:
C
C     PH = PHASE ESTIMATE
C     PV = PRINCIPLE VALUE OF PHASE AT FREQUENCY
C          OF PHASE ESTIMATE
C ISCONS = .FALSE. IF PHASE ESTIMATE NOT CONSISTENT;
C          .TRUE. IF PHASE ESTIMATE CONSISTENT,
C          PHASE RETURNED IN PH
C
      COMMON PI, TWOPI, THLINC, THLCON, NFFT, NPTS, N, L, H, H1, DVTMN2
      LOGICAL ISCONS
C
C FIND THE TWO ADMISSIBLE PHASE VALUES CLOSEST TO PH
C
      A0 = (PH-PV)/TWOPI
      A1 = FLOAT(IFIX(A0))*TWOPI + PV
      A2 = A1 + SIGN(TWOPI,A0)
      A3 = ABS(A1-PH)
      A4 = ABS(A2-PH)
C
C CHECK CONSISTENCY
C
      ISCONS = .FALSE.
      IF (A3.GT.THLCON .AND. A4.GT.THLCON) RETURN
      ISCONS = .TRUE.
C
C FIND THE CLOSEST UNWRAPPED PHASE ESTIMATE
C
      PH = A1
      IF (A3.GT.A4) PH = A2
      RETURN
      END
Return to Main Software Page