CTITLEMCOFGN -- COMPUTES 99 WIENER INTERPOLATION FILTERS 00010001 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR C. J. SICKING 00020001 CA DESIGNER C. J. SICKING 00030001 CA LANGUAGE FORTRAN 00040001 CA SYSTEM IBM AND CRAY 00050001 CA WRITTEN SEPTEMBER, 1980 00060001 C REVISED 03-05-82 ESN. FOR THE CRAY. 00070001 CA 00080001 CA 00090001 CA CALL MCOFGN (LFLTR, MXFREQ, ACOR, CCOR, FLTR, WORK, COEF) 00100000 CA 00110001 CA IN/OUT ARGUMENT TYPE DESCRIPTION 00120001 CA 00130001 CA IN LFLTR I4 NUMBER OF COEFICIENTS IN THE FILTERS 00140001 CA IN MXFREQ R4 MAXIMUM FREQUENCY AS A FRACTION OF NYQUIST00150001 CA IN ACOR R4 ARRAY FOR AUTOCORRELATION 00160001 CA IN CCOR R4 ARRAY FOR CROSS CORRELATION 00170001 CA IN FLTR R4 FILTER ARRAY USED BY MWIENR 00180001 CA IN WORK R4 WORK ARRAY USED BY MWIENR 00190001 CA OUT COEF R4 COEFICIENT ARRAY. 99 FILTERS OF LENGTH 00200000 CA LFLTR. 00210001 CA 00220001 CA 00230001 CA MCOFGN COMPUTES 99 INTERPOLATION FILTERS ON EACH CALL. EACH 00240001 CA FILTER IS LFLTR COEFICIENTS IN LENGTH. THESE 99 FILTERS CAN BE 00250001 CA USED TO COMPUTE THE VALUE OF A SIGNAL TO THE NEAREST 0.01 SAMPLE 00260001 CA INTERVAL. COMPUTING THE VALUE TO THE NEAREST 0.01 SAMPLE 00270001 CA INTERVAL IS ADEQUATE TO MAINTAIN 30 TO 40 DB FIDELITY IN THE 00280001 CA INTERPOLATION PROCESS. 00290001 CA 00300000 C 00310001 C 00320001 C IF WE HAVE A SIGNAL WITH SAMPLES AT LOCATIONS X1,X2,X3 AND X4 AND 00330001 C WE WANT THE VALUE OF THE SIGNAL AT LOCATION P, WE WOULD FIRST 00340001 C COMPUTE THE FRACTIONAL DISTANCE BY DIVIDING THE DISTANCE P-X2 BY 00350001 C THE DISTANCE X3-X2 AND THEN ROUND THE RESULT TO THE NEAREST 0.01 00360001 C AND MULTIPLY BY 100. THE RESULT IS A NUMBER BETWEEN 0 AND 100. 00370001 C THIS IS THE INDEX TO THE FILTER ARRAY AND IS USED TO DETERMINE 00380001 C WHICH OF THE 99 FILTERS TO USE IN THE INTERPOLATION. IF THE 00390001 C RESULT IS 0 OR 100, P IS WITHIN PLUS OR MINUS 0.005 * (SAMPLE 00400000 C INTERVAL) OF X2 OR X3 AND IS CLOSE ENOUGH TO USE THE 00410001 C KNOWN VALUE WITHOUT GOING THROUGH THE INTERPOLATION STEP. 00420001 C 00430001 C 00440001 C THE FILTERS COMPUTED ARE WIENER FILTERS. A CORRELATION FUNCTION 00450001 C IS REQUIRED IN ORDER TO COMPUTE THE FILTERS. THE FILTERS ARE 00460001 C COMPUTED USING THE MATRIX EQUATION: 00470001 C 00480001 C A F = C 00490001 C 00500000 C WHERE A IS THE AUTOCORRELATION MATRIX AND C IS THE CROSS COR- 00510001 C RELATION MATRIX. THE EQUATION IS SOLVED USING A MODIFIED VERSION 00520001 C OF SUBROUTINE EUREKA (RENAMED MWIENR) FROM ROBINSON'S BOOK. 00530001 C THE CORRELATION FUNCTION USED IS A BANDLIMITED SINC FUNCTION. 00540001 C FOR DETAILS OF THE INTERPOLATION PROCEDURE AND HOW THE 00550001 C PARAMETERS ARE RELATED SEE SICKING' PH.D. DISSERTATION. LFLTR 00560001 C AND MXFREQ ARE INTERRELATED AND THEIR VALUES SHOULD BE CHOSEN 00570001 C BASED ON KNOWGLEDGE OF THE MAXIMUM FREQUENCY PRESENT IN THE SIGNAL00580001 C AND THE CURVES PRESENTED IN SICKING'S DISSERTATION CHAPTER 4. 00590001 C 00600000 C 00610001 SUBROUTINE MCOFGN (LFLTR, MXFREQ, ACOR, CCOR, FLTR, WORK, COEF) 00620001 C 00630001 REAL ACOR ( 1) 00640001 REAL CCOR ( 1) 00650001 REAL COEF (LFLTR, 99) 00660001 REAL FLTR ( 1) 00670001 REAL WORK ( 1) 00680001 C 00690001 REAL MXFREQ 00700000 C 00710001 C COMPUTE AUTOCORRELATION VALUES 00720001 C 00730001 C THE AUTOCORRELATION MATRIX IS THE SAME FOR ALL FILTERS COMPUTED 00740001 C SO IT NEEDS TO BE COMPUTED ONLY ONCE. 00750001 C 00760001 C THE PARAMETER CORD IS THE LAG IN NUMBER OF SAMPLE INTERVALS. 00770001 C 00780001 PI = 3.141592654 00790001 CORD = 1.0 00800000 ACOR(1) = 1.0 00810001 C 00820001 DO 10 I = 2, LFLTR 00830001 TX = CORD * PI * MXFREQ 00840001 ACOR(I) = SIN(TX) / TX 00850001 C 00860001 10 CORD = CORD + 1.0 00870001 C 00880001 C COMPUTE EACH OF FIRST 50 FILTERS 00890001 C 00900000 C LCN IS INCREMENTED FROM 1 TO 50 WHILE FRAC IS INCREMENTED FROM 00910001 C O.01 TO 0.50. THE CROSS CORRELATION ARRAY IS FILLED ,THE FILTER 00920001 C IS COMPUTED, AND THE FILTER IS STORED IN THE OUTPUT COEF ARRAY. 00930001 C 00940001 C THE PARAMETER CORDIS IS THE CORRELATION LAG IN SAMPLE SPACINGS 00950001 C BETWEEN A GIVEN SAMPLE AND THE DESIRED SAMPLE LOCATION. 00960001 C 00970001 LCN = 1 00980001 FRAC = 0.01 00990001 C 01000000 20 CORDIS = -((LFLTR / 2) - 1) - FRAC 01010001 C 01020001 C COMPUTE CROSS CORRELATION VALUES 01030001 C 01040001 DO 50 I = 1, LFLTR 01050001 TX = CORDIS * PI * MXFREQ 01060001 IF (TX .EQ. 0.0) GO TO 30 01070001 VAL = SIN(TX) / TX 01080001 GO TO 40 01090001 C 01100000 30 VAL = 1.0 01110001 C 01120001 40 CCOR(I) = VAL 01130001 C 01140001 50 CORDIS = CORDIS + 1.0 01150001 C 01160001 C COMPUTE WIENER FILTER AND STORE IN COEF ARRAY 01170001 C 01180001 CALL MWIENR (LFLTR, ACOR, CCOR, WORK, FLTR) 01190001 C 01200000 DO 60 I = 1, LFLTR 01210001 60 COEF(I, LCN) = FLTR(I) 01220001 C 01230001 C INCREMENT FRAC AND LCN AND RETURN TO COMPUTE THE NEXT FILTER. 01240001 C 01250001 FRAC = FRAC + 0.01 01260001 LCN = LCN + 1 01270001 IF (LCN .LT. 51) GO TO 20 01280001 C 01290001 C THE FIRST 50 FILTERS HAVE BEEN COMPUTED. FILTER 50 IS THE 01300000 C FILTER TO BE USED FOR MIDPOINT INTERPOLATION. 01310001 C THERE IS SYMMETRY ABOUT MIDPOINT. REVERSE FIRST 49 FILTERS AND 01320001 C STORE IN LAST 49. 01330001 C 01340001 DO 80 I = 1, 49 01350001 J = 50 - I 01360001 K = 50 + I 01370001 C 01380001 DO 70 L = 1, LFLTR 01390001 M = LFLTR - L + 1 01400000 70 COEF(M, K) = COEF(L, J) 01410001 C 01420001 80 CONTINUE 01430001 C 01440001 RETURN 01450001 END 01460001