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 C FHIGH - UPPER FREQUENCY OF FILTER AS A FRACTION OF NYQUIST C FREQUENCY 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 C IFLAG - FLAG INDICATING WHETHER FILTER COEFFICIENTS ARE TO C BE CALCULATED AGAIN. IFLAG SHOULD BE SET EQUAL C TO ZERO ON THE FIRST CALL. 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 2.1 2 APR 79 C C********************************************************************** C DIMENSION YIN(1), YOUT(1), COEF(1), TEMP(1) DATA PI/3.14159265/ C C--------------------BRANCH IF COEFFICIENTS ALREADY CALCULATED--------- C IF (IFLAG .NE. 0) GO TO 30 C C--------------------CALCULATE KAISER WEIGHTS-------------------------- C CALL KAISER (COEF, NTERMS, A) C C--------------------CALCULATE FILTER COEFFICIENTS--------------------- C DO 20 I=1,NTERMS 20 COEF(I) = COEF(I)*(SIN(PI*I*FHIGH)-SIN(PI*I*FLOW))/(PI*I) C C--------------------FILTER DATA--------------------------------------- C 30 IF (NDATA .LT. NTERMS*2+1) RETURN !RETURN IF NDATA TOO SMALL STOP = 1. !IF BANDSTOP THEN STOP=1 IF (FHIGH .GE. FLOW) STOP = 0. ! ELSE STOP=0 DO 40 I=1,NTERMS !SAVE INITIAL POINTS 40 TEMP(I) = YIN(I) DO 60 I=NTERMS+1,NDATA-NTERMS !FILTER MIDDLE POINTS SUM = 0.0 DO 50 J=1,NTERMS 50 SUM = SUM + COEF(J)*(YIN(I-J)+YIN(I+J)) 60 YOUT(I-NTERMS) = SUM + (FHIGH-FLOW+STOP)*YIN(I) DO 70 I=NDATA-NTERMS+1,NDATA !MOVE FINAL UNFILTERED PNTS 70 YOUT(I) = YIN(I) DO 80 I=NDATA-NTERMS,NTERMS+1,-1 !SHIFT MIDDLE POINTS 80 YOUT(I) = YOUT(I-NTERMS) DO 90 I=1,NTERMS !RESTORE INITIAL POINTS 90 YOUT(I) = TEMP(I) RETURN END 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 31 MARCH 79 C C********************************************************************** C DIMENSION W(1) IF (A .LE. 21.) ALPHA = 0. IF (A .GE. 50.) ALPHA = 0.1102*(A-8.7) IF (A .LT. 50. .AND. A .GT. 21.) # ALPHA = 0.5842*(A-21.)**0.4 + 0.07886*(A-21.) DENOM = BESI0(ALPHA) DO 10 K=1,N ARG = FLOAT(K)/FLOAT(N) TEMP = BESI0(ALPHA*SQRT(1.-ARG*ARG))/DENOM W(K) = TEMP 10 CONTINUE RETURN END FUNCTION BESI0 (X) C C---------------------------------------------------------------------- C C COMPUTES THE ZERO-TH ORDER MODIFIED BESSEL FUNCTION I(X) C C---------------------------------------------------------------------- C T = X/3.75 IF (T.GT.1.) GO TO 10 T = T*T BESI0=1.+T*(3.5156229+T*(3.089424+T*(1.2067492+T*(.2659732 # +T*(.0360768+T*.0045813))))) RETURN 10 T = 1./T BESI0=SQRT(1./X)*EXP(X)*(.39894228+T*(.01328592+T*.00225319 # +T*(-.00157565+T*(.00916281+T*(-.02057706 # +T*(.02635537+T*(-.01647633+T*.00392377))))))) RETURN END 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 1.0: 8 MARCH 1979 C C---------------------------------------------------------------------- C DIMENSION YIN(1), YOUT(1), W(9), Z(9), COEF(9) DATA PI /3.14159265/ DATA W /.97546667,.90469038,.79570758,.66053343,.51335543, # .36854291,.23853482,.13233204,.05468065/ C IF (NPTS .LE. 18) RETURN !RETURN IF TOO FEW PTS DO 10 I=1,9 !MOVE INITIAL UNFILTERED PTS 10 Z(I) = YIN(I) DO 20 J=1,9 !CALCULATE COEFFICIENTS 20 COEF(J) = W(J) * SIN(PI*J*FCUTOF) / (PI*J) DO 40 I=10,NPTS-9 !LOOP ON DATA SUM=0.0 DO 30 J=1,9 !COMPUTE AND ADD 9 TERMS 30 SUM = SUM + COEF(J) * (YIN(I-J)+YIN(I+J)) 40 YOUT(I-9) = SUM + FCUTOF * YIN(I) !INSERT FIRST TERM DO 50 I=NPTS-8,NPTS !MOVE FINAL UNFILTERED PTS 50 YOUT(I) = YIN(I) DO 60 I=NPTS-9,10,-1 !SHIFT MIDDLE PTS BY 9 60 YOUT(I) = YOUT(I-9) DO 70 I=1,9 !PUT IN FIRST 9 PTS 70 YOUT(I) = Z(I) RETURN END 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 27 SEP 80 C C********************************************************************** C DIMENSION Y(1), COEF(50) DATA PI/3.14159265/ C CALL KAISER (COEF, NTERMS, A) DO 40 I=1,NTERMS 40 COEF(I) = COEF(I)*(SIN(PI*I*FHIGH)-SIN(PI*I*FLOW))/(PI*I) STOP = 0. IF (FHIGH .LT. FLOW) STOP = 1. CONST = PI / FLOAT(NDATA-1) DO 50 I=1,NDATA F = CONST * (I-1) Y(I) = FHIGH - FLOW + STOP DO 50 J=1,NTERMS 50 Y(I) = Y(I) + 2.*COEF(J)*COS(F*J) RETURN END 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. IF (ACHEB .LT. 0.) ACHEB = -ACHEB ITYPE = 1 IF (ACHEB .GT. 0.) ITYPE = 2 C CALL CONST (NORDER,FCUT1,FCUT2,ITYPE,ACHEB) N = (NTYPE+1)/2 N = N*NORDER DO 10 I=1,N A(I) = 0. 10 B(I) = 0. A(N+1) = 0. A(1) = 1 B(1) = 1. N = (NORDER+1)/2 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 SUBROUTINE CONST (NORDER,FCUT1,FCUT2,ITYPE,ACHEB) COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF DATA PI/3.14159265/ C GO TO (10,20),ITYPE 10 OMEGA = 3.**(0.5/FLOAT(NORDER)) GO TO 100 20 X = SQRT(3./(10.**(0.1*ACHEB)-1.)) X = ALOG(X+SQRT(X**2-1))/FLOAT(NORDER) OMEGA = 0.5*(EXP(X)+EXP(-X)) C 100 OMEGA = ATAN(OMEGA) D1 = PI*FCUT1/2. D2 = PI*FCUT2/2. IF (FCUT1 .GT. 0.) GO TO 200 CA = SIN(OMEGA-D2)/SIN(OMEGA+D2) NTYPE = 1 RETURN C 200 IF (FCUT2 .LT. 1.) GO TO 300 CA = COS(OMEGA-D1)/COS(OMEGA+D1) NTYPE = 2 RETURN C 300 CA = COS(D1+D2)/COS(D2-D1) IF (FCUT2 .LT. FCUT1) GO TO 400 CK = SIN(OMEGA)*COS(D2-D1)/(COS(OMEGA)*SIN(D2-D1)) NTYPE = 3 RETURN C 400 CK = SIN(OMEGA)*SIN(D1-D2)/(COS(OMEGA)*COS(D1-D2)) NTYPE = 4 RETURN END 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. 2) GO TO 10 THETA = FLOAT(NORDER+2*I-1)*PI/FLOAT(2*NORDER) PCREAL = COS(THETA) PREAL = PCREAL PCIMAG = SIN(THETA) PIMAG = PCIMAG RETURN C 10 IF (I .NE. 1) GO TO 20 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)) C 20 THETA = FLOAT(2*I-1)*PI/FLOAT(2*NORDER) PREAL = PCREAL * SIN(THETA) PIMAG = PCIMAG * COS(THETA) RETURN END 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 RETURN END SUBROUTINE TRNSFM (I,NORDER) COMMON /COEFS/NTYPE,CA,CK,PREAL,PIMAG,A0,B1,B0,C(5),D(5),NCOEF GO TO (10,20,30,40),NTYPE 10 C0=-CA C1 = 1. GO TO 25 20 C0 = CA C1 = -1. 25 C2 = 0. D1 = -CA D2 = 0. NCOEF = 3 GO TO 100 30 C0 = (1.-CK)/(1.+CK) C1 = 2.*CA*CK/(1.+CK) C2 = -1. D1 = -C1 D2 = -C0 NCOEF = 5 GO TO 100 40 C0 = (1.-CK)/(1.+CK) C1 = -2.*CA/(1.+CK) C2 = 1. D1 = C1 D2 = C0 NCOEF = 5 100 IF (2*I .GT. NORDER) GO TO 200 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 GO TO 300 200 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 300 DO 400 J=1,NCOEF C(J) = C(J)/D(NCOEF) 400 D(J) = D(J)/D(NCOEF) RETURN END SUBROUTINE NXT TRM (A, B, I, NORDER) 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 DO 20 K=J,1,-1 AK = 0. BK = 0. DO 10 L=1,NCOEF L1 = K-L+1 IF (L1 .GT. JOLD) GO TO 10 IF (L1 .LT. 1) GO TO 10 AK = AK + A(L1)*C(L) BK = BK + B(L1)*D(L) 10 CONTINUE A(K) = AK IF (K .LE. NORDER*M) B(K) = BK 20 CONTINUE RETURN END 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) DATA PI/3.14159265/ ZR = COS(PI*F) ZI = SIN(PI*F) ZNR = 1. ZNI = 0. HNR = 0. HNI = 0. HDR = 0. HDI = 0. 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 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 T = SQRT(HR*HR+HI*HI) HABS = T T = ATAN2(HI,HR) * 180./PI HPHASE = T RETURN END