Efficient Lattice Methods For Linear Prediction
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: EFFICIENT LATTICE PROGRAMS FOR LINEAR PREDICTION
C AUTHORS: R.VISWANATHAN AND J.MAKHOUL
C BOLT BERANEK AND NEWMAN INC
C 50 MOULTON STREET
C CAMBRIDGE, MA 02138
C
C INPUT: ALL INPUT DATA ARE PROVIDED IN DATA STATEMENTS
C-----------------------------------------------------------------------
C
DIMENSION PHI(20,20), A(20), SCR(20), SIG(100)
REAL K(20)
DATA NSTAGE, NMAX, NSIG, NFLAG /4,20,32,0/
C
C SIGNAL DATA FOR THE TEST EXAMPLES
C
DATA SIG(1), SIG(2), SIG(3), SIG(4) /-16.0,-18.0,-21.0,-24.0/
DATA SIG(5), SIG(6), SIG(7), SIG(8) /-27.0,-30.0,-31.0,-33.0/
DATA SIG(9), SIG(10), SIG(11), SIG(12) /-36.0,-39.0,-42.0,-44.0/
DATA SIG(13), SIG(14), SIG(15), SIG(16) /-46.0,-49.0,-51.0,-55.0/
DATA SIG(17), SIG(18), SIG(19), SIG(20) /-58.0,-61.0,-64.0,-67.0/
DATA SIG(21), SIG(22), SIG(23), SIG(24) /-71.0,-73.0,-68.0,-53.0/
DATA SIG(25), SIG(26), SIG(27), SIG(28) /-38.0,-36.0,-31.0,-13.0/
DATA SIG(29), SIG(30), SIG(31), SIG(32) /6.0,15.0,19.0,35.0/
C
C SET UP MACHINE CONSTANT (OUTPUT UNIT NUMBER)
C
NO = I1MACH(2)
C
C FIRST TEST EXAMPLE: NON-TOEPLITZ COVARIANCE MATRIX CASE
C
WRITE (NO,9999)
9999 FORMAT (1H1, 4X, 45H***FIRST TEST EXAMPLE. NON-TOEPLITZ COVARIANC,
* 16HE MATRIX CASE***)
C
C COVARIANCE MATRIX COMPUTATION USING METHOD 1
C
10 WRITE (NO,9998)
9998 FORMAT (//8X, 48HSOLUTION USING COVARIANCE MATRIX COMPUTED VIA ME,
* 6HTHOD 1)
WRITE (NO,9997)
9997 FORMAT (//44H REFLECTION COEFFICIENTS FROM PROGRAM CLHARM)
CALL COVAR1(SIG, NSIG, NSTAGE, PHI, NMAX)
DO 20 M=1,NSTAGE
CALL CLHARM(PHI, NMAX, M, A, K(M), ERROR, SCR)
20 CONTINUE
ERROR = ERROR/(2.*PHI(1,1))
WRITE (NO,9996) (K(I),I=1,NSTAGE)
9996 FORMAT (/1X, 4(E15.8, 2X))
WRITE (NO,9995) ERROR
9995 FORMAT (/30H FORWARD-PLUS-BACKWARD ERROR =, E15.8)
WRITE (NO,9994)
9994 FORMAT (//44H REFLECTION COEFFICIENTS FROM PROGRAM COVLAT,
* 22H (R=0, GEOMETRIC MEAN))
R = 0
IFLAG = 0
DO 30 M=1,NSTAGE
CALL COVLAT(PHI, NMAX, M, R, IFLAG, A, K(M), FERROR, BERROR,
* SCR)
30 CONTINUE
FERROR = FERROR/PHI(1,1)
BERROR = BERROR/PHI(1,1)
WRITE (NO,9996) (K(I),I=1,NSTAGE)
WRITE (NO,9993) FERROR, BERROR
9993 FORMAT (/16H FORWARD ERROR =, E15.8, 18H, BACKWARD ERROR =, E15.8)
WRITE (NO,9992)
9992 FORMAT (//39H REFLECTION COEFFS. FROM PROGRAM COVLAT, 9H (R=-INFI,
* 21HNITY, MINIMUM METHOD))
IFLAG = -1
DO 40 M=1,NSTAGE
CALL COVLAT(PHI, NMAX, M, R, IFLAG, A, K(M), FERROR, BERROR,
* SCR)
40 CONTINUE
FERROR = FERROR/PHI(1,1)
BERROR = BERROR/PHI(1,1)
WRITE (NO,9996) (K(I),I=1,NSTAGE)
WRITE (NO,9993) FERROR, BERROR
C
C COVARIANCE MATRIX COMPUTATION USING METHOD 2
C
WRITE (NO,9991)
9991 FORMAT (//8X, 48HSOLUTION USING COVARIANCE MATRIX COMPUTED VIA ME,
* 6HTHOD 2)
WRITE (NO,9997)
DO 50 M=1,NSTAGE
CALL COVAR2(SIG, NSIG, M, PHI, NMAX)
CALL CLHARM(PHI, NMAX, M, A, K(M), ERROR, SCR)
50 CONTINUE
ERROR = ERROR/(2.*PHI(1,1))
WRITE (NO,9996) (K(I),I=1,NSTAGE)
WRITE (NO,9995) ERROR
WRITE (NO,9994)
R = 0
IFLAG = 0
DO 60 M=1,NSTAGE
CALL COVAR2(SIG, NSIG, M, PHI, NMAX)
CALL COVLAT(PHI, NMAX, M, R, IFLAG, A, K(M), FERROR, BERROR,
* SCR)
60 CONTINUE
FERROR = FERROR/PHI(1,1)
BERROR = BERROR/PHI(1,1)
WRITE (NO,9996) (K(I),I=1,NSTAGE)
WRITE (NO,9993) FERROR, BERROR
WRITE (NO,9992)
IFLAG = -1
DO 70 M=1,NSTAGE
CALL COVAR2(SIG, NSIG, M, PHI, NMAX)
CALL COVLAT(PHI, NMAX, M, R, IFLAG, A, K(M), FERROR, BERROR,
* SCR)
70 CONTINUE
FERROR = FERROR/PHI(1,1)
BERROR = BERROR/PHI(1,1)
WRITE (NO,9996) (K(I),I=1,NSTAGE)
WRITE (NO,9993) FERROR, BERROR
IF (NFLAG.EQ.1) GO TO 110
C
C SECOND TEST EXAMPLE: STATIONARY OR TOEPLITZ MATRIX CASE
C
WRITE (NO,9990)
9990 FORMAT (1H1, 4X, 45H***SECOND TEST EXAMPLE. STATIONARY OR TOEPLIT,
* 16HZ MATRIX CASE***)
C
C AUGMENT THE SIGNAL VECTOR OF THE FIRST EXAMPLE WITH NSTAGE ZEROS
C SO THAT THE RESULTING COVARIANCE MATRIX IS TOEPLITZ
C
N1 = NSIG + NSTAGE
NSTP1 = NSTAGE + 1
N1P1 = N1 + 1
N = N1
80 NMNST = N - NSTAGE
SIG(N) = SIG(NMNST)
N = N - 1
IF (N.GE.NSTP1) GO TO 80
DO 90 N=1,NSTAGE
SIG(N) = 0
90 CONTINUE
N2 = N1 + NSTAGE
DO 100 N=N1P1,N2
SIG(N) = 0
100 CONTINUE
NSIG = N2
NFLAG = 1
GO TO 10
C
110 STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COVLAT
C GENERAL COVARIANCE LATTICE ROUTINE; IT COMPUTES THE REFLECTION
C COEFFICIENT OF STAGE M,GIVEN THE COVARIANCE MATRIX OF THE SIGNAL AND
C PREDICTOR COEFFICIENTS UP TO STAGE M-1.
C-----------------------------------------------------------------------
C
SUBROUTINE COVLAT(PHI, NMAX, M, R, IFLAG, A, K, FERROR, BERROR,
* SCR)
DIMENSION PHI(NMAX,NMAX), A(NMAX), SCR(NMAX)
REAL K, KF, KB
C
IF (M.GT.1) GO TO 20
C
C EXPLICIT COMPUTATION OF THE FIRST STAGE REFLECTION COEFFICIENT
C
F = PHI(1,1)
B = PHI(2,2)
C = PHI(1,2)
K = 0
IF (C.EQ.0.0) GO TO 10
KF = -C/B
KB = -C/F
K = GMEAN(KF,KB,R,IFLAG)
10 A(1) = K
GO TO 90
C
C RECURSIVE COMPUTATION OF THE M-TH STAGE (M.GE.2) REFLECTION
C COEFFICIENT
C
20 MP1 = M + 1
MM1 = M - 1
SUM1 = 0
SUM2 = 0
SUM3 = 0
SUM4 = 0
SUM5 = 0
SUM6 = 0
DO 30 I=1,MM1
IP1 = I + 1
SCR(I) = A(I)
SUM1 = SUM1 + A(I)*PHI(1,IP1)
MP1MI = MP1 - I
SUM2 = SUM2 + A(I)*PHI(MP1,MP1MI)
SUM3 = SUM3 + A(I)*(PHI(1,MP1MI)+PHI(IP1,MP1))
Y = A(I)**2
SUM4 = SUM4 + Y*PHI(IP1,IP1)
SUM5 = SUM5 + Y*PHI(MP1MI,MP1MI)
SUM6 = SUM6 + Y*PHI(IP1,MP1MI)
30 CONTINUE
SUM7 = 0
SUM8 = 0
SUM9 = 0
IF (M.EQ.2) GO TO 60
MM2 = M - 2
DO 50 I=1,MM2
IP1 = I + 1
MP1MI = MP1 - I
DO 40 J=IP1,MM1
Y = A(I)*A(J)
SUM7 = SUM7 + Y*PHI(IP1,J+1)
MP1MJ = MP1 - J
SUM8 = SUM8 + Y*PHI(MP1MI,MP1MJ)
SUM9 = SUM9 + Y*(PHI(IP1,MP1MJ)+PHI(J+1,MP1MI))
40 CONTINUE
50 CONTINUE
60 F = PHI(1,1) + 2.*(SUM1+SUM7) + SUM4
B = PHI(MP1,MP1) + 2.*(SUM2+SUM8) + SUM5
C = PHI(1,MP1) + SUM3 + SUM6 + SUM9
K = 0
IF (C.EQ.0.0) GO TO 70
KF = -C/B
KB = -C/F
K = GMEAN(KF,KB,R,IFLAG)
70 CONTINUE
C
C INSERT CALL TO USER'S OWN SUBROUTINE TO QUANTIZE THE REFLECTION
C COEFFICIENT K
C
C RECURSION TO CONVERT REFLECTION COEFFICIENTS TO PREDICTOR
C COEFFICIENTS
C
DO 80 I=1,MM1
MMI = M - I
A(I) = SCR(I) + K*SCR(MMI)
80 CONTINUE
A(M) = K
C
C COMPUTE FORWARD AND BACKWARD ERRORS AT THE OUTPUT OF THE M-TH
C LATTICE STAGE
C
90 X = K**2
Y = 2.*K*C
FERROR = F + X*B + Y
BERROR = B + X*F + Y
C
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: GMEAN
C FUNCTION TO COMPUTE THE R-TH GENERALIZED MEAN BETWEEN KF AND KB
C-----------------------------------------------------------------------
C
FUNCTION GMEAN(KF, KB, R, IFLAG)
REAL KF, KB
C
AKF = ABS(KF)
AKB = ABS(KB)
S = KF/AKF
IF (IFLAG) 10, 20, 50
C
C MINIMUM METHOD
C
10 GMEAN = AKF
IF (AKB.LT.AKF) GMEAN = AKB
GMEAN = GMEAN*S
RETURN
C
C GENERAL FINITE R-TH MEAN
C
20 IF (R.EQ.(-1.0)) GO TO 30
IF (R.EQ.0.0) GO TO 40
X = (AKF)**R + (AKB)**R
GMEAN = S*((X/2.)**(1./R))
RETURN
C
C HARMONIC MEAN METHOD
C
30 GMEAN = 2.*KF*KB/(KF+KB)
RETURN
C
C GEOMETRIC MEAN METHOD
C
40 GMEAN = S*SQRT(AKF*AKB)
RETURN
C
C MAXIMUM METHOD
C
50 GMEAN = AKF
IF (AKB.GT.AKF) GMEAN = AKB
GMEAN = GMEAN*S
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CLHARM
C COVARIANCE LATTICE ROUTINE FOR HARMONIC MEAN METHOD; IT COMPUTES THE
C REFLECTION COEFFICIENT OF STAGE M, GIVEN THE COVARIANCE MATRIX OF THE
C SIGNAL AND PREDICTOR COEFFICIENTS UP TO STAGE M-1.
C-----------------------------------------------------------------------
C
SUBROUTINE CLHARM(PHI, NMAX, M, A, K, ERROR, SCR)
DIMENSION PHI(NMAX,NMAX), A(NMAX), SCR(NMAX)
REAL K
C
IF (M.GT.1) GO TO 20
C
C EXPLICIT COMPUTATION OF THE FIRST STAGE REFLECTION COEFFICIENT
C
FPLUSB = PHI(1,1) + PHI(2,2)
C = PHI(1,2)
K = 0
IF (C.EQ.0.0) GO TO 10
K = -2.*C/FPLUSB
10 A(1) = K
GO TO 90
C
C RECURSIVE COMPUTATION OF THE M-TH STAGE (M.GE.2) REFLECTION
C COEFFICIENT
C
20 MP1 = M + 1
MM1 = M - 1
SUM1 = 0
SUM3 = 0
SUM4 = 0
SUM6 = 0
DO 30 I=1,MM1
IP1 = I + 1
SCR(I) = A(I)
MP1MI = MP1 - I
SUM1 = SUM1 + A(I)*(PHI(1,IP1)+PHI(MP1,MP1MI))
SUM3 = SUM3 + A(I)*(PHI(1,MP1MI)+PHI(IP1,MP1))
Y = A(I)**2
SUM4 = SUM4 + Y*(PHI(IP1,IP1)+PHI(MP1MI,MP1MI))
SUM6 = SUM6 + Y*PHI(IP1,MP1MI)
30 CONTINUE
SUM7 = 0
SUM9 = 0
IF (M.EQ.2) GO TO 60
MM2 = M - 2
DO 50 I=1,MM2
IP1 = I + 1
MP1MI = MP1 - I
DO 40 J=IP1,MM1
Y = A(I)*A(J)
MP1MJ = MP1 - J
SUM7 = SUM7 + Y*(PHI(IP1,J+1)+PHI(MP1MI,MP1MJ))
SUM9 = SUM9 + Y*(PHI(IP1,MP1MJ)+PHI(J+1,MP1MI))
40 CONTINUE
50 CONTINUE
60 FPLUSB = PHI(1,1) + PHI(MP1,MP1) + 2.*(SUM1+SUM7) + SUM4
C = PHI(1,MP1) + SUM3 + SUM6 + SUM9
K = 0
IF (C.EQ.0.0) GO TO 70
K = -2.*C/FPLUSB
70 CONTINUE
C
C INSERT CALL TO USER'S OWN SUBROUTINE TO QUANTIZE THE REFLECTION
C COEFFICIENT K
C
C RECURSION TO CONVERT REFLECTION COEFFICIENTS TO PREDICTOR
C COEFFICIENTS
C
DO 80 I=1,MM1
MMI = M - I
A(I) = SCR(I) + K*SCR(MMI)
80 CONTINUE
A(M) = K
C
C COMPUTE MINIMUM FORWARD-PLUS-BACKWARD ERROR AT THE OUTPUT OF
C THE M-TH LATTICE STAGE
C
90 ERROR = FPLUSB*(1.-K**2)
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COVAR1
C GIVEN THE SIGNAL, THIS ROUTINE GENERATES ITS COVARIANCE MATRIX
C USING THE RELATION: PHI(I,J)=SUM?SIG(K+1-I)*SIG(K+1-J)! OVER THE
C RANGE FROM K=NSTAGE+1 TO K=NSIG.
C-----------------------------------------------------------------------
C
SUBROUTINE COVAR1(SIG, NSIG, NSTAGE, PHI, NMAX)
DIMENSION SIG(NSIG), PHI(NMAX,NMAX)
C
NP1 = NSTAGE + 1
NP2 = NSTAGE + 2
NSP2 = NSIG + 2
DO 20 J=1,NP1
TEMP = 0
DO 10 K=NP1,NSIG
KP1MJ = K + 1 - J
TEMP = TEMP + SIG(K)*SIG(KP1MJ)
10 CONTINUE
PHI(1,J) = TEMP
20 CONTINUE
DO 50 I=2,NP1
IM1 = I - 1
DO 30 J=1,IM1
PHI(I,J) = PHI(J,I)
30 CONTINUE
NP2MI = NP2 - I
NSP2MI = NSP2 - I
DO 40 J=I,NP1
NP2MJ = NP2 - J
NSP2MJ = NSP2 - J
PHI(I,J) = PHI(I-1,J-1) + SIG(NP2MI)*SIG(NP2MJ) -
* SIG(NSP2MI)*SIG(NSP2MJ)
40 CONTINUE
50 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COVAR2
C GIVEN THE SIGNAL, THIS ROUTINE COMPUTES THE (M+1) X (M+1) COVARIANCE
C MATRIX CORRESPONDING TO THE LATTICE STAGE M. FOR M LARGER THAN
C 1, THIS COMPUTATION IS DONE EFFICIENTLY, USING THE COVARIANCE
C MATRIX FOR STAGE (M-1).
C-----------------------------------------------------------------------
C
SUBROUTINE COVAR2(SIG, NSIG, M, PHI, NMAX)
DIMENSION SIG(NSIG), PHI(NMAX,NMAX)
C
IF (M.GT.1) GO TO 20
C
C FIRST STAGE: M=1
C
TEMP1 = 0
TEMP2 = 0
DO 10 K=2,NSIG
TEMP1 = TEMP1 + SIG(K)**2
TEMP2 = TEMP2 + SIG(K)*SIG(K-1)
10 CONTINUE
PHI(1,1) = TEMP1
PHI(1,2) = TEMP2
PHI(2,1) = TEMP2
PHI(2,2) = TEMP1 + SIG(1)**2 - SIG(NSIG)**2
RETURN
C
C M-TH STAGE, M.GE.2
C
20 MP1 = M + 1
NSP1 = NSIG + 1
NSM = NSP1 - M
DO 30 J=2,MP1
NSJ = NSP1 + 1 - J
PHI(MP1,J) = PHI(M,J-1) - SIG(NSM)*SIG(NSJ)
PHI(J,MP1) = PHI(MP1,J)
30 CONTINUE
TEMP1 = 0
DO 40 K=MP1,NSIG
KMM = K - M
TEMP1 = TEMP1 + SIG(K)*SIG(KMM)
40 CONTINUE
PHI(MP1,1) = TEMP1
PHI(1,MP1) = TEMP1
DO 60 I=1,M
MP1MI = MP1 - I
DO 50 J=1,M
MP1MJ = MP1 - J
PHI(I,J) = PHI(I,J) - SIG(MP1MI)*SIG(MP1MJ)
PHI(J,I) = PHI(I,J)
50 CONTINUE
60 CONTINUE
C
RETURN
END
Return to Main Software Page