C 00010001 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CTITLESAALPHT -- SMOOTH THE ARRAY WITH A ALPHA TRIM FILTER 00020001 CA AUTHOR S. BICKEL 00030001 CA DESIGNER B. S BOK 00040001 CA LANGUAGE FORTRAN 00050001 CA SYSTEM IBM / CRAY 00060001 CA WRITTEN AUGUST, 1990 00070001 C REVISED 12-20-91 JJC - MODIFIED TO MEET CDP STANDARDS. C 00080001 CA 00090001 CA CALL SAALPHT(NY, Y, NF, X, IM, AL) 00100001 CA 00110001 CA IN/OUT Y = INPUT ARRAY TO BE SMOOTHED R4 00120001 CA INPUT NY = NUMBER OF ELEMENTS IN Y I4 00130001 CA INPUT NF = THE LENGTH OF THE FILTER I4 00140001 CA INPUT AL = FILTER PARAMETER ALPHA (0 <= AL <= 0.5) R4 00150001 CA (IF AL = .0 MOVING AVERAGE FILTER 00160001 CA AL = .5 MOVING MEDIAN FILTER ) 00170001 CA INPUT X = WORK ARRAY WITH NF ELEMENTS R4 00180001 CA INPUT IM = INDEX WORK ARRAY WITH NF ELEMENTS I4 00190001 CA 00200001 CA THIS ROUTINE IS TO SMOOTH THE ARRAY Y WITH A ALPHA TRIM FILTER 00210001 CA 00220001 C************************************************************** C * C SUBROUTINES AND FUNCTIONS CALLED FROM THIS ROUTINE * C * C SAHPSRT * C * C************************************************************** SUBROUTINE SAALPHT(NY,Y,NF,X,IM,AL) 00060000 C IMPLICIT INTEGER(A-Z) 00200001 C INTEGER NY,NF,IM(NF),LOCATE 00210000 INTEGER N,N2,N1,N3,J,L1,IP,K,IS,L,JL,JU,JM 00220000 INTEGER KK,T,TRIM 00230000 REAL Y(NY),X(NF),ALPHA,AL 00240000 C 00250001 N = MIN ( NF, NY ) ALPHA = MIN ( AL, .5 ) ALPHA = MAX ( ALPHA, .0 ) TRIM = ALPHA * N T = N - 2 * TRIM TRIM = T / 2 N2 = N / 2 IF (N2 .EQ. 0) RETURN DO 100 J = 1, N X(J) = Y(J) IM(J) = J 100 CONTINUE CALL SAHPSRT ( N, X, IM ) DO 120 J = 1, N Y(J) = IM(J) IM(J) = J 120 CONTINUE CALL SAHPSRT ( N, Y, IM ) N1 = N2 + 1 N3 = N2 + N1 IF (N3 .EQ. N) THEN Y(N - N2) = X(N1) IF (TRIM .GT. 0) THEN DO 140 J = 1, TRIM Y(N - N2) = Y(N - N2) + X(N1 + J) + X(N1 - J) 140 CONTINUE ENDIF ELSE Y(N - N2) = 0.0 IF (T .EQ. 0) THEN TRIM = 1 T = 2 ENDIF DO 160 J = 1, TRIM Y(N - N2) = Y(N - N2) + X(N2 - J + 1) + X(N1 + J - 1) 160 CONTINUE ENDIF Y(N - N2) = Y(N - N2) / T DO 180 J = 1, N - N1 Y(J) = Y(N - N2) 180 CONTINUE C 00620000 C INITIALIZATION COMPLETE; BEGIN MAIN LOOP. FIND THE INDEX FOY Y(J) 00630000 C 00640000 DO 320 J = N + 1, NY JL = 0 JU = N + 1 200 CONTINUE IF (JU - JL .GT. 1) THEN JM = ISHFT ( JU + JL, -1 ) IF (Y(J) .GT. X(JM)) THEN JL = JM ELSE JU = JM ENDIF GO TO 200 ENDIF L1 = JL + 1 C 00780000 C Y(J) HAS BEEN LOCATED. NOW UPDATA THE INDEX ARRAY 00790000 C 00800000 IP = IM(1) IS = IP IF (IP .GE. L1) IP = IP + 1 DO 220 K = 1, N - 1 IF (IM(K + 1) .GE. L1) THEN IM(K) = IM(K + 1) + 1 ELSE IM(K) = IM(K + 1) ENDIF IF (IM(K) .GT. IP) IM(K) = IM(K) - 1 220 CONTINUE IM(N) = L1 IF (L1 .GT. IP) IM(N) = L1 - 1 C 00930000 C IM(J) HAS BEEN UPDATED. NEXT UPDATA THE X ARRAY 00940000 C 00950000 L = L1 IF (IS .LT. L) L = L - 1 IF (IS .LT. L) THEN DO 240 K = IS, L - 1 X(K) = X(K + 1) 240 CONTINUE ELSE IF (IS .GT. L) THEN DO 260 K = IS, L + 1, -1 X(K) = X(K - 1) 260 CONTINUE ENDIF X(L) = Y(J) C 01060000 C X(J) HAS BEEN UPDATED. 01070000 C NOW SET THE OUTPUT EQUAL TO THE MEDIAN OF THE X ARRAY 01080000 C 01090000 IF (N3 .EQ. N) THEN Y(J - N2) = X(N1) C 01120000 C NOW FILL IN THE ALPHA TRIM FOR N ODD. 01130000 C 01140000 IF (TRIM .GT. 0) THEN DO 280 K = 1, TRIM Y(J - N2) = Y(J - N2) + X(N1 + K) + X(N1 - K) 280 CONTINUE ENDIF ELSE Y(J - N2) = 0.0 C 01210000 C NOW FILL IN THE ALPHA TRIM FOR N EVEN. 01220000 C 01230000 DO 300 K = 1, TRIM Y(J - N2) = Y(J - N2) + X(N1 - K) + X(N2 + K) 300 CONTINUE ENDIF 320 CONTINUE C 01280000 C MAIN LOOP COMPLETE; SCALE THE ARRAY - THIS LOOP SHOULD VECTORIZE 01290000 C 01300000 IF (T .GT. 1) THEN DO 340 J = N + 1, NY Y(J - N2) = Y(J - N2) / T 340 CONTINUE ENDIF C 01350000 C FILL IN THE END OF THE ARRAY 01360000 C 01370000 DO 360 J = NY - N2 + 1, NY Y(J) = Y(NY - N2) 360 CONTINUE IF (N3 .EQ. N) RETURN C 01410000 C INTERPOLATE POINTS FOR EVEN VALUES OF N 01420000 C 01430000 DO 380 J = NY, 2, -1 Y(J) = (Y(J) + Y(J - 1) ) / 2. 380 CONTINUE C C RETURN END