CTITLESAWEST -- CALCULATE AVERAGE OPERATOR AND ESTIMATE WAVELET 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D. NYMAN 00000200 CA DESIGNER D. NYMAN 00000300 CA LANGUAGE S/370 FORTRAN H 00000400 CA SYSTEM IBM AND CRAY CA WRITTEN 05/19/78 00000500 C REVISED 10/15/86 REM. FIX TO ASSURE R*8 BOUNDARY FOR ARY(WAV) 00000600 C REVISED 10/15/90 MAA. MODIFIED FOR CRAY USE. CA 00000700 CA 00000800 CA CALL SAWEST 00000900 CA (PCC, PAR, NOPRS, SMPINT, WAVLEN, STABLZ, ARY, NUSM, IER) 00001000 CA 00001100 CA INPUT PCC = ARRAY OF SUMS OF PARTIAL CORRELATION I2,R8 00001200 CA COEFFICIENTS OF ALL TRACES DECONED, FOR 00001300 CA EACH WINDOW. FIRST 8 BYTES ARE TRACE 00001400 CA COUNTERS FOR UP TO FOUR WINDOWS 00001500 CA PAR(1,I) = I-TH OPERATOR LENGTH I4 00001600 CA PAR(2,I) = I-TH PREDICTION DISTANCE (>0), OR NUMBER 00001700 CA OF ANTICIPATORY TERMS IN OPERATOR (<0) 00001800 CA NOPRS = NUMBER OF OPERATORS AND WINDOWS I4 00001900 CA SMPINT = SAMPLING INTERVAL IN MS I4 00002000 CA WAVLEN = LENGTH OF ESTIMATED WAVELET IN MS I4 00002100 CA STABLZ = STABILIZATION FACTORS, FOR INTERPOLATION I4 00002200 CA DECON WAVELETS ONLY 00002300 CA NUSM = TOTAL AVAILABLE LENGTH OF ARRAY ARY I4 00002400 CA 00002500 CA OUTPUT PCC = AVERAGE PARTIAL CORRELATION COEFFICIENTS I2,I4 00002600 CA PAR(1,I) = START TIME IN MS OF EACH SUB-ARRAY I4 00002700 CA PAR(2,I) = LENGTH IN POINTS OF EACH SUB-ARRAY 00002800 CA PAR(3,I) = TIME IN MS OF FIRST OPERATOR TERM 00002900 CA PAR(4,I) = NUMBER OF TERMS IN OPERATOR 00003000 CA PAR(5,I) = NUMBER OF OPERATOR AND WAVELET SUB-ARRAYS 00003100 CA FOR EACH WINDOW 00003200 CA ARY = ARRAY OF OPERATOR(S) AND WAVELET(S) R4 00003300 CA NUSM = REMAINING UNUSED LENGTH OF ARRAY ARY I4 00003400 CA IER = 0: NO ERRORS I4 00003500 CA 1: ARRAY ARY TOO SHORT 00003600 CA 2: FFT LENGTH TOO LONG 00003700 CA 3: NO TRACES PROCESSED FOR ONE OR MORE WINDOWS 00003800 CA 4: NO WAVELET ESTIMATION FOR PREDICTION DECON 00003900 CA 00004000 CA FOR EACH TIME WINDOW, SAWEST CALCULATES THE AVERAGE OPERATOR, 00004100 CA AND 'INVERTS' THAT OPERATOR TO OBTAIN AN ESTIMATE OF THE SEISMIC 00004200 CA WAVELET. FOR EACH WINDOW, THE SAME OPERATOR PARAMETERS MUST BE 00004300 CA USED FOR ALL TRACES. 00004400 CA 00004500 CAEND 00004600 C=======================================================================00004700 C EJECT 00004800 C=======================================================================00004900 C 00005000 C 00005100 C LOCAL VARIABLES AND CONSTANTS 00005200 C 00005300 C FFTL = FFT LENGTH I4 00005400 C FFTX = BASE 2 LOGARITHM OF FFTL: FFTL = 2**FFTX I4 00005500 C IC = ARRAY ARY INDEX POINTER I4 00005600 C IPC = ARRAY PCC INDEX POINTER I4 00005700 C LENG = SUB-ARRAY LENGTH I4 00005800 C NANT = NUMBER OF ANTICIPATORY TERMS IN OPERATOR I4 00005900 C NAVE = NUMBER OF TRACES DECONED IN EACH WINDOW, FROM PCC(1) I2 00006000 C NMEM = NUMBER OF MEMORY TERMS IN OPERATOR I4 00006100 C OLEN = TOTAL LENGTH OF OPERATOR I4 00006200 C OPL = NUMBER OF OPERATOR COEFFICIENTS I4 00006300 C OPR = ARRAY ARY OPERATOR INDEX I4 00006400 C OST = OPERATOR START TIME I4 00006500 C PASS = INTERPOLATION DECON WAVELET STABILIZATION INDEX I4 00006600 C PRED = PREDICTION DISTANCE / ANTICIPATION LENGTH I4 00006700 C RAVE = EQUIVALENCED TO NAVE R8 00006800 C WAV = ARRAY ARY WAVELET INDEX I4 00006900 C WST = WAVELET START TIME I4 00007000 C 00007100 C=======================================================================00007200 C EJECT 00007300 C 00007400 SUBROUTINE SAWEST 00007500 * (PCC, PAR, NOPRS, SMPINT, WAVLEN, STABLZ, ARY, NUSM, IER) 00007600 C 00007700 IMPLICIT INTEGER (A-Z) 00007800 REAL ARY, X, Y 00007900 C REAL*8 PCC, Z, RAVE 00008000 DOUBLE PRECISION PCC, Z, RAVE 00008000 INTEGER NAVE(4) 00008100 EQUIVALENCE (RAVE,NAVE(1)) 00008200 DIMENSION PCC(1), PAR(5,1), STABLZ(1), ARY(1) 00008300 EXTERNAL S1ATP 00008400 C 00008500 RAVE=PCC(1) 00008600 C....................................................................... WRITE(6,1) RAVE, NAVE(1) 1 FORMAT('<* (SAWEST)BEFORE: RAVE;NAVE *> ', I15, I10) C....................................................................... C CALL USEQVL (RAVE, NAVE) C....................................................................... WRITE(6,2) RAVE, NAVE(1) 2 FORMAT('<* (SAWEST)AFTER : RAVE;NAVE *> ', I15, I10) C....................................................................... C 00010590 IPC=2 00008700 IC=1 00008800 C 00008900 C GENERATE AVERAGE OPERATOR AND WAVELET 00009000 C ARRAYS FOR EACH TIME WINDOW 00009100 C 00009200 DO 100 I=1,NOPRS 00009300 C 00009400 OPL=PAR(1,I) 00009500 PRED=PAR(2,I) 00009600 OLEN=OPL+MAX0(PRED,1) 00009700 OST=MIN0(PRED,0)*SMPINT 00009800 LENG=MAX0(OLEN,WAVLEN/SMPINT+1) 00009900 WST=((OST*LENG-OLEN/2)/OLEN-SMPINT/2)/SMPINT*SMPINT 00010000 PAR(1,I)=WST 00010100 PAR(2,I)=LENG 00010200 PAR(3,I)=OST 00010300 PAR(4,I)=OLEN 00010400 PAR(5,I)=0 00010500 OPR=IC+(OST-WST)/SMPINT 00010600 WAV=IC+LENG 00010700 DWAV = WAV 00010710 C 00010800 C ASSURE DOUBLE WORD BOUNDARY FOR ARY(WAV) 00010900 C 00011000 IF (MOD(LOC(ARY(DWAV)),8) .NE. 0) DWAV = DWAV + 1 00011100 IF (DWAV+LENG-1.GT.NUSM) GO TO 910 00011200 CALL ARSET (ARY(IC),LENG,0.) 00011300 C 00011400 C DIVIDE PCC SUM BY NUMBER OF TRACES PROCESSED 00011500 C 00011600 N=NAVE(I) 00011700 IF (N.EQ.0 .AND. OPL.GE.2) GO TO 930 00011800 X=1./AMAX0(N,1) 00011900 K=IPC+OPL-1 00012000 C 00012100 DO 5 J=IPC,K 00012200 5 PCC(J)=PCC(J)*X 00012300 C 00012400 IF (PRED.GT.1) GO TO 940 00012500 C 00012600 C SPIKING / INTERPOLATION DECON COEFFICIENTS 00012700 C 00012800 NANT=-MIN0(PRED,0) 00012900 NMEM=OPL-NANT 00013000 CALL SAKCOF (PCC(IPC),NANT,NMEM,ARY(OPR),ARY(DWAV)) 00013100 IPC=IPC+OPL 00013200 IF (OPL.LT.2) GO TO 100 00013300 C 00013400 IF (NANT.EQ.0) GO TO 60 00013500 C 00013600 C INTERPOLATION DECON WAVELET FROM 00013700 C FREQUENCY DOMAIN INVERSION 00013800 C 00013900 FFTX=ALOG(LENG*2-1.)/.693147+1. 00014000 FFTL=2**FFTX 00014100 PASS=1 00014200 C 00014300 C PERFORM INVERSION FOR EACH VALUE OF THE ARRAY 00014400 C STABLZ UNTIL A NEGATIVE VALUE IS ENCOUNTERED 00014500 C 00014600 10 IF (WAV+FFTL+1.GT.NUSM) GO TO 910 00014700 CALL ARSET (ARY(WAV),FFTL,0.) 00014800 K=FFTL-NANT-OPR 00014900 M=OPR+OLEN-1 00015000 X=0. 00015100 C 00015200 DO 20 J=OPR,M 00015300 L=MOD(J+K,FFTL) 00015400 ARY(WAV+L)=ARY(J) 00015500 20 X=ARY(J)*ARY(J)+X 00015600 C 00015700 C FOURIER TRANSFORM OPERATOR 00015800 C 00015900 CALL S2DFT2 (FFTX,ARY(WAV),*920) 00016000 C 00016100 C INVERT FOURIER TRANSFORM 00016200 C STABILIZE INVERSE IF STABLZ(PASS) > 0 00016300 C 00016400 X=10.**(STABLZ(PASS)*.5-10.5)*X/(FFTL*FFTL) 00016500 IF (STABLZ(PASS).EQ.0) X=0. 00016600 K=WAV+FFTL 00016700 C 00016800 DO 30 J=WAV,K,2 00016900 Y=ARY(J)*ARY(J)+ARY(J+1)*ARY(J+1) 00017000 IF (Y.LT.X) Y=SQRT(Y*X) 00017100 ARY(J)=ARY(J)/Y 00017200 30 ARY(J+1)=-ARY(J+1)/Y 00017300 ARY(WAV)=0. 00017400 C 00017500 C INVERSE FOURIER TRANSFORM 00017600 C 00017700 CALL S2DFI2 (FFTX,ARY(WAV),*920) 00017800 C 00017900 C REPOSITION WAVELET IN ARRAY AND NORMALIZE 00018000 C 00018100 K=WAV+LENG 00018200 L=LENG+WST/SMPINT+FFTL 00018300 X=0. 00018400 C 00018500 DO 40 J=1,LENG 00018600 M=MOD(L-J,FFTL) 00018700 ARY(K-J)=ARY(WAV+M) 00018800 40 X=AMAX1(ABS(ARY(K-J)),X) 00018900 C 00019000 X=1./X 00019100 K = K - 1 00019200 DO 50 J=WAV,K 00019300 50 ARY(J)=ARY(J)*X 00019400 C 00019500 PASS=PASS+1 00019600 PAR(5,I)=PASS 00019700 IF (STABLZ(PASS).LT.0) GO TO 100 00019800 WAV=WAV+LENG 00019900 GO TO 10 00020000 C 00020100 C SPIKING DECON WAVELET ESTIMATION 00020200 C FROM INVERSION OF MEMORY TYPE OPERATOR 00020300 C 00020400 60 ARY(WAV)=1. 00020500 N=LENG-1 00020600 M=OLEN-1 00020700 X=1. 00020800 C 00020900 DO 80 J=1,N 00021000 C 00021100 L=MIN0(J,M) 00021200 Z=0. 00021300 C 00021400 DO 70 K=1,L 00021500 70 Z=ARY(WAV+J-K)*ARY(OPR+K)+Z 00021600 C 00021700 ARY(WAV+J)=-Z 00021800 80 X=AMAX1(ABS(ARY(WAV+J)),X) 00021900 C 00022000 C NORMALIZE WAVELET 00022100 C 00022200 X=1./X 00022300 K=WAV+LENG-1 00022400 DO 90 J=WAV,K 00022500 90 ARY(J)=ARY(J)*X 00022600 C 00022700 PAR(5,I)=2 00022800 C 00022900 100 IC=IC+LENG*PAR(5,I) 00023000 C 00023100 C NORMAL TERMINATION 00023200 C 00023300 NUSM=NUSM-IC+1 00023400 IER=0 00023500 GO TO 1000 00023600 C 00023700 C INSUFFICIENT MEMORY 00023800 C 00023900 910 IER=1 00024000 GO TO 1000 00024100 C 00024200 C FFT LENGTH TOO LONG 00024300 C 00024400 920 IER=2 00024500 GO TO 1000 00024600 C 00024700 C NO TRACES PROCESSED 00024800 C 00024900 930 IER=3 00025000 GO TO 1000 00025100 C 00025200 C PREDICTION DECON OPERATOR NOT ALLOWED 00025300 C 00025400 940 IER=4 00025500 C 00025600 1000 RETURN 00025700 C 00025800 END 00025900