CTITLEMAHLBT - TIME DOMAIN HILBERT TRANSFORMER 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 SYSTEMS IBM/CRAY CA LANGUAGE VS FORTRAN VERSION 2.2 CA WRITTEN 12-03-87 C REVISED 05-23-91 JJC - CHANGED INTEGER NONE TO (A-Z). CA CA CALLING SEQUENCE: CA CALL MAHLBT(A,H,NT,IS1,IS2,FTRLEN) CA CA IN/OUT ARGUMENT TYPE DESCRIPTION CA ------ -------- ---- ----------- CA CA IN A (R4) THE INPUT VECTOR TO HILBERT TRANSFORM CA OUT H (R4) THE TRANSFORMED VECTOR CA IN NT I4 THE DIMENSIONED LENGTHS OF 'A' AND 'H'. CA IN IS1 I4 THE INDEX OF THE FIRST SAMPLE TO CA TRANSFORM CA IN IS2 I4 THE INDEX OF THE LAST SAMPLE TO CA TRANSFORM CA IN FTRLEN I4 THE LENGTH OF THE HILBERT TRANSFORMER CA (VALID LENGTHS ARE: 0, 31, 43, 63, 95) CA (A LENGTH OF 0 SUPPRESSES THE HILBERT CA TRANSFORMATION) CA CA PURPOSE: CA CA THIS SUBROUTINE USES THE REAL PART OF COMPLEX ARRAY "A" TO CA COMPUTE ITS HILBERT TRANSFORM, WHICH IT PLACES IN THE CA IMAGINARY PART OF THE ARRAY. CA CA CA SUBROUTINES CALLED: NONE CA CA CA METHOD: CA CA THE DATA IS CONVOLVED WITH A FIR HILBERT TRANSFORMER SHOWN IN CA RABINER AND SCHAFER, "ON THE BEHAVIOR OF MINIMAX FIR DIGITAL HILBERT CA TRANSFORMERS", BELL SYST. TECH. J. VOL. 53 NO. 2, P 380, (1974). C C THIS SUBROUTINE MUST BE VECTORIZED IN ORDER TO RUN EFFICIENTLY. C SUBROUTINE MAHLBT(A,H,NT,IS1,IS2,FTRLEN) C CJJ IMPLICIT NONE IMPLICIT INTEGER (A-Z) C INTEGER NLEN, MAXLEN, LEN1, LEN2, LEN3, LEN4 INTEGER L1, L2, L3, L4, MAXL, NT, I, NDX, FTSAVE INTEGER FTRLEN, L, LT1, LT2, M, M2, LIM1, LIM2, N REAL A, H, HLTR, FILTER, H1, H2, H3, H4 PARAMETER (NLEN=4, MAXLEN=95) PARAMETER (LEN1=95,LEN2=63,LEN3=43,LEN4=31) PARAMETER (L1=((LEN1+1)/4)) PARAMETER (L2=((LEN2+1)/4)) PARAMETER (L3=((LEN3+1)/4)) PARAMETER (L4=((LEN4+1)/4)) PARAMETER (MAXL=((MAXLEN+1)/4)) DIMENSION A(NT), H(NT) DIMENSION HLTR(-MAXLEN:MAXLEN) DIMENSION FILTER(MAXL,NLEN) DIMENSION H1(L1), H2(L2), H3(L3), H4(L4) EQUIVALENCE (FILTER(1,1),H1) EQUIVALENCE (FILTER(1,2),H2) EQUIVALENCE (FILTER(1,3),H3) EQUIVALENCE (FILTER(1,4),H4) INTEGER IS1, IS2, VALID(NLEN) SAVE FTSAVE, HLTR EXTERNAL ARSET DATA VALID/LEN1, LEN2, LEN3, LEN4/ C C COEFFICIENTS FOR 95 POINT FILTER, FL=0.01 C DATA H1/-0.0130099, -0.0045718, -0.0053689, -0.0062800, * -0.0072616, -0.0083873, -0.0096455, -0.0110350, * -0.0126022, -0.0143770, -0.0163895, -0.0186922, * -0.0213465, -0.0244424, -0.0281175, -0.0325665, * -0.0380864, -0.0451608, -0.0546233, -0.0680547, * -0.0888468, -0.1258168, -0.2112989, -0.6363167/ C C COEFFICIENTS FOR 63 POINT FILTER, FL=0.02 C DATA H2/-0.0055706, -0.0044618, -0.0062078, -0.0083775, * -0.0110475, -0.0143195, -0.0183266, -0.0232716, * -0.0294397, -0.0373177, -0.0477262, -0.0622273, * -0.0841979, -0.1224340, -0.2092445, -0.6356280/ C C COEFFICIENTS FOR 43 POINT FILTER, FL=0.02 C DATA H3/-0.0216528, -0.0146215, -0.0196329, -0.0259590, * -0.0340917, -0.0448233, -0.0597557, -0.0821807, * -0.1209598, -0.2083547, -0.6353344/ C C COEFFICIENTS FOR 31 POINT FILTER, FL=0.05 C DATA H4/-0.0041956, -0.0092821, -0.0188358, -0.0344010, * -0.0595516, -0.1030376, -0.1968315, -0.6313536/ C C INITIALIZE THE OUTPUT ARRAY. C C CALL ARSET(H,NT,0.0) DO 10 I=1,NT 10 H(I) = 0.0 C C IS THIS FTRLEN A VALID LENGTH? C IF (FTSAVE .NE. FTRLEN) THEN DO 20 I=1,NLEN NDX = I FTSAVE = VALID(NDX) 20 IF(FTSAVE .EQ. FTRLEN) GO TO 30 FTSAVE = 0 RETURN C C INITIALIZE THE FILTER FUNCTION. C 30 L = (FTSAVE + 1)/4 DO 40 I=0,L-1 HLTR(I) = -FILTER(L-I,NDX) HLTR(-I-1) = -HLTR(I) 40 CONTINUE ENDIF C C PERFORM THE CONVOLUTION C L = (FTRLEN + 1)/4 LT1 = MAX0(1,IS1) LT2 = MIN0(NT,IS2) IF(LT1 .GT. LT2) RETURN DO 100 M=-L,L-1 M2 = M*2 C LIM1 = MAX0(1,IS1,M2+2) C LIM2 = MIN0(NT,IS2,NT+M2+1) LIM1 = M2 + 2 LIM2 = NT + M2 + 1 IF(LIM1 .LT. LT1) LIM1 = LT1 IF(LIM2 .GT. LT2) LIM2 = LT2 DO 150 N=LIM1, LIM2 H(N) = H(N) + HLTR(M)*A(N-M2-1) 150 CONTINUE 100 CONTINUE RETURN C END C================================================= C C HERE IS A SIMPLE PROGRAM TO TEST MAHLBT. C C C PARAMETER (IPUNIT=2) C INTEGER FTRLEN C PARAMETER (FTRLEN=43) C PARAMETER (NS=1024) C DIMENSION A(NS), H(NS) C C DO 100 I=1,NS C J = I - NS/2 C A(I) = COS(J/10.)*EXP(-(J/40.)**2) C100 CONTINUE C C C CALL MAHLBT(A,H,NS,1,NS,FTRLEN) C C CALL USTLGF(IPUNIT,'REAL',A,NS,1.0) C CALL USTLGF(IPUNIT,'IMAG',H,NS,1.0) C CALL USTLND C STOP END