C
C-----------------------------------------------------------------------
C  MAIN PROGRAM: TEST PROGRAM FOR COVAR AND AUTO
C  AUTHORS:      A. H. GRAY, JR.  AND J. D. MARKEL
C      GRAY -    UNIV. OF CALIF., SANTA BARBARA, CA 93109
C      BOTH -    SIGNAL TECHNOLOGY, INC., 15 W. DE LA GUERRA STREET
C                SANTA BARBARA, CA 93101
C-----------------------------------------------------------------------
C
      DIMENSION X(100), A(21), RC(21)
      IOUTD = I1MACH(2)
      X(1) = 1.
      X(2) = 1.44
      X(3) = 1.44*X(2) - 1.26
      DO 10 J=4,40
        X(J) = 1.44*X(J-1) - 1.26*X(J-2) + .81*X(J-3)
  10  CONTINUE
9999  FORMAT (1X, I3, 2E15.6)
9998  FORMAT (1X, E12.6/)
      N = 40
      DO 20 M=2,3
        CALL COVAR(N, X, M, A, ALPHA, RC)
        MP = M + 1
        WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
        WRITE (IOUTD,9998) ALPHA
        CALL AUTO(N, X, M, A, ALPHA, RC)
        WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
        WRITE (IOUTD,9998) ALPHA
  20  CONTINUE
      DO 30 JB=1,40
        J = 44 - JB
        X(J) = X(J-3)
  30  CONTINUE
      DO 40 J=1,3
        JB = 43 + J
        X(J) = 0.
        X(JB) = 0.
  40  CONTINUE
      M = 3
      N = 46
      CALL COVAR(N, X, M, A, ALPHA, RC)
      MP = M + 1
      WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
      WRITE (IOUTD,9998) ALPHA
      CALL AUTO(N, X, M, A, ALPHA, RC)
      WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
      WRITE (IOUTD,9998) ALPHA
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  COVAR
C A SUBROUTINE FOR IMPLEMENTING THE COVARIANCE
C METHOD OF LINEAR PREDICTION ANALYSIS
C-----------------------------------------------------------------------
C
      SUBROUTINE COVAR(N, X, M, A, ALPHA, GRC)
C
C   INPUTS:   N - NO. OF DATA POINTS
C             X(N) - INPUT DATA SEQUENCE
C             M - ORDER OF FILTER (M<21, SEE NOTE*)
C  OUTPUTS:   A - FILTER COEFFICIENTS
C             ALPHA - RESIDUAL "ENERGY"
C             A - FILTER COEFFICIENTS
C             GRC - "GENERALIZED REFLECTION COEFFICIENTS",
C
C *PROGRAM LIMITED TO M=20, BECAUSE OF THE DIMENSIONS
C B(M*(M-1)/2), BETA(M), AND CC(M+1)
C
      DIMENSION X(1), A(1), GRC(1)
      DIMENSION B(190), BETA(20), CC(21)
      MP = M + 1
      MT = (MP*M)/2
      DO 10 J=1,MT
        B(J) = 0.
  10  CONTINUE
      ALPHA = 0.
      CC(1) = 0.
      CC(2) = 0.
      DO 20 NP=MP,N
        NP1 = NP - 1
        ALPHA = ALPHA + X(NP)*X(NP)
        CC(1) = CC(1) + X(NP)*X(NP1)
        CC(2) = CC(2) + X(NP1)*X(NP1)
  20  CONTINUE
      B(1) = 1.
      BETA(1) = CC(2)
      GRC(1) = -CC(1)/CC(2)
      A(1) = 1.
      A(2) = GRC(1)
      ALPHA = ALPHA + GRC(1)*CC(1)
      MF = M
      DO 130 MINC=2,MF
        DO 30 J=1,MINC
          JP = MINC + 2 - J
          N1 = MP + 1 - JP
          N2 = N + 1 - MINC
          N3 = N + 2 - JP
          N4 = MP - MINC
          CC(JP) = CC(JP-1) + X(N4)*X(N1) - X(N2)*X(N3)
  30    CONTINUE
        CC(1) = 0.
        DO 40 NP=MP,N
          N1 = NP - MINC
          CC(1) = CC(1) + X(N1)*X(NP)
  40    CONTINUE
        MSUB = (MINC*MINC-MINC)/2
        MM1 = MINC - 1
        N1 = MSUB + MINC
        B(N1) = 1.
        DO 80 IP=1,MM1
          ISUB = (IP*IP-IP)/2
          IF (BETA(IP)) 150, 150, 50
  50      GAM = 0.
          DO 60 J=1,IP
            N1 = ISUB + J
            GAM = GAM + CC(J+1)*B(N1)
  60      CONTINUE
          GAM = GAM/BETA(IP)
          DO 70 JP=1,IP
            N1 = MSUB + JP
            N2 = ISUB + JP
            B(N1) = B(N1) - GAM*B(N2)
  70      CONTINUE
  80    CONTINUE
        BETA(MINC) = 0.
        DO 90 J=1,MINC
          N1 = MSUB + J
          BETA(MINC) = BETA(MINC) + CC(J+1)*B(N1)
  90    CONTINUE
        IF (BETA(MINC)) 150, 150, 100
 100    S = 0.
        DO 110 IP=1,MINC
          S = S + CC(IP)*A(IP)
 110    CONTINUE
        GRC(MINC) = -S/BETA(MINC)
        DO 120 IP=2,MINC
          M2 = MSUB + IP - 1
          A(IP) = A(IP) + GRC(MINC)*B(M2)
 120    CONTINUE
        A(MINC+1) = GRC(MINC)
        S = GRC(MINC)*GRC(MINC)*BETA(MINC)
        ALPHA = ALPHA - S
        IF (ALPHA) 150, 150, 130
 130  CONTINUE
 140  RETURN
 150  CONTINUE
C
C WARNING - SINGULAR MATRIX
C
      IOUTD = I1MACH(2)
      WRITE (IOUTD,9999)
9999  FORMAT (34H WARNING - SINGULAR MATRIX - COVAR)
      GO TO 140
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  AUTO
C A SUBROUTINE FOR IMPLEMENTING THE AUTOCORRELATION
C METHOD OF LINEAR PREDICTION ANALYSIS
C-----------------------------------------------------------------------
C
      SUBROUTINE AUTO(N, X, M, A, ALPHA, RC)
C
C   INPUTS:   N - NO. OF DATA POINTS
C             X(N) - INPUT DATA SEQUENCE
C             M - ORDER OF FILTER (M<21, SEE NOTE*)
C  OUTPUTS:   A - FILTER COEFFICIENTS
C             ALPHA - RESIDUAL "ENERGY"
C             RC - REFLECTION COEFFICIENTS
C
C *PROGRAM LIMITED TO M<21 BECAUSE OF DIMENSIONS OF R(.)
C
      DIMENSION X(1), A(1), RC(1)
      DIMENSION R(21)
      MP = M + 1
      DO 20 K=1,MP
        R(K) = 0.
        NK = N - K + 1
        DO 10 NP=1,NK
          N1 = NP + K - 1
          R(K) = R(K) + X(NP)*X(N1)
  10    CONTINUE
  20  CONTINUE
      RC(1) = -R(2)/R(1)
      A(1) = 1.
      A(2) = RC(1)
      ALPHA = R(1) + R(2)*RC(1)
      DO 50 MINC=2,M
        S = 0.
        DO 30 IP=1,MINC
          N1 = MINC - IP + 2
          S = S + R(N1)*A(IP)
  30    CONTINUE
        RC(MINC) = -S/ALPHA
        MH = MINC/2 + 1
        DO 40 IP=2,MH
          IB = MINC - IP + 2
          AT = A(IP) + RC(MINC)*A(IB)
          A(IB) = A(IB) + RC(MINC)*A(IP)
          A(IP) = AT
  40    CONTINUE
        A(MINC+1) = RC(MINC)
        ALPHA = ALPHA + RC(MINC)*S
        IF (ALPHA) 70, 70, 50
  50  CONTINUE
  60  RETURN
  70  CONTINUE
C
C WARNING - SINGULAR MATRIX
C
      IOUTD = I1MACH(2)
      WRITE (IOUTD,9999)
9999  FORMAT (33H WARNING - SINGULAR MATRIX - AUTO)
      GO TO 60
      END
C