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