C-------------------------------------------------------------     ************
C                                                                     FILTERG
C                                                                  ************
      SUBROUTINE FILTERG(A,M,D,N,O)
C
C     Computes the convolution of two vectors.
C     Martin J. McBride.  9/28/85
C     General Electric CRD, Information System Operation.
C
      DOUBLE PRECISION A(1),D(1),O(1)
      INTEGER M,N,LOOP,I,J

C  If M and N are improper, then give error message and abort run.
      IF (N .LT. M) THEN
         PRINT*
         PRINT*,'FILTERG called with N less than M.'
         PRINT*
         STOP
      ENDIF
      IF (M .LE. 0) THEN
         PRINT*
         PRINT*,'FILTERG called with M less than or equal to 0.'
         PRINT*
         STOP
      ENDIF

C  Initialization of an upper limit for some of the loops and the output
C  vector to zero.
      LOOP = N - M + 1
      DO 10 I = 1,LOOP
   10    O(I) = 0.0

C  Loops to compute the filtered output.
      DO 20 J = 1,M
         DO 15 I = 1,LOOP
   15       O(I) = O(I) + A(J)*D(I+J-1)
   20 CONTINUE
      RETURN
      END

C-------------------------------------------------------------     ************
C                                                                     FILTERS
C                                                                  ************
      SUBROUTINE FILTERS(A,M,D,N,O)
C
C     Computes the convolution of two vectors and assumes that
C          the filter coefficient vector is symmetric.
C     Martin J. McBride.  9/28/85
C     General Electric CRD, Information System Operation.
C
      DOUBLE PRECISION A(1),D(1),O(1)
      INTEGER M,N,LOOP,I,J,MIDL

C  If M and N are improper, then give error message and abort run.
      IF (N .LT. M) THEN
         PRINT*
         PRINT*,'FILTERS called with N less than M.'
         PRINT*
         STOP
      ENDIF
      IF (M .LE. 0) THEN
         PRINT*
         PRINT*,'FILTERS called with M less than or equal to 0.'
         PRINT*
         STOP
      ENDIF

C  Initialization of an upper limit for some of the loops.
      LOOP = N - M + 1
      IF (MOD(M,2) .NE. 0) GO TO 30

C  Filtering for an even number of filter coefficients.
      DO 10 I = 1,LOOP
         O(I) = 0.0
         DO 20 J = 1,M/2
   20       O(I) = O(I) + A(J) * (D(I+J-1) + D(M+I-J))
   10 CONTINUE
      RETURN

C  Filtering for an odd number of filter coefficients.
   30 MIDL = (M+1)/2
      DO 40 I = 1,LOOP
         O(I) = A(MIDL)*D(I+MIDL-1)
         DO 50 J = 1,MIDL-1
   50       O(I) = O(I) + A(J) * (D(I+J-1) + D(M+I-J))
   40 CONTINUE
      RETURN
      END

C-------------------------------------------------------------     ************
C                                                                     OPFILT
C                                                                  ************
      SUBROUTINE OPFILT(N,A,B,C,R)
C
C  OPFILT IS A RECURSIVE SOLUTION TO THE WIENER-LEVINSON
C  OPTIMAL FILTER EQUATION: R*A=B, WHERE R IS THE AUTO-
C  CORRELATION FUNCTION OF THE INPUT SIGNAL, B IS THE
C  AUTO-CORRELATION FUNCTION FOR THE INFORMATION (USUALLY
C  SHARP: NAMELY B = (1.,0,0,...,0))
C
C        N=NUMBER OF FILTER COEFFICIENTS
C        R=AN N DIMENSIONED ARRAY, SIGNAL AUTO-CORRELATION
C        B=AN N DIMENSIONED ARRAY, INFORMATION AUTO-CORRELATION
C        A=OUTPUT OF FILTER COEFFICIENTS
C        C=SCRATCH ARRAY OF DIMENSION 2*N
C
C  REFERENCE: N. LEVINSON, 'JOURNAL OF MATHEMATICS AND PHYSICS',
C             VOL.XXV, NO.4, JANUARY 1947, PP.261-278
C
C  MODIFIED BY MIKE ESS MARCH 20,1980 BY REMOVING SUBROUTINE CALLS TO
C   CSOLV, ASOLV AND SAXPY AND SUBSTITUTING INLINE CODE.
C
      INTEGER N,M,I,MP1
      DOUBLE PRECISION R(1),B(1),A(1),C(1),C1,SA,SDOT
C
C  FIRST RECURSIVE STEP
C
      A(1)=B(1)*(1.0/R(1))
      C(1)=R(2)*(1.0/R(1))
C
C  NEXT N-1 STEPS ARE ITERATIONS ON THE FILTER LENGTH
C
      DO 80 M = 0,N-2
C
         IF(MOD(M,2) .EQ. 0) GO TO 40
C
         IF(M .EQ. 0) GO TO 20
            C(N+1) = (R(M+2) - SDOT(M,C,1,R(2),1))
     *                 /(R(1) - SDOT(M,C,1,R(2),-1))
            C1 = C(N+1)
CDIR$ IVDEP
            DO 10 I = 2,M+1
               C(N+I) = (C(I-1) * 1.0) - C1 * C(M+2-I)
  10        CONTINUE
  20     CONTINUE
         MP1 = M + 1
         A(M+2) = (B(M+2) - SDOT(MP1,A,1,R(2),-1))
     *              /(R(1) - SDOT(MP1,C(N+1),1,R(2),-1))
         SA = - A(M+2)
         DO 30 I = 1,MP1
            A(I) = (A(I) * 1.0) + SA * C(N+I)
  30     CONTINUE
         GO TO 80
C
  40  CONTINUE
         IF(M .LE. 0) GO TO 60
            C(1) = (R(M+2) - SDOT(M,C(N+1),1,R(2),1))
     *             /(R(1) - SDOT(M,C(N+1),1,R(2),-1))
            C1 = C(1)
CDIR$ IVDEP
            DO 50 I = 2,M+1
               C(I) = (C(N+I-1) * 1.0) - C1 * C(N+M+2-I)
  50        CONTINUE
  60     CONTINUE
         MP1 = M + 1
         A(M+2) = (B(M+2) - SDOT(MP1,A,1,R(2),-1))
     *              /(R(1) - SDOT(MP1,C,1,R(2),-1))
         SA = - A(M+2)
         DO 70 I = 1,MP1
            A(I) = (A(I) * 1.0) + SA * C(I)
  70     CONTINUE
C
  80  CONTINUE
      RETURN
      END