CTITLESAKDCN -- MINIMUM VARIANCE SPIKING & INTERPOLATION DECON 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D. NYMAN 00000020 CA DESIGNER D. NYMAN 00000030 CA LANGUAGE S/370 FORTRAN H 00000040 CA SYSTEM IBM AND CRAY CA WRITTEN 05-11-78 00000050 C REVISED 10-06-90 MAA. MODIFIED FOR CRAY-YMP USE. 00000060 C REVISED 00000070 CA 00000080 CA 00000090 CA CALL SAKDCN (X, NPTS, NANT, NMEM, VAR, PCC, F, B, D) 00000100 CA 00000110 CA INPUT X = TRACE TIME WINDOW ARRAY R4 00000120 CA NPTS = NUMBER OF POINTS IN ARRAY X I4 00000130 CA NANT = NUMBER OF ANTICIPATORY OPERATOR COEFFICIENTS I4 00000140 CA NMEM = NUMBER OF MEMORY OPERATOR COEFFICIENTS I4 00000150 CA 00000160 CA OUTPUT X = DECONVOLVED TRACE TIME WINDOW ARRAY R4 00000170 CA VAR = RELATIVE OUTPUT TRACE VARIANCE R4 00000180 CA PCC = PARTIAL CORRELATION COEFFICIENT ARRAY R8 00000190 CA 00000200 CA F, B, D = TEMPORARY WORK ARRAYS, AT LEAST OF LENGTH R4 00000210 CA NPTS+NANT+NMEM 00000220 CA 00000230 CA 00000240 CA SAKDCN PERFORMS MINIMUM VARIANCE SPIKING (NANT = 0) AND 00000250 CA INTERPOLATION (NANT > 0) DECON ON THE TRACE SEGMENT X. THE 00000260 CA DECON OPERATOR, WITH NANT ANTICIPATORY TERMS AND NMEM MEMORY 00000270 CA TERMS CAN BE OBTAINED FROM THE PARTIAL CORRELATION COEFFICIENTS 00000280 CA WITH THE ROUTINE 'SAKCOF'. 00000290 CA 00000300 CA 00000310 C=======================================================================00000320 C EJECT 00000330 C=======================================================================00000340 C 00000350 C THE MINIMUM VARIANCE SPIKING AND INTERPOLATION DECON METHOD 00000360 C IMPLEMENTED HERE IS AN EXTENTION OF THE GRAM-SCHMIDT 00000370 C AUTO-REGRESSIVE ORTHO-NORMAL DECONVOLUTION FORMALISM 00000380 C DEVELOPED BY A. ALAM (ARCO R&D, PLANO, TEXAS). 00000390 C 00000400 C LOCAL VARIABLES AND CONSTANTS 00000410 C 00000420 C B = DELAYED BACKWARD PREDICTION ERROR R4 00000430 C E = INTERPOLATION ERROR (EPSILON) R4 00000440 C F = FORWARD PREDICTION ERROR R4 00000450 C NP = TOTAL NUMBER OF POINTS IN X, F, B, E I4 00000460 C RD = RHO(DELTA) R4 00000470 C RE = RHO(EPSILON) R4 00000480 C RO = RHO R4 00000490 C SE = SQUARED NORM OF E R4 00000500 C SFB = SQUARED NORM OF F AND B R4 00000510 C SIG = SIGMA R4 00000520 C SX0 = INPUT TRACE VARIANCE R4 00000530 C 00000540 C AS IMPLEMENTED, THE STATIONARITY ASSUMPTION ABOUT THE DATA IS 00000550 C EXPLICITLY EXPLOITED: SINCE UNDER THIS ASSUMPTION THE INNER 00000560 C PRODUCT VALUES ARE STATISTICALLY INDEPENDENT OF THE PARTICULAR 00000570 C DATA SEGMENT, AS THE FORWARD AND BACKWARD PREDICTION ERRORS ARE 00000580 C SHIFTED AT EACH ITERATION THE OPERATION WINDOW IS SHORTENED TO 00000590 C INCLUDE ONLY THE OVERLAPING TERMS. 00000600 C 00000610 C DURING CALCULATION OF THE 'ANTICIPATION' PART OF INTERPOLATION 00000620 C DECON, THE INTERPOLATION ERROR IS ADVANCED ONE TIME UNIT AT 00000630 C EACH ITERATION IN ORDER TO RETAIN THE INDEXING CONVENTIONS OF 00000640 C THE MEMORY/SPIKING DECON CODE. 00000650 C 00000660 C=======================================================================00000670 C EJECT 00000680 C 00000690 SUBROUTINE SAKDCN (X, NPTS, NANT, NMEM, VAR, PCC, F, B, E) 00000700 C 00000710 DIMENSION X(1), PCC(1), F(1), B(1), E(1) 00000720 DOUBLE PRECISION PCC 00000730 C 00000740 C CHECK FOR ZERO VALUED DATA 00000750 C 00000760 L=NMEM+1 00000770 M=NPTS-NANT 00000780 SX=0. 00000790 C 00000800 DO 5 I=L,M 00000810 5 SX=X(I)*X(I)+SX 00000820 C 00000830 IF (SX.EQ.0. .OR. NANT+NMEM.LT.2) GO TO 110 00000840 C 00000850 C INITIALIZE F, B, SF, SB, AND RFB 00000860 C 00000870 SB=0. 00000880 RFB=0. 00000890 C 00000900 DO 10 I=2,NPTS 00000910 F(I)=X(I) 00000920 B(I)=X(I-1) 00000930 SB=B(I)*B(I)+SB 00000940 10 RFB=F(I)*B(I)+RFB 00000950 C 00000960 SF=SB-X(1)*X(1)+X(NPTS)*X(NPTS) 00000970 NC=NMEM 00000980 IF (NMEM.EQ.0) NC=NANT 00000990 NP=NPTS-1 00001000 N1=NPTS+1 00001010 C 00001020 C PERFORM SPIKING DECON OF ORDER 00001030 C NMEM (NANT IF NMEM = 0) 00001040 C 00001050 DO 20 J=1,NC 00001060 IF (SF*SB.EQ.0) GO TO 110 00001070 RF=-RFB/SB 00001080 RB=-RFB/SF 00001090 PCC(J)=(RF+RB)*.5 00001100 IF (DABS(PCC(J)).GT.1.) GO TO 110 00001110 SF=0. 00001120 SB=0. 00001130 RFB=0. 00001140 NP=NP-1 00001150 B(N1)=B(NPTS)+RB*F(NPTS) 00001160 C 00001170 C OBTIN F AND B OF ORDER N FROM F AND B OF ORDER 00001180 C N-1, AND CALCULATE NEW VALUES FOR SF, SB AND RFB 00001190 C 00001200 DO 20 I=1,NP 00001210 L=N1-I 00001220 F(L)=F(L)+RF*B(L) 00001230 B(L)=B(L-1)+RB*F(L-1) 00001240 SF=F(L)*F(L)+SF 00001250 SB=B(L)*B(L)+SB 00001260 20 RFB=F(L)*B(L)+RFB 00001270 C 00001280 L=NPTS-NP 00001290 F(L)=F(L)+RF*B(L) 00001300 C 00001310 IF (NANT*NMEM.EQ.0) GO TO 70 00001320 C 00001330 C PERFORM ANTICIPATION PART OF INTERPOLATION DECON 00001340 C 00001350 SIG=0. 00001360 X(NPTS+1)=0. 00001370 C 00001380 DO 30 I=L,NPTS 00001390 E(I)=F(I) 00001400 30 SIG=E(I)*X(I+1)+SIG 00001410 C 00001420 SE=SF+E(L)*E(L)-E(NPTS)*E(NPTS) 00001430 C 00001440 C AT EACH ITERATION ON J, FIRST F AND B OF ORDER 00001450 C (NMEM+J) ARE CALCULATED FROM F AND B OF ORDER 00001460 C (NMEM+J-1), AND THEN E OF ORDER (J,NMEN) IS 00001470 C FOUND FROM E OF ORDER (J-1,NMEM) AND F OF ORDER 00001480 C (NMEM+J). 00001490 C 00001500 DO 50 J=1,NANT 00001510 IF (SF*SB*SE.EQ.0) GO TO 110 00001520 RF=-RFB/SB 00001530 RB=-RFB/SF 00001540 PCC(NMEM+J)=(RF+RB)*.5 00001550 IF (DABS(PCC(NMEM+J)).GT.1.) GO TO 110 00001560 SF=(1.-RF*RF)*SF 00001570 A=SF+SIG*SIG/SE 00001580 RD=-SIG/A 00001590 RE=SF/A 00001600 SF=0. 00001610 SB=0. 00001620 SE=0. 00001630 RFB=0. 00001640 SIG=0. 00001650 NP=NP-1 00001660 C 00001670 DO 40 I=1,NP 00001680 L=N1-I 00001690 M=L-1 00001700 F(L)=F(L)+RF*B(L) 00001710 B(L)=B(M)+RB*F(M) 00001720 E(L)=RE*E(M)+RD*F(L) 00001730 SF=F(L)*F(L)+SF 00001740 SB=B(L)*B(L)+SB 00001750 SE=E(L)*E(L)+SE 00001760 RFB=F(L)*B(L)+RFB 00001770 40 SIG=E(L)*X(L+1)+SIG 00001780 C 00001790 F(M)=F(M)+RF*B(M) 00001800 E(M)=E(M-1)+RD*F(M) 00001810 SE=E(M)*E(M)-E(NPTS)*E(NPTS)+SE 00001820 50 SIG=E(M)*X(L)+SIG 00001830 C 00001840 C CALCULATE RELATIVE OUTPUT TRACE VARIANCE 00001850 C 00001860 VAR=SQRT((SE+E(NPTS)*E(NPTS))/SX) 00001870 C 00001880 C PLACE NORMALIZED INTERPOLATION ERROR IN ARRAY X 00001890 C 00001900 L=NMEM+1 00001910 M=NPTS-NANT 00001920 C 00001930 DO 60 I=L,M 00001940 60 X(I)=E(I+NANT)/VAR 00001950 C 00001960 GO TO 140 00001970 C 00001980 C NANT OR NMEM IS ZERO 00001990 C 00002000 70 IF (NMEM.EQ.0) GO TO 90 00002010 C 00002020 C PURELY MEMORY OPERATOR 00002030 C 00002040 VAR=SQRT((SF+F(L)*F(L))/SX) 00002050 C 00002060 DO 80 I=L,NPTS 00002070 80 X(I)=F(I)/VAR 00002080 C 00002090 GO TO 140 00002100 C 00002110 C PURELY ANTICIPATORY OPERATOR 00002120 C 00002130 90 VAR=SQRT((SB+B(N1)*B(N1))/SX) 00002140 L=NANT+1 00002150 M=NPTS-NANT 00002160 C 00002170 DO 100 I=1,M 00002180 100 X(I)=B(I+L)/VAR 00002190 C 00002200 GO TO 140 00002210 C 00002220 C TRACE SEGEMENT IS ZERO 00002230 C 00002240 110 VAR=0. 00002250 NC=NANT+NMEM 00002260 C 00002270 DO 120 I=1,NC 00002280 120 PCC(I)=0. 00002290 C 00002300 DO 130 I=1,NPTS 00002310 130 X(I)=0. 00002320 C 00002330 140 RETURN 00002340 C 00002350 END 00002360