CTITLESARKFT - DESIGN NMO DESTRETCH INVERSION FILTERS FOR RICKER WAVELET C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** C @PROCESS VECTOR(LEVEL(2) REPORT(XLIST)) C CA AUTHOR H. W. SWAN CA DESIGNER H. W. SWAN CA SYSTEM IBM/CRAY CA LANGUAGE VS FORTRAN VERSION 2.2 CA WRITTEN 12-03-87 CA LAST REVISED 07-19-88 CA C REVISED HWS 04-26-88 CHANGED SWLEV CALL TO OPFILT C HWS 07-19-88 FIXED ERROR IN H1, H2 INITIALIZATION C REVISED JJC 05-23-91 CHANGED IMPLICIT TO (A-Z). C REVISED JJC 12-17-91 MODIFIED TO MEET SPARC STANDARDS. CA CA CALLING SEQUENCE: CA CALL SARKFT(H1, H2, CORR, CVECT, WORK, MXH, NH, CA W0T, PCTWNS, PCTCNS, A, KS, IER) CA CA IN/OUT ARGUMENT TYPE DESCRIPTION CA ------ -------- ---- ----------- CA CA IN MXH I4 THE DIMENSIONED SIZE OF THE OUTPUT INVERSION CA FILTERS CA CA IN NH I4 THE ACTUAL SIZE OF THE OUTPUT INVERSION CA FILTERS CA CA IN W0T R4 OMEGA 0 TIMES THE TIME STEP SIZE (UNITLESS) CA CA IN PCTWNS R4 THE ADDITIONAL WHITE NOISE PERCENTAGE CA CA IN PCTCNS R4 THE ADDITIONAL COLORED NOISE PERCENTAGE CA CA IN A R4 THE AUTOREGRESSIVE NOISE COEFFICIENT. CA RN(N) = PCTNS * 1E-2 * RW(0) * A ** N CA CA IN KS I4 SECOND DIMENSION OF THE WORK ARRAY CA CA OUT H1, H2 (R4) THE INVERSION FILTERS. CA DIMENSIONED: (-MXH:MXH) CA CA OUT CORR (R4) THE AUTOCORRELTATION OF THE WAVELET. CA 1-D ARRAY DIMENSIONED: (0:2*MAXH) CA CA OUT CVECT (R4) THE CROSS CORRELATION VECTOR BETWEEN CA W(T) AND TW'(T). DIMENSIONED: CA (-MAXH:MAXH) CA CA OUT WORK (R8) A WORK ARRAY. 2-D ARRAY, CA DIMENSIONED: (-MXH:MXH,KS), CA WHERE "KS" IS 4 CA FOR THE IBM MAINFRAME VERSION. CA KS=1: THE CROSS CORRELATION BETWEEN W(T) AND CA TW'(T) CA KS>1: A WORK AREA FOR SUBROUTINE SWLEV. CA CA WARNINGS: 1) THIS ARRAY MUST BE ALIGNED ON A DOUBLE WORD BOUNDARY. CA 2) SINCE THIS ARRAY IS PASSED FROM SPAVEL AS 'XCOM(WORK)' CA IT IS DECLARED TO BE SINGLE PRECISION WITHIN SARKFT CA IN ORDER TO PREVENT AN INTER-COMPILATION CONFLICT. CA THIS IS ACCEPTABLE, INDIVIDUAL ELEMENTS OF 'WORK' ARE CA NEVER EXPLICITLY REFERENCED, ONLY PASSED TO 'SWLEV'. CA BEWARE THAT IT SHOULD NEVER BE EXPLICITLY REFERENCED! CA 3) ON THE CRAY, THIS ARRAY IS USED AS A SINGLE-PRECISION CA WORK ARRAY, THROUGH SCILIB ROUTINE OPFILT. CA CA CA CA OUT IER I4 AN ERROR FLAG: 0 = NORMAL COMPLETION, CA 1 = AUTOCORRELATION MATRIX WAS SINGULAR. CA 2 = MXH WAS NOT LARGE ENOUGH. CA CA CA PURPOSE: CA CA THIS SUBROUTINE GENERATES THE H1(T) AND H2(T) INVERSION CA FILTERS FOR NMO DESTRETCHING, SPECIFICALLY FOR A RICKER WAVELET. CA CA CA CA METHOD: CA CA USING THE FUNCTIONAL FORM FOR A RICKER WAVELET: CA CA 2 2 CA W(T) = (1 - 0.5*W0T ) * EXP (-0.25*W0T ) , CA CA THE AUTOCORRELATION OF W(T) AND CROSS CORRELATION BETWEEN CA W(T) AND {TW'(T)} ARE EXPLICITLY EVALUATED. THE INVERSION CA FILTERS ARE THEN OBTAINED THROUGH CONVENTIONAL LEVINSON CA INVERSIONS. CA CA RESTRICTIONS: (CHECKED BY IER FLAG) CA CA 1) PCTNS > 0.0 CA 2) NH <= MXH CA CA CA SUBROUTINES CALLED: CA CA SCOPY - ESSL/SCILIB ROUTINE TO COPY A VECTOR CA OPFILT - SCILIB ROUTINE TO PERFORM A LEVINSON INVERSION. CA CA NOTE: AN IBM INTERFACE TO ESSL ROUTINE 'SWLEV' IS PART OF CA PRODUCTION SPARC. CA SUBROUTINE SARKFT(H1, H2, CORR, CVECT, WORK, MXH, NH, * W0T, PCTWNS, PCTCNS, A, KS, IER) CJJ IMPLICIT NONE IMPLICIT INTEGER (A-Z) C REAL THRESH REAL C0, C1, CORR, CVECT, H1, H2, PCTCNS, PCTWNS REAL T, T2, TEMP, W0T, WORK INTEGER I, IER, KS, LAG, MAXLAG, MXH, NH, NWORK DOUBLE PRECISION A, AN DIMENSION CORR(0:2*MXH), WORK(-MXH:MXH,KS) DIMENSION H1(-MXH:MXH), H2(-MXH:MXH), CVECT(-MXH:MXH) PARAMETER(THRESH=200.0) EXTERNAL SCOPY, OPFILT C C CHECK THE ARRAY SIZES. C IER = 0 MAXLAG = 2*NH NWORK = KS * (MAXLAG + 1) IF(PCTWNS .LE. 0.0) IER = 1 IF(NH .GT. MXH) IER = 2 IF(ABS(A) .GE. 1.0) IER = 3 IF(IER .NE. 0) RETURN C C ZERO OUT THE FULL H1 AND H2 ARRAYS. C DO 100 LAG = -MXH, MXH H1(LAG) = 0.0 100 H2(LAG) = 0.0 C C COMPUTE THE WAVELET'S AUTOCORRELATION FUNCTION. C DO 200 LAG = 0, MAXLAG T = LAG * W0T T2 = T*T CORR(LAG) = 0.0 IF(T2 .LT. THRESH) THEN CORR(LAG) = (48 - (24 - T2)*T2) * EXP(-T2/8) ENDIF 200 CONTINUE C C COMPUTE THE H1(T) INVERSION FILTER. C CALL SCOPY(NH+1, CORR, 1, CVECT(0), 1) CALL SCOPY(NH+1, CORR, 1, CVECT(-NH), -1) C0 = CORR(0) * PCTWNS * 1E-2 C1 = CORR(0) * PCTCNS * 1E-2 CORR(0) = CORR(0) + C0 AN = C1 DO 350 I=0, MAXLAG CORR(I) = CORR(I) + AN 350 AN = AN * A C C CALL SWLEV(CORR, 1, CVECT(-NH), 1, H1(-NH), 1, C * MAXLAG+1, WORK, NWORK) C CALL OPFILT(MAXLAG+1, H1(-NH), CVECT(-NH), WORK, CORR) C C COMPUTE THE CROSS CORRELATION BETWEEN W(T) AND TW'(T). C DO 500 LAG = -NH, NH T = LAG * W0T T2 = T*T TEMP = 0.0 IF(T2 .LE. THRESH) THEN TEMP = -(((T2 - 36)*T2 + 144)*T2 + 192) * * EXP(-T2/8)/8 ENDIF 500 CVECT(LAG) = TEMP C C COMPUTE THE H2(T) FILTER. C C CALL SWLEV(CORR, 1, CVECT(-NH), 1, H2(-NH), 1, C * MAXLAG+1, WORK, NWORK) C CALL OPFILT(MAXLAG+1, H2(-NH), CVECT(-NH), WORK, CORR) AN = C1 CORR(0) = CORR(0) - C0 DO 600 I=0, MAXLAG CORR(I) = CORR(I) - AN 600 AN = AN * A RETURN END