Make your own free website on Tripod.com
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