CTITLESAGEND -- DETERMINES DECONVOLUTION FILTER FOR GEND 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D. R. MARTINEZ 00000200 CA DESIGNER R. D. KNIGHT 00000300 CA LANGUAGE FORTRAN H 00000400 CA SYSTEM S/370 00000500 CA WRITTEN MAR 1982 00000600 C REVISED JUN 82 BY RDK FOR IMPLEMENTATION IN SPARC. 00000700 C REVISED 7/21/82 BY RDK ; CHANGED CALL TO M1EURK. PADDED WAVE-00000800 C LET TO FILTER LENGTH IF NECESSARY. 00000900 C REVISED 8/05/82 BY RDK ; ADDED REC, TRC TO CALL LISTS FOR 00001000 C SAGEND, M1EURK, AND SAGFLT. 00001100 C REVISED 3/15/83 BY RDK ; ADDED SCALING FOR MULTIPLE WAVELETS 00001200 C REVISED 5/24/83 BY RDK ; INHIBIT FILTERING FOR ALL PASS 00001300 C FREQUENCY LIMITS. 00001400 CA 00001500 CA 00001600 CA CALL SAGEND ( AUTOCR, DATFLT, DCNFLT, DT, FH, FL, INVFLT,IPR, 00001700 CA IPRNT, LAUTO, LFLT, LFLTR, LFL2, LSH, LWAV, XN, 00001800 CA OUTFLT, SCR, WAVLET, WORK, S1, S2, ICDP, ITRC, 00001900 CA NOFILT ) 00002000 CA 00002100 CA IN/OUT ARGUMENT TYPE DESCRIPTION 00002200 CA 00002300 CA IN AUTOCR R4 ONE-SIDED AUTOCORRELATION FUNCTION 00002400 CA OUT DATFLT R4 ARRAY FOR LEAST SQUARES INVERSE FILTER 00002500 CA OUT DCNFLT R4 ARRAY FOR ERROR FUNCTION 00002600 CA IN DT R4 SAMPLE RATE OF DATA IN MILLISECONDS 00002700 CA IN FH R4 HIGH FREQUENCY TO PASS (HERTZ) 00002800 CA IN FL R4 LOW FREQUENCY TO PASS (HERTZ) 00002900 CA OUT INVFLT R4 BUFFER FOR CORRECTION FILTER 00003000 CA IN IPR I4 PRINT UNIT NUMBER 00003100 CA IN/OUT IPRNT I4 PRINT SWITCH/ERROR FLAG 00003200 CA IN LAUTO I4 LENGTH OF 1-SIDED AUTOCORRELATION FUNCTION00003300 CA IN/OUT LFLT I4 LENGTH OF OUTFLT OPERATOR IN SAMPLES 00003400 CA IN LFLTR I4 LENGTH OF MINIMUM PHASE FILTER 00003500 CA IN/OUT LFL2 I4 LENGTH OF ERROR FUNCTION 00003600 CA OUT LSH I4 PHASE SHIFT IN SAMPLES INTRODUCED BY DCON 00003700 CA IN LWAV I4 INPUT REFERENCE WAVELET LENGTH IN SAMPLES 00003800 CA IN XN R4 WHITE NOISE FACTOR 00003900 CA OUT OUTFLT R4 DECON FILTER TO BE APPLIED TO INPUT TRACE 00004000 CA IN SCR R4 SCRATCH ARRAY USE BY SUBROUTINE 00004100 CA IN WAVLET R4 INPUT REFERENCE WAVELET 00004200 CA IN WORK R4 SCRATCH ARRAY USED BY PROCESS 00004300 CA IN S1 R4 BUFFER FOR 3838 OPERATIONS LENGTH 8200 00004400 CA MUST BE SAME ARRAY ON EACH CALL TO SAGEND 00004500 CA IN S2 R4 BUFFER FOR 3838 OPERATIONS LENGTH 8200 00004600 CA MUST BE SAME ARRAY ON EACH CALL TO SAGEND 00004700 CA IN ICDP I4 CURRENT CDP OR SHOT NUMBER 00004800 CA IN ITRC I4 CURRENT TRACE NUMBER 00004900 CA IN NOFILT I4 BANDPASS FILTERING INHIBITOR 00005000 CA = 0 --> FILTER 00005100 CA = 1 --> DON'T FILTER 00005200 CA 00005300 CA THIS SUBROUTINE DETERMINES THE CORRECT DECONVOLUTION OPERATOR 00005400 CA TO BE APPLIED TO THE CURRENT INPUT TRACE BY PROCESS GEND. 00005500 CA 00005600 CA 00005700 CA SUBROUTINES CALLED: ARDVFC 00005800 CA ARFLP 00005900 CA ARMVE 00006000 CA ARSET 00006100 CA SAGFLT 00006200 CA SAGACC 00006300 CA SAGAC3 00006400 CA M1EURK 00006500 C 00006600 C 00006700 SUBROUTINE SAGEND ( AUTOCR, DATFLT, DCNFLT, DT, FH, FL, 00006800 * INVFLT, IPR, IPRNT, LAUTO, LFLT, LFLTR, 00006900 * LFL2, LSH, LWAV, XN, OUTFLT, SCR, WAVLET, 00007000 * WORK, S1, S2, ICDP, ITRC, NOFILT ) 00007100 C 00007200 DIMENSION AUTOCR ( 1) 00007300 DIMENSION DATFLT ( 1) 00007400 DIMENSION DCNFLT ( 1) 00007500 DIMENSION INVFLT ( 1) 00007600 DIMENSION OUTFLT ( 1) 00007700 DIMENSION SCR ( 1) 00007800 DIMENSION S1 ( 1) 00007900 DIMENSION S2 ( 1) 00008000 DIMENSION WAVLET ( 1) 00008100 DIMENSION WORK ( 1) 00008200 C 00008300 EXTERNAL SAGACC 00008400 C 00008500 REAL INVFLT 00008600 REAL*8 SUMSQ1 00008700 REAL*8 SUMSQ2 00008800 C 00008900 DATA ICDPN /999999/ 00009000 DATA ITRCN /999999/ 00009100 C 00009200 IER = 0 00009300 LAUTO2 = LAUTO*2 00009400 C 00009500 C ==================================================================== 00009600 C DEVELOP THE MINIMUM-PHASE FILTER 00009700 C ==================================================================== 00009800 C 00009900 CALL ARMVE(AUTOCR,SCR(LAUTO),LAUTO ) 00010000 CALL ARFLP(AUTOCR,SCR( 1),LAUTO ) 00010100 SCR(LAUTO2) = 0.0 00010200 C 00010300 IF ( NOFILT.EQ.0 ) 00010400 *CALL SAGFLT(LAUTO2,SCR,FL,FH,DT,S1,ICDP,ITRC) 00010500 C 00010600 SCR(LAUTO) = (1.+XN)*SCR(LAUTO) 00010700 WORK(0001) = 1.0 00010800 CALL ARSET (WORK(2),LFLTR-1,0.0) 00010900 CALL M1EURK ( SCR(LAUTO),WORK,DATFLT,LFLTR,S1,S2,1,ICDP,ITRC) 00011000 C 00011100 C LWAV > 0 ===> CALCULATE PHASE CORRECTION OPERATOR 00011200 C LWAV = 0 ===> STANDARD DECON 00011300 C LWAV < 0 ===> PHASE CORRECTION OPERATOR AVAILABLE 00011400 C 00011500 IF ( LWAV.LE.0 ) GO TO 70 00011600 C 00011700 C ==================================================================== 00011800 C DESIGN LEAST-SQUARES FILTER 00011900 C ==================================================================== 00012000 C 00012100 LWAV2=2*LWAV 00012200 IF(LWAV2.LT.LFLTR) LWAV2 = LFLTR 00012300 LWH = MAX0((LWAV2-LWAV)/2,1) 00012400 C 00012500 CALL ARSET(SCR,LWAV2,0.0) 00012600 CALL ARMVE(WAVLET,SCR(LWH),LWAV) 00012700 C 00012800 AMP1 = 1.E70 00012900 AMP2 = -AMP1 00013000 C 00013100 DO 10 I=1,LWAV2 00013200 IF(AMP1.GT.SCR(I)) AMP1=SCR(I) 00013300 IF(AMP2.LT.SCR(I)) AMP2=SCR(I) 00013400 IF(SCR(I).EQ.AMP1) LSLE=I 00013500 IF(SCR(I).EQ.AMP2) LSHE=I 00013600 10 CONTINUE 00013700 C 00013800 SAC =ABS(SCR( LSLE)) 00013900 C IF (SAC.GT.SCR( LSHE)) WRITE ( IPR, 9000 ) 00014000 IF (SAC.GT.SCR( LSHE)) LSHE=LSLE 00014100 LPT = LSHE 00014200 C 00014300 C FILTER WAVELET AND DESIGN ERROR FUNCTION 00014400 C 00014500 IF ( NOFILT.EQ.0 ) 00014600 *CALL SAGFLT(LWAV2,SCR,FL,FH,DT,S1,ICDPN,ITRCN) 00014700 C 00014800 SUMSQ1 = 0.D00 00014900 DO 20 I=1,LWAV2 00015000 SUMSQ1 = SUMSQ1 + SCR(I)*SCR(I) 00015100 20 CONTINUE 00015200 C 00015300 CALL SAGAC3 (SCR,LWAV2,WORK,LFLTR, S1,S2) 00015400 C 00015500 WORK(1) = (1.+XN)*WORK(1) 00015600 WORK(LFLTR+1)=1.0 00015700 CALL ARSET (WORK(LFLTR+2),LFLTR-1,0.0) 00015800 CALL M1EURK (WORK,WORK(LFLTR+1),DCNFLT,LFLTR,S1,S2,1,ICDPN,ITRCN) 00015900 C 00016000 LFL2=LWAV2+LFLTR-1 00016100 LSH=LFL2-LPT+1 00016200 C 00016300 CALL SAGACC( SCR,LWAV2, DCNFLT,LFLTR, SCR(LWAV2+1),LFL2,1, S1,S2 )00016400 C 00016500 AMP1 = 1.E70 00016600 AMP2 = -AMP1 00016700 C 00016800 DO 30 I=1,LFL2 00016900 IF(AMP1.GT.SCR(LWAV2+I)) AMP1=SCR(LWAV2+I) 00017000 IF(AMP2.LT.SCR(LWAV2+I)) AMP2=SCR(LWAV2+I) 00017100 IF(SCR(LWAV2+I).EQ.AMP1) LSLE=I 00017200 IF(SCR(LWAV2+I).EQ.AMP2) LSHE=I 00017300 30 CONTINUE 00017400 C 00017500 SACX=ABS(SCR(LWAV2+LSLE)) 00017600 IF(SACX.LT.ABS(SCR(LWAV2+LSHE))) SACX=ABS(SCR(LWAV2+LSHE)) 00017700 C 00017800 CALL ARDVFC(SCR(LWAV2+1),SCR(LWAV2+1),SACX,LFL2) 00017900 C 00018000 C PLOT INPUT WAVELET AND MINIMUM-PHASE FILTER 00018100 C 00018200 ST = (1-LPT)*DT 00018300 C 00018400 WRITE ( IPR, 9010) ICDP 00018500 CALL USFPLT (SCR,LWAV2,ST,DT,0,IPR) 00018600 C 00018700 WRITE ( IPR, 9020) 00018800 CALL USFPLT (DCNFLT,LFLTR,0.,DT,0,IPR) 00018900 C 00019000 C DESIGN CORRECTION FILTER AND NORMALIZE FILTER 00019100 C 00019200 CALL SAGAC3 ( SCR(LWAV2+1),LFL2, WORK,LFL2, S1,S2 ) 00019300 C 00019400 WORK(1) = 1000.*WORK(1) 00019500 CALL ARMVE(SCR(LWAV2+1),WORK(LFL2+1),LFL2) 00019600 CALL ARFLP(WORK(LFL2+1),WORK(LFL2+1),LFL2) 00019700 CALL M1EURK (WORK,WORK(LFL2+1),INVFLT,LFL2,S1,S2,0,ICDPN,ITRCN) 00019800 C 00019900 AMP1 = 1.E70 00020000 AMP2 = -AMP1 00020100 C 00020200 DO 40 I=1,LFL2 00020300 IF(AMP1.GT.INVFLT(I)) AMP1=INVFLT(I) 00020400 IF(AMP2.LT.INVFLT(I)) AMP2=INVFLT(I) 00020500 IF(INVFLT(I).EQ.AMP1) LSLE=I 00020600 IF(INVFLT(I).EQ.AMP2) LSHE=I 00020700 40 CONTINUE 00020800 C 00020900 SAC =ABS (INVFLT(LSLE) ) 00021000 IF (SAC.GT.INVFLT(LSHE)) LSHE=LSLE 00021100 C 00021200 CALL ARDVFC(INVFLT,INVFLT,ABS(INVFLT(LSHE)),LFL2) 00021300 C 00021400 LFL4 = 2*LFL2 - 1 00021500 C 00021600 CC CALL SAGACC(INVFLT,LFL2, DCNFLT,LFLTR, WORK,LFL3, 1,S1,S2) 00021700 CC CALL SAGACC(WORK, LFL3, SCR ,LWAV2, WORK,LFL4, 1,S1,S2) 00021800 CALL ARMPFC(SCR(LWAV2+1),SCR(LWAV2+1),SACX,LFL2 ) 00021900 CALL SAGACC(SCR(LWAV2+1),LFL2, INVFLT,LFL2, WORK,LFL4, 1,S1,S2) 00022000 C 00022100 AMP1 = 1.E70 00022200 AMP2 = -AMP1 00022300 C 00022400 DO 50 I=1,LFL4 00022500 IF(AMP1.GT.WORK (I)) AMP1=WORK (I) 00022600 IF(AMP2.LT.WORK (I)) AMP2=WORK (I) 00022700 IF(WORK (I).EQ.AMP1) LSLE=I 00022800 IF(WORK (I).EQ.AMP2) LSHE=I 00022900 50 CONTINUE 00023000 C 00023100 SAC =ABS (WORK (LSLE) ) 00023200 IF (SAC.GT.WORK (LSHE)) LSHE=LSLE 00023300 C 00023400 I = LWAV2/2 - 1 00023500 I1 = MAX0(LSHE-I, 1) 00023600 I2 = I1 + LWAV2 - 1 00023700 SUMSQ2 = 0.D00 00023800 C 00023900 DO 60 I=I1,I2 00024000 SUMSQ2 = SUMSQ2 + WORK(I)*WORK(I) 00024100 60 CONTINUE 00024200 C 00024300 SCF = SNGL(DSQRT(SUMSQ1/SUMSQ2)) 00024400 C 00024500 CALL ARMPFC(INVFLT,INVFLT,SCF,LFL2 ) 00024600 C CALL ARMPFC( WORK , WORK ,SCF,LFL4 ) 00024700 C 00024800 WRITE ( IPR, 9030) 00024900 CALL USFPLT (WORK(LSH),2*LPT-1,ST,DT,0,IPR) 00025000 C 00025100 70 IF(LWAV.NE.0) GO TO 80 00025200 C 00025300 LFLT=LFLTR 00025400 CALL ARMVE (DATFLT,OUTFLT,LFLTR) 00025500 RETURN 00025600 C 00025700 80 LFL3=LFL2+LFLTR-1 00025800 C 00025900 CALL SAGACC( DATFLT,LFLTR, INVFLT,LFL2, OUTFLT,LFL3,1, S1,S2 ) 00026000 IF (NOFILT.EQ.0) 00026100 *CALL SAGFLT(LFL3,OUTFLT,FL,FH,DT,S1,ICDP,ITRC) 00026200 IF(IER.NE.0) GO TO 90 00026300 C 00026400 LFLT=LFL3 00026500 RETURN 00026600 C 00026700 90 IPRNT = -999 00026800 RETURN 00026900 C 00027000 C9000 FORMAT(//,' ### WARNING ### MAXIMUM AMPLITUDE IN WAVELET CORRESPO00027100 C *NDS TO TROUGH RATHER THAN PEAK ', /, 00027200 C * ' A POLARITY REVERSAL OR TIME SHIFT MAY00027300 C * OCCUR IN DECONVOLVED TRACES ' ) 00027400 C 00027500 9010 FORMAT(1H1,//, 10X, 'FILTERED INPUT WAVELET TO BE USED', 00027600 * /, 10X, ' STARTING WITH SHOT/CDP ',I5 , 00027700 * /, 10X, '---------------------------------',// ) 00027800 C 00027900 9020 FORMAT(1H1,//, 10X, 'MINIMUM PHASE FILTER FOR CURRENT WAVELET', 00028000 * /, 10X, '------- ----- ------ --- ------- -------',//)00028100 C 00028200 9030 FORMAT(1H1,//, 10X, 'FINAL (PHASE CORRECTED) DECONVOLUTION ', 00028300 * /, 10X, ' RESULT FOR CURRENT WAVELET ', 00028400 * /, 10X, '--------------------------------------',//) 00028500 C 00028600 END 00028700