CTITLESAOPTO -- OPTIMUM OPERATOR CALCULATION 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA 00020000 CA AUTHOR D. NYMAN 00030000 CA DESIGNER D. NYMAN 00040000 CA LANGUAGE FORTRAN 00050000 CA SYSTEM IBM AND CRAY 00060000 CA WRITTEN 06/07/78 00070000 CA REVISED 09/12/79 DCN CALCULATE ONLY REQUESTED 00080000 CA OPERATOR WHEN KOPR > 0 00090000 CA REVISED 12/07/83 DCN: STOP FOR NEGATIVE ERROR 00100000 CA REVISED 12/16/86 ESN: FOR THE CRAY. 00110000 CA 00120000 CA 00130000 CA CALL SAOPTO (R, C, OP, OPL, KOPR, H, A, D, E) 00140000 CA 00150000 CA INPUT R = OPL+1 AUTO-CORRELATION VALUES OF INPUT R8 00160000 CA WAVELET. R(K) CORRESPONDS TO A LAG OF K-1. 00170000 CA C = 2*OPL-1 CROSS-CORRELATION VALUES OF INPUT R8 00180000 CA WAVELET AND DESIRED OUTPUT WAVELET. C(K) 00190000 CA CORRESPONDS TO INPUT WAVELET LAGGED BY K-OPL. 00200000 CA OPL = OPERATOR LENGTH IN SAMPLES I4 00210000 CA KOPR = 0: FIND OPTIMUM NUMBER OF ANTICIPATORY TERMS I4 00220000 CA >0: FIND OPERATOR STARTING AT C(KOPR), 00230000 CA WHICH HAS OPL-KOPR ANTICIPATORY TERMS. 00240000 CA H(1) = SUM OF SQUARES OF OUTPUT WAVELET VALUES R8 00250000 CA A = WORK ARRAY OF LENGTH OPL+1 R8 00260000 CA D = WORK ARRAY OF LENGTH OPL+1 R8 00270000 CA E = WORK ARRAY OF LENGTH OPL+1 R8 00280000 CA 00290000 CA OUTPUT OP = OPERATOR ARRAY R8 00300000 CA KOPR = INDEX OF ARRAY C AT WHICH OPERATOR CALCU- I4 00310000 CA LATION WAS BEGUN. NUMBER OF ANTICIPATORY 00320000 CA TERMS IN OPERATOR IS OPL-KOPR. 00330000 CA H = OPERATOR ERROR VS. CALCULATION STARTING R8 00340000 CA INDEX OF ARRAY C. IF A STARTING POSITION 00350000 CA IS SPECIFIED (INPUT KOPT > 0), ONLY THAT 00360000 CA OPERATOR ERROR IS RETURNED. A NEGATIVE ERROR 00370000 CA INDICATES AN UNSTABLE OPERATOR CALCULATION, 00380000 CA AND A VALUE LESS THAN ABOUT 1E-5 MAY BE 00390000 CA CAUSE FOR CONCERN. 00400000 CA 00410000 CA 00420000 CA THIS SUBROUTINE CALCULATES THE OPERATOR OF LENGTH OPL WHICH BEST 00430000 CA TRANSFORMS THE INPUT WAVELET INTO THE DESIRED OUTPUT WAVELET. 00440000 CA 00450000 C=======================================================================00460000 C EJECT 00470000 C=======================================================================00480000 C 00490000 C SAOPTO FINDS THE SOLUTION TO THE MATRIX EQUATION 00500000 C 00510000 C | R(1) R(2) . . R(N-1) R(N) | | OP(1) | | C(K) | 00520000 C | R(2) R(1) . . R(N-2) R(N-1)| | OP(2) | | C(K+1) | 00530000 C | . . . . . . | * | . | = | . | 00540000 C | . . . . . . | | . | | . | 00550000 C |R(N-1) R(N-2) . . R(1) R(2) | |OP(N-1)| |C(K+N-2)| 00560000 C | R(N) R(N-1) . . R(2) R(1) | | OP(N) | |C(K+N-1)| 00570000 C 00580000 C THE VALUE OF K (KOPR) IS FOUND WHICH MINIMIZES THE ERROR IN 00590000 C TRANSFORMING THE INPUT WAVELET INTO THE DESIRED OUTPUT WAVELET. 00600000 C 00610000 C DENOTE THE ABOVE MATRIX EQUATION BY R * OP = CK . IN THIS 00620000 C NOTATION, THE ARRAYS D AND E HAVE THE PROPERTIES: 00630000 C 00640000 C | XXXX | | EPS1 | 00650000 C | 0 | | 0 | 00660000 C R * D = | . | AND R * E = | . | , 00670000 C | 0 | | 0 | 00680000 C | SIGN | | EPSN | 00690000 C 00700000 C WHERE XXXX IS 0 WHEN BUILDING THE OPERATOR FOR K = 1, AND 00710000 C XXXX IS NOT DETERMINED WHEN TRANSLATING THE OPERATOR FROM 00720000 C K = J TO K = J+1 . THE BUILDING RECURSION ALGORITHM YIELDS 00730000 C 00740000 C OP(LENGTH=N) = OP(LENGTH=N-1) + D , 00750000 C 00760000 C WHERE D = U * E + V * E# , AND E# REPRESENTS E IN 00770000 C REVERSE ORDER. THE VALUES U AND V SATISFY 00780000 C U * EPS1 + V * EPSN = 0 AND U * EPSN + V * EPS1 = SIGN. 00790000 C THE TRANSLATION RECURSION ALGORITHM YIELDS 00800000 C 00810000 C OP(POSITION K) = OP(POSITION K-1) + D , 00820000 C 00830000 C WHERE D = U * E + V * E# , AND U AND V SATISFY 00840000 C V * E(N) + OP(1) = 0 AND U * EPSN + V * EPS1 = SIGN. 00850000 C 00860000 C 00870000 C LOCAL VARIABLES AND CONSTANTS OTHER THAN DEFINED ABOVE. 00880000 C 00890000 C BEST = LOGICAL TRUE IF KOPR = 0 L4 00900000 C ERK = OPERATOR ERROR FOR POSITION K R8 00910000 C ER0 = MINIMUM ERROR THROUTH POSITION K R8 00920000 C N = INTERMEDIATE OPERATOR LENGTH I4 00930000 C P0 = DESIRED OUTPUT WAVELET POWER, FOR ERROR NORMALIZATION R8 00940000 C 00950000 C=======================================================================00960000 C EJECT 00970000 C 00980000 SUBROUTINE SAOPTO (R, C, OP, OPL, KOPR, H, A, D, E) 00990000 C 01000000 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 01010000 C 01020001 COMMON /SYSTEM/ SYSTEM 01030001 C 01040001 INTEGER OPL 01050000 INTEGER S1CPCH 01060001 INTEGER SYSTEM 01070001 C 01080001 LOGICAL BEST 01090000 C 01100001 DIMENSION R(1), C(1), OP(1), H(1), A(1), D(1), E(1) 01110000 C 01120000 C SET BEST AND P0 01130000 C 01140000 BEST=KOPR.LE.0 01150000 P0=H(1) 01160000 M0=MAX0(KOPR-1,0) 01170000 J=OPL*2-1-M0 01180000 C 01190000 C FIND FIRST NON-ZERO C(K) 01200000 C 01210000 DO 10 M=1,J 01220000 IF (C(M+M0).NE.0.0) GO TO 20 01230000 10 CONTINUE 01240000 C 01250000 M=J 01260000 C 01270000 20 A(1)=C(M+M0)/R(1) 01280000 EPSN=C(M+M0) 01290000 IF (M.EQ.1) GO TO 50 01300000 C 01310000 C BUILD OPERATOR OF LENGTH M FOR WHICH 01320000 C C(M0+1) THROUGH C(M0+M-1) ARE ZERO 01330000 C 01340000 J=MAX0(M-OPL,0) 01350000 M=M-J 01360000 M0=M0+J 01370000 EPSN2=EPSN*EPSN 01380000 C 01390000 DO 40 N=2,M 01400000 N1=N-1 01410000 C 01420000 C CALCULATE EPS1 AND MOVE OPERATOR ARRAY 01430000 C 01440000 J=N+1 01450000 EPS1=0. 01460000 DO 30 I=1,N1 01470000 J=J-1 01480000 A(J)=A(J-1) 01490000 30 EPS1=A(J)*R(J)+EPS1 01500000 C 01510000 C CALCULATE OPERATOR OF LENGTH N 01520000 C FROM OPERATOR OF LENGTH N-1 01530000 C 01540000 A(1)=0. 01550000 U=EPSN2/(EPSN2-EPS1*EPS1) 01560000 V=-U*EPS1/EPSN 01570000 J=N+1 01580000 K=J/2 01590000 C 01600000 DO 40 I=1,K 01610000 J=J-1 01620000 Z=U*A(I)+V*A(J) 01630000 A(J)=U*A(J)+V*A(I) 01640000 40 A(I)=Z 01650000 C 01660000 C BUILD OPERATOR OF LENGTH OPL 01670000 C FROM OPERATOR OF LENGTH M 01680000 C 01690000 50 DO 60 I=1,M 01700000 60 D(I)=A(I) 01710000 C CALL DUMP ( 60, D, -500, 6 ) 01720005 C 01730000 IF (M.EQ.OPL) GO TO 100 01740000 M=M+1 01750000 E(1)=0. 01760000 C 01770000 DO 90 N=M,OPL 01780000 N1=N-1 01790000 C 01800000 C CALCULATE EPS1 AND SIGN , AND SET ARRAY E 01810000 C 01820000 EPS1=0. 01830000 SIGN=C(N+M0) 01840000 C 01850000 DO 70 I=1,N1 01860000 SIGN=SIGN-R(I+1)*A(N-I) 01870000 EPS1=EPS1+R(I+1)*D(I) 01880000 70 E(I+1)=D(I) 01890000 C 01900000 C CALCULATE D AND OPERATOR ARRAYS OF LENGTH N 01910000 C 01920000 U=SIGN*EPSN/(EPSN*EPSN-EPS1*EPS1) 01930000 V=-U*EPS1/EPSN 01940000 A(N)=0. 01950000 J=N+1 01960000 C 01970000 DO 80 I=1,N 01980000 J=J-1 01990000 D(I)=U*E(I)+V*E(J) 02000000 80 A(I)=A(I)+D(I) 02010000 C 02020000 90 EPSN=SIGN 02030000 C CALL DUMP ( 90, D, -500, 6 ) 02040005 C 02050000 C TRANSLATE OPERATOR TO POSITIONS C(K) , K = 2 02060000 C TO OPL, USING OPERATOR FOR POSITION C(1) 02070000 C 02080000 C CALCULATE EPS1 AND ER0 , AND SET ARRAY E 02090000 C 02100000 100 ER0=P0 02110000 EPS1=0. 02120000 N=OPL+1 02130000 SIGN=C(N+M0) 02140000 J=N+1 02150000 E(1)=0. 02160000 C 02170000 DO 110 I=1,OPL 02180000 J=J-1 02190000 OP(I)=A(I) 02200000 ER0=ER0-A(I)*C(I+M0) 02210000 EPS1=R(I+1)*D(I)+EPS1 02220000 SIGN=SIGN-A(I)*R(J) 02230000 110 E(I+1)=D(I) 02240000 C 02250000 KOPR=M0+1 02260000 C WRITE (6,7123) KOPR 02270007 C7123 FORMAT (5X,'SAOPTO -- KOPR =',5I10) 02280007 H(KOPR)=ER0/P0 02290000 IF (.NOT. BEST) GO TO 1000 02300000 C 02310000 IF (M0.GT.0) THEN 02320001 IF (S1CPCH(SYSTEM,1,'IBM',1,3) .EQ. 0) THEN 02330001 CALL ARSET (H,M0*2,0.0) 02340001 ELSE 02350001 CALL ARSET (H,M0 ,0.0) 02360001 ENDIF 02370001 ENDIF 02380001 IF (KOPR .GE. OPL) GO TO 1000 02390001 C 02400000 EN=D(OPL) 02410000 A(N)=0. 02420000 L=KOPR+1 02430000 C 02440000 DO 140 K=L,OPL 02450000 C 02460000 C TRANSLATE OPERATOR TO POSITION K FROM POSITION K-1.02470000 C CALCULATE CURRENT ERK , AND NEXT SIGN. 02480000 C 02490000 V=-A(1)/EN 02500000 U=(SIGN-V*EPS1)/EPSN 02510000 SIGN=C(K+OPL) 02520000 ERK=P0 02530000 J=N 02540000 C 02550000 DO 120 I=2,N 02560000 J=J-1 02570000 A(I-1)=A(I)+U*E(I)+V*E(J) 02580000 ERK=ERK-A(I-1)*C(K+I-2) 02590000 120 SIGN=SIGN-A(I-1)*R(J+1) 02600000 C 02610000 H(K)=ERK/P0 02620000 C 02630000 IF (ERK.GT.ER0) GO TO 140 02640000 C 02650000 C CURRENT ERROR IS LESS THAN PREVIOUS MINIMUM 02660000 C 02670000 DO 130 I=1,OPL 02680000 130 OP(I)=A(I) 02690000 C 02700000 KOPR=K 02710000 C WRITE (6,7123) KOPR,KOPR 02720007 ER0=ERK 02730000 IF (ERK .LT. 0.0) GO TO 1000 02740000 C 02750000 140 CONTINUE 02760000 C 02770000 1000 RETURN 02780000 C 02790000 C END 02800005 C SUBROUTINE DUMP ( N, X, L, IPR ) 02810005 C REAL X(2) 02820005 C M = L 02830005 C IF (M .GT. 100 ) M = 100 02840005 C IF (L .LT. 0 ) M = -L 02850005 C WRITE (IPR, 90) N,L 02860005 C 90 FORMAT (2X,2I6,'------------------------------------------') 02870005 C WRITE (IPR,100) (X(I),I=1,M) 02880005 C100 FORMAT (2X,10E13.5) 02890005 C RETURN 02900005 END 02910004