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 REAL 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 REAL 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 REAL R(1),B(1),A(1),C(1),C1,SA 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