CTITLESATCORF -- LAG, WEIGHT, AND CORRELATION OF TRACE AND PILOT C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR DANIEL POLAK CA DESIGNER DANIEL POLAK CA LANGUAGE VS FORTRAN CA SYSTEM IBM CA WRITTEN 10-31-80 C REVISED MO-DA-YR BY PROGRAMMER FOR REASON C REVISED 07-13-81 POLAK - CHANGED THE DIMENSION OF 'CCW' ARRAY C REVISED 10-05-82 NELAN - ADDED IN COMMON/APSTAT/ FOR AP C STATUS AND CHECK FOR RETRY. C REVISED 01-06-89 TRULOCK - CHANGE ARRAY PROCESSOR CALLS TO C FORTRAN CODE. C REVISED 10-24-89 TRULOCK - RENAME SUBROUTINE TO SATCORF. C REVISED 01-23-91 CLJ - ADDED SECONDARY PEAK OPTION CA CA CALL SATCORF (BULK, WORK, NSMPLE, FFTLN, LAG, PRLAG, NXCOR, CA * MXNXCR, CORFN, SHFT, AMP, DOT, SPEAK) CA CA IN/OUT ARGUMENT TYPE DESCRIPTION CA CA IN BULK R4 PILOT TRACE AND 'NXCOR' TRACES FOR CA CORRELATION WITH ZEROES ADDED FOR FFT CA IN WORK R4 WORK ARRAY FOR FFT'S CA IN NSMPLE I4 LENGTH OF EACH TRACE WITHOUT ADDED ZEROES CA IN FFTLN I4 LENGTH OF FFT CA IN LAG I4 NUMBER OF SAMPLES OF LAG TO EXAMINE IN CA EACH DIRECTION (ASSUMED TO BE POSITIVE) CA FOR THE SECONDARY PEAK CA IN PRLAG I4 NUMBER OF SAMPLES OF LAG TO EXAMINE IN CA EACH DIRECTION (ASSUMED TO BE POSITIVE) CA FOR THE PRIMARY PEAK CA IN NXCOR I4 NUMBER OF TRACES TO CORRELATE WITH PILOT CA IN MXNXCR I4 MAXIMUM NUMBER OF CORRELATIONS PER CALL CA IN CORFN R4 WORK AREA FOR CORRELATION FUNCTIONS CA OUT SHFT R4 COMPUTED TIME SHIFTS IN SAMPLES CA OUT AMP R4 AMPLITUDE VALUE OF TRACES CA OUT DOT R4 DOT PRODUCT OF TRACES AND PILOT FOR OPTIMAL CA LAG CA IN SPEAK I4 COMPUTED TIME IN SAMPLES FROM THE PRIMARY CA TO THE SECONDARY PEAK CA CA THIS ROUTINE COMPUTES LAGS, AMPLITUDES, AND CORRELATIONS OF CA 'NXCOR' INPUT TRACES RELATIVE TO PILOT TRACE. CA *** FREQUENCY DOMAIN CORRELATION *** CA CA C DEBUG UNIT(6),TRACE,INIT(R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12, C * R13) C END DEBUG SUBROUTINE SATCORF (BULK, WORK, NSMPLE, FFTLN, LAG, PRLAG, NXCOR, * MXNXCR, CORFN, SHFT, AMP, DOT, SPEAK) C IMPLICIT INTEGER (A-Z) EXTERNAL S1ATP C C REAL ARRAYS IN PARAMETER LIST C REAL AMP (MXNXCR) REAL CORFN ( 1) REAL DOT (MXNXCR) REAL SPEAK (MXNXCR) REAL SHFT (MXNXCR) REAL BULK ( 1) REAL WORK ( 1) C EQUIVALENCE (TRCNDX , R1) EQUIVALENCE (FFTNDX , R2) EQUIVALENCE (AMPNDX , R3) EQUIVALENCE (FFTLEN , R4) EQUIVALENCE (LENFFT , R5) EQUIVALENCE (NSMPL , R6) EQUIVALENCE (NCOR , R7) EQUIVALENCE (CORNDX , R9) EQUIVALENCE (CORLEN , R10) C C INTEGER ARRAYS -- LOCAL C INTEGER KY ( 2) C C REAL VARIABLES/ARRAYS -- LOCAL C REAL A REAL B REAL FACT REAL SCALE REAL T REAL XTEMP REAL XVPSS REAL SPCOR REAL SPCORL REAL SPCORR C C INITIALIZE C DO 10 I = 1, MXNXCR AMP(I) = 0. DOT(I) = 0. SHFT(I) = 0. SPEAK(I)= 0. 10 CONTINUE= 0. C C CALCULATE INDICES AND LENGTHS FOR THE INDEX REGISTERS (POINTERS) C TRCNDX = 0 NSMPL = NSMPLE FFTLEN = FFTLN + 2 FFTNDX = (NXCOR + 1) * FFTLEN LENFFT = FFTLEN / 2 AMPNDX = (NXCOR + 2) * FFTLEN FACT = 1.0 / (4.0 * FFTLN) CORNDX = NSMPL - 1 CORLEN = 2 * LAG + 1 NCOR = NXCOR C PCORLN = 2 * PRLAG + 1 IF (PCORLN .GT. CORLEN) PCORLN = CORLEN C C C INITIALIZE WORK REGISTERS (POINTERS) C R11 = R4 R11 = R11 - R10 R13 = R7 C C EXECUTE FFT ON THE PILOT TRACE C NXARRV = 0 NICTX = -2 + R4 C **NOTE** NICTX MUST EQUAL 2**M FOR SOME POSITIVE INTEGER M CALL ARMVE(BULK(NXARRV), BULK(R2), NICTX) CALL S1FMAG(NICTX, MAG, LFOUR) CALL S2DFT2(MAG, BULK(R2), *9010) SCALE = SQRT(NICTX*2.0) DO 1000 IVPSS = R2, R2+R4-1, 2 1000 BULK(IVPSS) = SCALE * BULK(IVPSS) SCALE = -SCALE DO 2000 IVPSS = R2+1, R2+R4-1, 2 2000 BULK(IVPSS) = SCALE * BULK(IVPSS) DO 3000 IVPSS = 0, R4-1 3000 BULK(IVPSS) = BULK(IVPSS + R2) C C TOP OF LOOP FOR CORRELATION OF 'MXNXCR' TRACES C 4000 CONTINUE C C UPDATE INDICES C R1 = R1 + R4 R3 = R3 + 1 C C COMPUTE THE SUM OF THE SQUARES OF THE GATHER TRACE ELEMENTS C XVPSS = 0.0 DO 5000 IVPSS = R1, R1+R4-1 5000 XVPSS = XVPSS + BULK(IVPSS) * BULK(IVPSS) BULK(R3) = XVPSS C C REVERSE THE GATHER TRACE FOR CORRELATION C IR6 = R6+R1-1 DO 6000 IVPSS = R1, R1+R6/2-1 XTEMP = BULK(IVPSS) BULK(IVPSS) = BULK(IR6) BULK(IR6) = XTEMP 6000 IR6 = IR6 - 1 C C EXECUTE FFT ON THE GATHER TRACE C NICTX = -2 + R4 C **NOTE** NICTX MUST EQUAL 2**M FOR SOME POSITIVE INTEGER M CALL ARMVE(BULK(R1), BULK(R2), NICTX) CALL S1FMAG(NICTX, MAG, LFOUR) CALL S2DFT2(MAG, BULK(R2), *9010) SCALE = SQRT(NICTX*2.0) DO 7000 IVPSS = R2, R2+R4-1, 2 7000 BULK(IVPSS) = SCALE * BULK(IVPSS) SCALE = -SCALE DO 7100 IVPSS = R2+1, R2+R4-1, 2 7100 BULK(IVPSS) = SCALE * BULK(IVPSS) C C MULTIPLY FFT OUTPUTS (COMPLEX MULTIPLY WITH CONJUGATE OUTPUT) C NXARRV = 0 CALL ARMPC (BULK(NXARRV), BULK(R2), BULK(R2), R5) C C EXECUTE AN INVERSE FFT ON THE CONJUGATE PRODUCT C NICTY = -2 + R4 C **NOTE** NICTX MUST EQUAL NICTY/2 + 1 DO 8000 IVPSS = 1, R4, 2 8000 WORK(IVPSS) = BULK(R2 + IVPSS - 1) DO 8100 IVPSS = 2, R4, 2 8100 WORK(IVPSS) = -BULK(R2 + IVPSS - 1) CALL S1FMAG(NICTY, MAG, LFOUR) CALL S2DFI2(MAG, WORK, *9010) SCALE = SQRT(NICTY*2.0) DO 8200 IVPSS = 1, NICTY 8200 BULK(R1 + IVPSS - 1) = SCALE * WORK(IVPSS) C C NORMALIZE THE OUTPUT C DO 9000 IVPSS = R1, R1+R4-1 9000 BULK(IVPSS) = BULK(IVPSS) * FACT C C MOVE CORRELATION FUNCTION FOR COMPRESSED RETURN OF ALL C CORRELATION FUNCTIONS C R11 = R11 + R10 R12 = R1 R12 = R12 + R9 C DO 10000 IVPSS = 1, R10 10000 BULK(IVPSS - 1 + R11) = BULK(IVPSS - 1 + R12) C C HAS 'NXCOR' TRACES BEEN CORRELATED? C R7 = R7 - 1 IF(R7 .GT. 0 ) GO TO 4000 C C PUT RESULTS INTO OUTPUT ARRAYS C NGET = MXNXCR * (2 * LAG + 1) C C MOVE CORRELATION FUNCTIONS INTO CORFN C CALL ARMVE (BULK(R4), CORFN, NGET) R3 = R3 - R13 R3 = R3 + 1 C DO 11000 IVPSS = 1, R13 11000 BULK(IVPSS - 1) = SQRT(ABS(BULK(IVPSS - 1 + R3))) C C MOVE AMPLITUDE DATA INTO AMP C NXARRV = 0 CALL ARMVE (BULK(NXARRV), AMP, MXNXCR) C C DETERMINE THE OPTIMAL RELATIVE LAGS C C POFSET = (CORLEN - PCORLN) / 2 C DO 24000 I = 1, NXCOR BEG = (I - 1) * CORLEN + 1 END = BEG + CORLEN - 1 C C LOCATE MAX CORRELATION SAMPLE C CCLJCCCCCCDO 12000 J = BEG, END C BEGPR = BEG + POFSET ENDPR = END - POFSET DO 12000 J = BEGPR, ENDPR IF (CORFN(J) .LT. DOT(I)) GO TO 12000 JJ = J DOT (I) = CORFN(J) 12000 CONTINUE C C IF CORRELATION NOT ZERO PERFORM QUADRATIC INTERPOLATION AROUND C PEAK SAMPLE TO OBTAIN PEAK CORRELATION AND TIME SHIFT. C IIPR = JJ - BEGPR + 1 II = JJ - BEG + 1 IF (DOT(I) .EQ. 0.) GO TO 14000 SHFT(I) = II - LAG - 1 CCLJ IF (II .EQ. 1 .OR. II .EQ. CORLEN) GO TO 13000 IF (IIPR .EQ. 1 .OR. IIPR .EQ. PCORLN) GO TO 13000 A = .5 * (CORFN(JJ+1) - 2. * DOT(I) + CORFN(JJ-1)) IF (A .EQ. 0.) GO TO 14000 B = .5 * (CORFN(JJ+1) - CORFN(JJ-1)) T = -.5 * B / A DOT(I) = (A * T + B) * T + DOT(I) SHFT(I) = SHFT(I) + T GO TO 14000 C 13000 DOT(I) = 0. C 14000 CONTINUE C C C LOCATE MAX CORRELATION SAMPLE ON SECONDARY PEAK C IF (DOT(I) .EQ. 0.0) GO TO 17000 C C FIND SECONDARY PEAK ON THE LEFT SIDE C KY(1) = 0 IF (JJ-BEG+1 .GT. 1) * CALL SIZC (CORFN(BEG), -1, JJ-BEG+1, 2, KY(1) ) IF (KY(1) .EQ. 0) THEN LTM = 0 SPCORL = 0.0 GO TO 13220 ELSE IFZC = KY(1)-1 + BEG END IF CALL SIZC (CORFN(BEG), -1, JJ-BEG+1, 3, KY(1) ) IF (KY(1) .EQ. 0) THEN LTM = 0 SPCORL = 0.0 GO TO 13220 ELSE ILZC = KY(1) + BEG END IF C LTM = 0 SPCORL = 0.0 DO 13200 J = IFZC, ILZC, -1 IF (CORFN(J) .GT. SPCORL) THEN LTM = J SPCORL = CORFN(J) END IF 13200 CONTINUE C 13220 CONTINUE C C FIND SECONDARY PEAK ON THE RIGHT SIDE C KY(1) = 0 IF (END-JJ+1 .GT. 1) * CALL SIZC (CORFN(JJ), 1, END-JJ+1, 2, KY(1) ) IF (KY(1) .EQ. 0) THEN RTM = 0 SPCORR = 0.0 GO TO 13420 ELSE IFZC = KY(1)-1 + JJ END IF CALL SIZC (CORFN(JJ), 1, END-JJ+1, 3, KY(1) ) IF (KY(1) .EQ. 0) THEN RTM = 0 SPCORR = 0.0 GO TO 13420 ELSE ILZC = KY(1)-2 + JJ END IF C RTM = 0 SPCORR = 0.0 DO 13400 J = IFZC, ILZC IF (CORFN(J) .GT. SPCORR) THEN RTM = J SPCORR = CORFN(J) END IF 13400 CONTINUE C 13420 CONTINUE C C PICK THE LARGER OF THE TWO CORRELATIONS IF (SPCORL .GT. SPCORR) THEN SPCOR = SPCORL SPTM = LTM ELSE SPCOR = SPCORR SPTM = RTM END IF C C C INTERPOLATION FOR SECONDARY PEAK C IF CORRELATION NOT ZERO PERFORM QUADRATIC C INTERPOLATION AROUND PEAK SAMPLE TO OBTAIN C PEAK CORRELATION AND TIME SHIFT. C II = SPTM - BEG + 1 IF (SPCOR .EQ. 0) GO TO 17000 SPEAK(I) = II - LAG - 1 IF (II .EQ. 1 .OR. II .EQ. CORLEN) GO TO 17000 A = .5 * (CORFN(SPTM+1) - 2. * SPEAK(I) + CORFN(SPTM-1)) IF (A .EQ. 0.) GO TO 18000 B = .5 * (CORFN(SPTM+1) - CORFN(SPTM-1)) T = -.5 * B / A SPEAK(I) = SPEAK(I) + T SPEAK(I) = ABS (SPEAK(I) - SHFT(I)) GO TO 18000 C 17000 CONTINUE SPEAK(I) = 0.0 18000 CONTINUE C C 24000 CONTINUE C RETURN 9010 PRINT *,' SATCORF ERROR 9010' RETURN END