C FLTLIB.FOR C SUBROUTINE FILTER ( YIN, YOUT, NDATA, FLOW, FHIGH, A, COEF, # TEMP, NTERMS, IFLAG ) C C---------------------------------------------------------------------- C C DESCRIPTION: C This subroutine may be used as a lowpass, highpass, bandpass C or bandstop non-recursive filter for evenly spaced data. C The k-th filtered data point is calculated from the C (k-NTERMS)-th through (k+NTERMS)-th data points using a C Kaiser window. NTERMS data points at each end of the data C remain unfiltered due to the filter algorithm. Frequencies C are expressed in terms of the Nyquist frequency, 1/2T, where C T is the time between data samples. C C PARAMETERS: C YIN - Data array to be filtered C YOUT - Data array after filtering (may be the same as YIN) C NDATA - Number of data points in YIN C FLOW - Lower frequency of filter as a fraction of Nyquist C frequency, 1/2T C FHIGH - Upper frequency of filter as a fraction of Nyquist C frequency, 1/2T C A - Size of Gibbs phenomenon wiggles in -DB C (50 is a good choice) C COEF - Array of calculated filter coefficients C (dimension=NTERMS) C TEMP - Temporary array used by subroutine (dimension=NTERMS) C NTERMS - Number of terms in the filter formula (order of filter) C IFLAG - Flag indicating whether filter coefficients are to be C calculated again. IFLAG should be set equal to zero C each time new values of FLOW and FHIGH are used. C C USAGE: C The following conditions are necessary for various types of C filters: C NO FILTERING - FLOW = 0, FHIGH = 1 C LOW PASS - FLOW = 0, 0 < FHIGH < 1 C HIGH PASS - 0 < FLOW < 1, FHIGH = 1 C BAND PASS - 0 < FLOW < FHIGH < 1 C BAND STOP - 0 < FHIGH < FLOW < 1 C C Programmer: Robert Walraven Version 3.0 1 DEC 84 C C---------------------------------------------------------------------- C PARAMETER (PI=3.14159265) DIMENSION YIN(1), YOUT(1), COEF(1), TEMP(1) C C.......If IFLAG = 0, compute Kaiser weights and coefficients C IF (IFLAG .EQ. 0) THEN CALL KAISER (COEF, NTERMS, A) DO 10 I=1,NTERMS COEF(I) = COEF(I)*(SIN(PI*I*FHIGH)-SIN(PI*I*FLOW))/(PI*I) 10 CONTINUE ENDIF C C.......If NDATA is too small, return immediately C IF (NDATA .LT. NTERMS*2+1) RETURN !RETURN IF NDATA TOO SMALL C C.......Compute bandstop factor C IF (FHIGH .LT. FLOW) THEN STOP = 1. ELSE STOP = 0. ENDIF C C.......Save initial points C DO 20 I=1,NTERMS TEMP(I) = YIN(I) 20 CONTINUE C C.......Filter middle points C DO 40 I=NTERMS+1,NDATA-NTERMS SUM = 0.0 DO 30 J=1,NTERMS SUM = SUM + COEF(J)*(YIN(I-J)+YIN(I+J)) 30 CONTINUE YOUT(I-NTERMS) = SUM + (FHIGH-FLOW+STOP)*YIN(I) 40 CONTINUE C C.......Move final unfiltered points C DO 50 I=NDATA-NTERMS+1,NDATA YOUT(I) = YIN(I) 50 CONTINUE C C.......Shift middle points C DO 60 I=NDATA-NTERMS,NTERMS+1,-1 YOUT(I) = YOUT(I-NTERMS) 60 CONTINUE C C.......Restore initial points C DO 70 I=1,NTERMS YOUT(I) = TEMP(I) 70 CONTINUE C RETURN END C********************************************************************** SUBROUTINE KAISER (W, N, A) C C---------------------------------------------------------------------- C C Computes Kaiser weights W(N,K) for digital filters. C C PARAMETERS: C W - Calculated array of Kaiser weights C N - Value of N in W(N,K), i.e., number of terms C A - Size of Gibbs phenomenon wiggles in -DB. C C Programmer: Robert Walraven Version 3.0 1 Dec 84 C C---------------------------------------------------------------------- C DIMENSION W(1) C IF (A .LE. 21.) THEN ALPHA = 0. ELSE IF (A .GE. 50.) THEN ALPHA = 0.1102 * (A-8.7) ELSE ALPHA = 0.5842*(A-21.)**0.4 + 0.07886*(A-21.) ENDIF C DENOM = BESI0(ALPHA) C DO 10 K=1,N ARG = FLOAT(K)/FLOAT(N) W(K) = BESI0(ALPHA*SQRT(1.-ARG*ARG))/DENOM 10 CONTINUE RETURN END C********************************************************************** FUNCTION BESI0 (X) C C---------------------------------------------------------------------- C C Computes the zero-th order modified bessel function I(x) C C---------------------------------------------------------------------- C PARAMETER (A1=3.5156229) PARAMETER (A2=3.0894240) PARAMETER (A3=1.2067492) PARAMETER (A4=0.2659732) PARAMETER (A5=0.0360768) PARAMETER (A6=0.0045813) C PARAMETER (B1=0.39894228) PARAMETER (B2=0.01328592) PARAMETER (B3=0.00225319) PARAMETER (B4=-.00157565) PARAMETER (B5=0.00916281) PARAMETER (B6=-.02057706) PARAMETER (B7=0.02635537) PARAMETER (B8=-.01647633) PARAMETER (B9=0.00392377) C T = X/3.75 C IF (T .LE. 1.) THEN T = T*T BESI0=1.+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+T*A6))))) ELSE T = 1./T BESI0=SQRT(1./X)*EXP(X)*(B1+T*(B2+T*B3+T*(B4+T*(B5+T*(B6 # +T*(B7+T*(B8+T*B9))))))) ENDIF C RETURN END C********************************************************************** SUBROUTINE LOWPAS ( YIN, YOUT, NPTS, FCUTOF ) C C---------------------------------------------------------------------- C C THIS SUBROUTINE REMOVES HIGH FREQUENCY COMPONENTS FROM EVENLY C SAMPLED DATA BY PASSING THE DATA THROUGH A 9-TERM LOW-PASS DIGITAL C FILTER WITH A 50 DB KAISER WINDOW. THE CUTOFF FREQUENCY MAY BE C VARIED WITH THE VARIABLE 'FCUTOF'. NINE DATA POINTS AT EACH C END OF THE DATA REMAIN UNFILTERED DUE TO THE FILTER ALGORITHM. C C PARAMETERS: C YIN - DATA ARRAY TO BE FILTERED C YOUT - DATA ARRAY AFTER FILTERING (MAY BE THE SAME AS YIN) C NPTS - NUMBER OF POINTS IN DATA ARRAY C FCUTOF - CUTOFF FREQUENCY OF FILTER AS A FRACTION OF C THE NYQUIST FREQUENCY, 1/2T, WHERE T IS THE C TIME BETWEEN DATA SAMPLES. C C PROGRAMMER: ROBERT WALRAVEN, APPLIED SCIENCE, UCD C VERSION 3.0: 1 DEC 84 C C---------------------------------------------------------------------- C DIMENSION YIN(1), YOUT(1), W(9), Z(9), COEF(9) PARAMETER (PI=3.14159265) DATA W /.97546667,.90469038,.79570758,.66053343,.51335543, # .36854291,.23853482,.13233204,.05468065/ C C.......Return if too few points C IF (NPTS .LE. 18) RETURN C C.......Move initial unfiltered points C DO 10 I=1,9 Z(I) = YIN(I) 10 CONTINUE C C.......Calculate coefficients C DO 20 J=1,9 COEF(J) = W(J) * SIN(PI*J*FCUTOF) / (PI*J) 20 CONTINUE C C.......Filter middle data C DO 40 I=10,NPTS-9 SUM=0.0 DO 30 J=1,9 SUM = SUM + COEF(J) * (YIN(I-J)+YIN(I+J)) 30 CONTINUE YOUT(I-9) = SUM + FCUTOF * YIN(I) 40 CONTINUE C C.......Move final unfiltered points C DO 50 I=NPTS-8,NPTS YOUT(I) = YIN(I) 50 CONTINUE C C.......Shift middle points by 9 C DO 60 I=NPTS-9,10,-1 YOUT(I) = YOUT(I-9) 60 CONTINUE C C.......Put in first 9 points C DO 70 I=1,9 YOUT(I) = Z(I) 70 CONTINUE C RETURN END C********************************************************************** SUBROUTINE H OF Z NR (Y, NDATA, NTERMS, A, FLOW, FHIGH) C C---------------------------------------------------------------------- C C COMPUTES TRANSFER FUNCTION OF NON-RECURSIVE FILTER. C C FILTER DESIGN PARAMETERS: C Y - ARRAY WHERE TRANSFER FUNCTION IS STORED C NDATA - NUMBER OF ELEMENTS IN Y C NTERMS - NUMBER OF FILTER TERMS C A - SIZE OF GIBBS PHENOMENON RIPPLES IN -DB C FLOW - LOWER CUTOFF OF FILTER AS MULTIPLE OF NYQUIST FREQ. C FHIGH - UPPER CUTOFF OF FILTER AS MULTIPLE OF NYQUIST FREQ. C C PROGRAMMER: ROBERT WALRAVEN VERSION 3.0 1 DEC 84 C C********************************************************************** C DIMENSION Y(1), COEF(50) PARAMETER (PI=3.14159265) C C.......Compute Kaiser weights C CALL KAISER (COEF, NTERMS, A) C C.......Compute coefficients C DO 10 I=1,NTERMS COEF(I) = COEF(I)*(SIN(PI*I*FHIGH)-SIN(PI*I*FLOW))/(PI*I) 10 CONTINUE C C.......Compute bandstop factor C IF (FHIGH .LT. HFLOW) THEN STOP = 1. ELSE STOP = 0. ENDIF C C.......Compute transfer function C CONST = PI / FLOAT(NDATA-1) C DO 30 I=1,NDATA F = CONST * (I-1) Y(I) = FHIGH - FLOW + STOP DO 20 J=1,NTERMS Y(I) = Y(I) + 2.*COEF(J)*COS(F*J) 20 CONTINUE 30 CONTINUE C RETURN END C************************************************************************* SUBROUTINE FLT REC (N,FLOW,FHIGH,IKIND,X,Y,Z,NPTS) C C------------------------------------------------------------------------- C C Passes data through an N order recursive digital filter. C C N - Desired order of filter C FLOW - Lower frequency of filter (range = 0 to 1) C FHIGH - Upper frequency of filter (range = 0 to 1) C IKIND - 1 for Butterworth, 2 for 0.1 DB Chebychev C X - Input data array C Y - Output data array C Z - Dummy array (dimensioned the same as X and Y) C NPTS - Number of data points in X C C Programmer: Robert Walraven Version 3.0 1 Dec 84 C C------------------------------------------------------------------------- C DIMENSION X(1), Y(1), Z(1), A(3), B(2) C IF (FLOW .EQ. 0.) THEN F = FHIGH ITYPE = 1 ELSE IF (FHIGH .EQ. 0.) THEN F = FLOW ITYPE = 2 ELSE F = FLOW IF (FHIGH .GE. FLOW) THEN ITYPE = 3 ELSE ITYPE = 4 ENDIF ENDIF C DO 10 I=1,NPTS Z(I) = X(I) 10 CONTINUE C CALL R FLT LW (N,F,IKIND,Y,Z,NPTS) C IF (ITYPE .EQ. 1) RETURN C IF (ITYPE .EQ. 2) THEN DO 20 I=1,NPTS Y(I) = X(I) - Y(I) 20 CONTINUE RETURN ENDIF C DO 30 I=1,NPTS Z(I) = X(I) - Y(I) 30 CONTINUE CALL R FLT LW (N,FHIGH,IKIND,Y,Z,NPTS) C IF (ITYPE .EQ. 4) THEN DO 40 I=1,NPTS Y(I) = X(I) - Y(I) 40 CONTINUE ENDIF C RETURN END C************************************************************************* SUBROUTINE R FLT LW (N,F,IKIND,Y,Z,NPTS) C C------------------------------------------------------------------------- C C Passes data through a second order low pass recursive filter C that will eventually be applied N/2 times to produce N order C filtering. C C N - Order of ultimate filter C F - Cutoff frequency desired (range = 0 to 1) C IKIND - 1 for Butterworth, 2 for 0.1 DB Chebychev C Y - Output data array C Z - Input data array (= output on return) C NPTS - Number of data points in Y,Z C C Programmer: Robert Walraven Version 3.0 1 Dec 84 C C------------------------------------------------------------------------ C DIMENSION Y(1), Z(1), CONST(10), COEF(8) PARAMETER (HALFPI=1.57079633) C DATA CONST /1., 1.3160740, 1.5033033, 1.6404924, 1.7504911, 1 1., 1.2787175, 1.4320239, 1.5385184, 1.6202971/ C DATA COEF /0.92103004, 1.3722269, 0.29289321, 0.6976604, 1 0.00000000, 1.2900624, 0.17157287, 0.5005791/ C NTIMES = (N+1) / 2 IF (NTIMES .LT. 1) NTIMES = 1 IF (NTIMES .GT. 5) NTIMES = 5 IF (NKIND .LT. 1) NKIND = 1 IF (NKIND .GT. 2) NKIND = 2 C C = CONST( NTIMES + 5*(IKIND-1) ) OMEGA = COEF(IKIND) A0 = COEF (IKIND+2) B1 = COEF (IKIND+4) B0 = COEF (IKIND+6) C D = ATAN (C*SIN(HALFPI*F)/COS(HALFPI*F)) IF (D .GT. HALFPI) D = HALFPI C CA = SIN(OMEGA-D) / SIN(OMEGA+D) CA2 = CA*CA C1 = 1. + B0*CA2 - B1*CA C2 = A0*(1.-CA)**2/C1 C3 = (B1*(1.+CA2)-2.*CA*(1.+B0))/C1 C4 = (CA2+B0-B1*CA)/C1 C DO 30 J=1,NTIMES Y(1) = Z(1) Y(2) = Z(2) DO 10 I=3,NPTS Y(I) = C2*(Z(I)+2.*Z(I-1)+Z(I-2))-C3*Y(I-1)-C4*Y(I-2) 10 CONTINUE DO 20 I=1,NPTS Z(I) = Y(I) 20 CONTINUE 30 CONTINUE C RETURN END C************************************************************************* SUBROUTINE COEFS (A, B, NORDER, FCUT1, FCUT2, ACHEB) C C COMPUTES RECURSIVE DIGITAL FILTER COEFFICIENTS. C C A,B ARE ARRAYS OF COEFFICIENTS TO BE RETURNED. A AND B C SHOULD BE DIMENSIONED GREATER THAN OR EQUAL TO C NORDER+1 AND NORDER, RESPECTIVELY, IN THE MAIN PROGRAM. C NORDER IS THE ORDER OF THE FILTER. C FCUT1,FCUT2 ARE LOWER AND UPPER CUTOFFS OF FILTER, I.E., THE C FILTER PASSES FREQUENCIES BETWEEN FCUT1 AND FCUT2. C FOR DESIRED TYPE OF FILTER CHOOSE FCUT1 AND FCUT2 C AS FOLLOWS: C FCUT1 FCUT2 C LOWPASS 0 0 TO 1 C HIGHPASS 0 TO 1 1 C BANDPASS 0 TO 1 >FCUT1 C BANDSTOP 0 TO 1 0, ACHEB IS ATTENUATION FACTOR OF CHEBYCHEV FILTER C IN DB WHEN ITYPE = 2 C DIMENSION A(1),B(1) COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF C IF (NORDER .LT. 1) NORDER=1 IF (FCUT1 .EQ. FCUT2) FCUT2 = FCUT1 + .2 IF (FCUT1.LT.0. .OR. FCUT1.GT.1.) FCUT1=0. IF (FCUT2.LT.0. .OR. FCUT2.GT.1.) FCUT2=1. C ACHEB = ABS (ACHEB) IF (ACHEB .EQ. 0.) THEN ITYPE = 1 ELSE ITYPE = 2 ENDIF C CALL CONST (NORDER,FCUT1,FCUT2,ITYPE,ACHEB) N = (NTYPE+1)/2 N = N*NORDER DO 10 I=1,N A(I) = 0. B(I) = 0. 10 CONTINUE A(N+1) = 0. A(1) = 1 B(1) = 1. N = (NORDER+1)/2 C DO 100 I=1,N II = I CALL POLE(II,NORDER,ITYPE,ACHEB) IF (2*I .LE. NORDER) CALL S TO Z CALL TRNSFM (II,NORDER) CALL NXT TRM (A,B,II,NORDER) 100 CONTINUE RETURN END C*********************************************************************** SUBROUTINE CONST (NORDER,FCUT1,FCUT2,ITYPE,ACHEB) C C Computes constants for recursive digital filter C COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF DATA PI/3.14159265/ C IF (ITYPE .EQ. 1) THEN OMEGA = 3.**(0.5/FLOAT(NORDER)) ELSE X = SQRT(3./(10.**(0.1*ACHEB)-1.)) X = ALOG(X+SQRT(X**2-1))/FLOAT(NORDER) OMEGA = 0.5*(EXP(X)+EXP(-X)) ENDIF C OMEGA = ATAN(OMEGA) D1 = PI*FCUT1/2. D2 = PI*FCUT2/2. IF (FCUT1 .LE. 0.) THEN CA = SIN(OMEGA-D2)/SIN(OMEGA+D2) NTYPE = 1 RETURN ENDIF C IF (FCUT2 .GE. 1.) THEN CA = COS(OMEGA-D1)/COS(OMEGA+D1) NTYPE = 2 RETURN ENDIF C CA = COS(D1+D2)/COS(D2-D1) IF (FCUT2 .GE. FCUT1) THEN CK = SIN(OMEGA)*COS(D2-D1)/(COS(OMEGA)*SIN(D2-D1)) NTYPE = 3 RETURN ENDIF C CK = SIN(OMEGA)*SIN(D1-D2)/(COS(OMEGA)*COS(D1-D2)) NTYPE = 4 RETURN END C************************************************************************* SUBROUTINE POLE (I, NORDER, ITYPE, ACHEB) C C CALCULATES POLES OF RECURSIVE FILTER OF DESIRED TYPE C COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF DATA PI,PCREAL,PCIMAG/3.14159265,0.,0./ C IF (ITYPE .EQ. 1) THEN THETA = FLOAT(NORDER+2*I-1)*PI/FLOAT(2*NORDER) PCREAL = COS(THETA) PREAL = PCREAL PCIMAG = SIN(THETA) PIMAG = PCIMAG RETURN ENDIF C IF (I .EQ. 1) THEN X = 1./SQRT(10.**(0.1*ACHEB)-1.) X = ALOG(X+SQRT(X**2+1.))/FLOAT(NORDER) PCREAL = -0.5*(EXP(X)-EXP(-X)) PCIMAG = 0.5*(EXP(X)+EXP(-X)) ENDIF C THETA = FLOAT(2*I-1)*PI/FLOAT(2*NORDER) PREAL = PCREAL * SIN(THETA) PIMAG = PCIMAG * COS(THETA) C RETURN END C*********************************************************************** SUBROUTINE S TO Z C C PERFORMS VARIABLE TRANSFORMATION FROM S TO Z C COMMON/COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF C P2 = PREAL**2 + PIMAG**2 DENOM = 1 + P2 -2.*PREAL QREAL = (1.-P2)/DENOM QIMAG = 2.*PIMAG/DENOM A0 = P2/DENOM B1 = -2.*QREAL B0 = QREAL**2 + QIMAG**2 C RETURN END C*********************************************************************** SUBROUTINE TRNSFM (I,NORDER) C C Computes recursive digital filter transformations C COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF C IF (NTYPE .EQ. 1) THEN C0=-CA C1 = 1. C2 = 0. D1 = -CA D2 = 0. NCOEF = 3 ELSE IF (NTYPE .EQ. 2) THEN C0 = CA C1 = -1. C2 = 0. D1 = -CA D2 = 0. NCOEF = 3 ELSE IF (NTYPE .EQ. 3) THEN C0 = (1.-CK)/(1.+CK) C1 = 2.*CA*CK/(1.+CK) C2 = -1. D1 = -C1 D2 = -C0 NCOEF = 5 ELSE IF (NTYPE .EQ. 4) THEN C0 = (1.-CK)/(1.+CK) C1 = -2.*CA/(1.+CK) C2 = 1. D1 = C1 D2 = C0 NCOEF = 5 ENDIF C IF (2*I .LE. NORDER) THEN C(5) = A0*(C2+D2)**2 C(4) = 2.*A0*(C1+D1)*(C2+D2) C(3) = A0*((C1+D1)**2+2.*(1.+C0)*(C2+D2)) C(2) = 2.*A0*(1.+C0)*(C1+D1) C(1) = A0*(1.+C0)**2 D(5) = C2**2+B1*C2*D2+B0*D2**2 D(4) = 2.*C1*C2+B1*(C1*D2+C2*D1)+2.*B0*D1*D2 D(3) = 2.*C0*C2+C1**2+B1*(C2+C1*D1+C0*D2)+B0*(2.*D2+D1**2) D(2) = 2.*C0*C1+B1*(C1+C0*D1)+2.*B0*D1 D(1) = C0**2+B1*C0+B0 ELSE NCOEF = (NCOEF+1)/2 A0= PREAL/(PREAL-1.) B0 = (PREAL+1.)/(PREAL-1.) C(3) = A0*(C2+D2) C(2) = A0*(C1+D1) C(1) = A0*(C0+1.) D(3) = C2+B0*D2 D(2) = C1+B0*D1 D(1) = C0+B0 ENDIF C DO 100 J=1,NCOEF C(J) = C(J)/D(NCOEF) D(J) = D(J)/D(NCOEF) 100 CONTINUE C RETURN END C************************************************************************* SUBROUTINE NXT TRM (A, B, I, NORDER) C C Computes next term of recursive digital filter C COMMON/COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF DIMENSION A(1),B(1) C N = NCOEF M = (NTYPE+1)/2 IF (2*I .GT. NORDER) N=2*NCOEF-1 J = (N-1)*I+1 JOLD = J-N+1 C DO 20 K=J,1,-1 AK = 0. BK = 0. DO 10 L=1,NCOEF L1 = K-L+1 IF (L1.LE.JOLD .AND. L1.GE.1) THEN AK = AK + A(L1)*C(L) BK = BK + B(L1)*D(L) ENDIF 10 CONTINUE A(K) = AK IF (K .LE. NORDER*M) B(K) = BK 20 CONTINUE C RETURN END C************************************************************************ SUBROUTINE H OF Z R (A, B, NORDER, F, HABS, HPHASE) C C COMPUTES TRANSFER FUNCTION FOR RECURSIVE DIGITAL FILTER. C C A,B ARE ARRAYS OF FILTER COEFFICIENTS. C NORDER IS ORDER OF FILTER. C F IS FREQUENCY AT WHICH TRANSFER FUNCTION IS TO BE COMPUTED. C (F IS IN UNITS OF 1/2T, SO RANGE IS 0. TO 1.) C HABS,HPHASE ARE ABSOLUTE VALUE AND PHASE OF TRANSFER FUNCTION. C DIMENSION A(1),B(1) PARAMETER (PI=3.14159265) C ZR = COS(PI*F) ZI = SIN(PI*F) ZNR = 1. ZNI = 0. HNR = 0. HNI = 0. HDR = 0. HDI = 0. C DO 10 I=1,NORDER HNR = HNR + A(I)*ZNR HNI = HNI + A(I)*ZNI HDR = HDR + B(I)*ZNR HDI = HDI + B(I)*ZNI ZTEMP = ZNR ZNR = ZNR*ZR - ZNI*ZI ZNI = ZTEMP*ZI + ZNI*ZR 10 CONTINUE C HNR = HNR + A(NORDER+1)*ZNR HNI = HNI + A(NORDER+1)*ZNI HDR = HDR + ZNR HDI = HDI + ZNI T = HDR*HDR + HDI*HDI HR = (HNR*HDR + HNI*HDI)/T HI = (HNI*HDR - HNR*HDI)/T HABS = SQRT(HR*HR+HI*HI) HPHASE = ATAN2(HI,HR) * 180./PI C RETURN END