C @(#)WFIR.F 1.2 1/11/91 C MAIN PROGRAM: WINDOW DESIGN OF LINEAR PHASE, LOWPASS, HIGHPASS C BANDPASS, AND BANDSTOP FIR DIGITAL FILTERS C C AUTHORS: LAWRENCE R. RABINER AND CAROL A. MCGONEGAL C BELL LABORATORIES, MURRAY HILL, NEW JERSEY, 07974 C C MODIFIED JAN. 1978 BY DOUG PAUL, MIT LINCOLN LABORATORIES C TO INCLUDE SUBROUTINES FOR OBTAINING FILTER BAND EDGES AND RIPPLES C C MODIFIED JULY. 1990 BY EDWARD LEE, UCB, FOR GABRIEL COMPATIBILITY C C INPUT: NF IS THE FILTER LENGTH IN SAMPLES C 3 <= NF <= 1024 C C ITYPE IS THE WINDOW TYPE C ITYPE = 1 RECTANGULAR WINDOW C ITYPE = 2 TRIANGULAR WINDOW C ITYPE = 3 HAMMING WINDOW C ITYPE = 4 GENERALIZED HAMMING WINDOW C ITYPE = 5 HANNING WINDOW C ITYPE = 6 KAISER (I0-SINH) WINDOW C ITYPE = 7 CHEBYSHEV WINDOW C C JTYPE IS THE FILTER TYPE C JTYPE = 1 LOWPASS FILTER C JTYPE = 2 HIGHPASS FILTER C JTYPE = 3 BANDPASS FILTER C JTYPE = 4 BANDSTOP FILTER C C FC IS THE NORMALIZED CUTOFF FREQUENCY C 0 <= FC <= 0.5 C C FL AND FH ARE THE NORMALIZED FILTER CUTOFF FREQUENCIES C 0 <= FL <= FH <= 0.5 C C IWP OPTIONALLY PRINTS OUT THE WINDOW VALUES C IWP = 0 DO NOT PRINT C IWP = 1 PRINT C C IMD REQUESTS ADDITIONAL RUNS C IMD = 1 NEW RUN C IMD = 0 TERMINATES PROGRAM C C----------------------------------------------------------------------- C DIMENSION W(512), G(512) INTEGER OTCD1, OTCD2 CHARACTER FNAME*100 CHARACTER ANSWER*1 C PI = 4.0*ATAN(1.0) TWOPI = 2.0*PI C C DEFINE I/O DEVICE CODES C INPUT: INPUT TO THIS PROGRAM IS USER-INTERACTIVE C THAT IS - A QUESTION IS WRITTEN ON THE USER C TERMINAL (OTCD1) AND THE USER TYPES IN THE ANSWER. C C OUTPUT: STANDARD OUTPUT IS OTCD1 AND OTCD2 C C OTCD1 = I1MACH(4) OTCD1 = 6 C OTCD2 = I1MACH(2) OTCD2 = 6 C WRITE (*,'(,A)') ' PROVISIONAL WINDOW FIR FILTER DESIGN ' WRITE (*,'(,A)') ' USE AT YOUR OWN RISK --------------- ' C C INPUT THE FILTER LENGTH(NF), WINDOW TYPE(ITYPE) AND FILTER TYPE(JTYPE) C 10 WRITE (*,'(/,A)') ' ENTER NAME OF INPUT COMMAND FILE (PRESS FOR MANUAL ENTRY, ' WRITE (*,'(,A\)') ' SORRY, NO TILDE-EXPANSION. GIVE PATH RELATIVE + TO YOUR HOME OR STARTUP DIRECTORY): ' READ (*,'(A12)') FNAME IF (FNAME .EQ. ' ') THEN INCOD = 5 ELSE INCOD = 3 OPEN (3, FILE=FNAME) ENDIF C 15 WRITE (*,'(/,A)') ' ENTER FILTER LENGTH: ' READ (INCOD,*) NF IF (NF.LE.1024) GO TO 30 20 WRITE (OTCD2,9997) NF 9997 FORMAT (16H FILTER LENGTH =, I4, 17H IS OUT OF BOUNDS) GO TO 15 30 WRITE (*,'(,A,/,10X,A)') ' ENTER WINDOW TYPE (1=RECTANGULAR, 2=TRI +ANGULAR,', '3=HAMMING, 4=GEN. HAMMING, 5=HANNING, 6=KAISER, 7=CHEB +YSHEV): ' READ (INCOD,*) ITYPE 35 IF (ITYPE.NE.7 .AND. NF.LT.3) GO TO 20 IF (ITYPE.EQ.7 .AND. (NF.EQ.1 .OR. NF.EQ.2)) GO TO 20 WRITE (*,'(,A,/,20X,A)') ' ENTER FILTER TYPE (1=LOWPASS, 2=HIGHPAS +S,', '3=BANDPASS, 4=BANDSTOP): ' READ (INCOD,*) JTYPE C C N IS HALF THE LENGTH OF THE SYMMETRIC FILTER C N = (NF+1)/2 IF (JTYPE.NE.1 .AND. JTYPE.NE.2) GO TO 50 C C FOR THE IDEAL LOWPASS OR HIGHPASS DESIGN - INPUT FC C 40 WRITE (*,'(,A)') ' ENTER CUTOFF FREQUENCY (AS A FRACTION OF SAMPLE + FREQUENCY): ' READ (INCOD,*) FC IF (FC.GT.0.0 .AND. FC.LT.0.5) GO TO 60 WRITE (OTCD1,9995) FC 9995 FORMAT (13H CUTOFF FREQ=, F14.7, 29H IS OUT OF BOUNDS, REENTER DA, * 2HTA) GO TO 40 C C FOR THE IDEAL BANDPASS OR BANDSTOP DESIGN - INPUT FL AND FH C 50 WRITE (*,'(,A)') ' ENTER LOWER CUTOFF FREQUENCY (AS A FRACTION OF +SAMPLE FREQUENCY): ' READ (INCOD,*) FL WRITE (*,'(,A)') ' ENTER UPPER CUTOFF FREQUENCY (AS A FRACTION OF +SAMPLE FREQUENCY): ' READ (INCOD,*) FH IF (FL.GT.0.0 .AND. FL.LT.0.5 .AND. FH.GT.0.0 .AND. FH.LT.0.5 * .AND. FH.GT.FL) GO TO 60 IF (FL.LT.0. .OR. FL.GT.0.5) WRITE (OTCD1,9995) FL IF (FH.LT.0. .OR. FH.GT.0.5) WRITE (OTCD1,9995) FH IF (FH.LT.FL) WRITE (OTCD1,9992) FH, FL 9992 FORMAT (4H FH=, F14.7, 20H IS SMALLER THAN FL=, F14.7, 8H REENTER, * 5H DATA) GO TO 50 60 IF (ITYPE.NE.7) GO TO 70 C C INPUT FOR CHEBYSHEV WINDOW--2 OF THE 3 PARAMETERS NF, DPLOG, AND DF C MUST BE SPECIFIED, WHERE DPLOG IS THE DESIRED FILTER RIPPLE(DB SCALE), C DF IS THE TRANSITION WIDTH (NORMALIZED) OF THE FILTER, C AND NF IS THE FILTER LENGTH. THE UNSPECIFIED PARAMETER C IS READ IN WITH THE ZERO VALUE. C WRITE (*,'(,A)') ' ENTER CHEBYSHEV RIPPLE (IN DB), AND OPTIONAL TR +ANSITION WIDTH,' WRITE (*,'(,A)') ' SEPARATED BY A COMMA, USING DECIMAL POINTS: ' READ (INCOD,9993) DPLOG, DF 9993 FORMAT (2F14.7) DP = 10.0**(-DPLOG/20.0) CALL CHEBC(NF, DP, DF, N, X0, XN) C C IEO IS AN EVEN, ODD INDICATOR, IEO = 0 FOR EVEN, IEO = 1 FOR ODD C 70 IEO = MOD(NF,2) C C THE FOLLOWING SEEMED WRONG: HAD 1 INSTEAD OF 4. EVEN LENGTH C LOWPASS FILTERS COULD RESULT IN THE ORIGINAL. C IF (IEO.EQ.1 .OR. JTYPE.EQ.4 .OR. JTYPE.EQ.3) GO TO 80 WRITE (*,'(,A)') ' ORDER MUST BE ODD FOR HIGHPASS OF LOWPASS FILTE +RS -- BEING INCREASED BY 1. ' NF = NF + 1 N = (1+NF)/2 IEO = 1 80 CONTINUE C C COMPUTE IDEAL (UNWINDOWED) IMPULSE RESPONSE FOR FILTER C C1 = FC IF (JTYPE.EQ.3 .OR. JTYPE.EQ.4) C1 = FH - FL IF (IEO.EQ.1) G(1) = 2.*C1 I1 = IEO + 1 DO 90 I=I1,N XN = I - 1 IF (IEO.EQ.0) XN = XN + 0.5 C = PI*XN C3 = C*C1 IF (JTYPE.EQ.1 .OR. JTYPE.EQ.2) C3 = 2.*C3 G(I) = SIN(C3)/C IF (JTYPE.EQ.3 .OR. JTYPE.EQ.4) G(I) = G(I)*2.*COS(C*(FL+FH)) 90 CONTINUE C C COMPUTE A RECTANGULAR WINDOW C IF (ITYPE.EQ.1) WRITE (OTCD2,9989) NF 9989 FORMAT (23H RECTANGULAR WINDOW-NF=, I4) DO 100 I=1,N W(I) = 1. 100 CONTINUE C C DISPATCH ON WINDOW TYPE C GO TO (200, 110, 120, 140, 150, 160, 170), ITYPE C C TRIANGULAR WINDOW C 110 CALL TRIANG(NF, W, N, IEO) WRITE (OTCD2,9988) NF 9988 FORMAT (22H TRIANGULAR WINDOW-NF=, I4) GO TO 180 C C HAMMING WINDOW C 120 ALPHA = 0.54 WRITE (OTCD2,9987) NF 9987 FORMAT (23H HAMMING WINDOW LENGTH=, I4) 130 BETA = 1. - ALPHA CALL HAMMIN(NF, W, N, IEO, ALPHA, BETA) WRITE (OTCD2,9986) ALPHA 9986 FORMAT (7H ALPHA=, F14.7) GO TO 180 C C GENERALIZED HAMMING WINDOW C FORM OF WINDOW IS W(M)=ALPHA+BETA*COS((TWOPI*M)/(NF-1)) C BETA IS AUTOMATICALLY SET TO 1.-ALPHA C READ IN ALPHA C 140 WRITE (OTCD1,9985) 9985 FORMAT (45H SPECIFY ALPHA FOR GENERALIZED HAMMING WINDOW) READ (INCOD,9993) ALPHA WRITE (OTCD2,9984) NF 9984 FORMAT (35H GENERALIZED HAMMING WINDOW LENGTH=, I4) GO TO 130 C C HANNING WINDOW C 150 ALPHA = 0.5 WRITE (OTCD2,9983) NF 9983 FORMAT (23H HANNING WINDOW LENGTH=, I4) C C INCREASE NF BY 2 AND N BY 1 FOR HANNING WINDOW SO ZERO C ENDPOINTS ARE NOT PART OF WINDOW C NF = NF + 2 N = N + 1 GO TO 130 C C KAISER (I0-SINH) WINDOW C NEED TO SPECIFY PARAMETER ATT=STOPBAND ATTENUATION IN DB C 160 WRITE (OTCD1,9982) 9982 FORMAT (33H SPECIFY ATTENUATION IN DB(F14.7), 16H FOR KAISER WIND, * 2HOW) READ (INCOD,9993) ATT IF (ATT.GT.50.) BETA = 0.1102*(ATT-8.7) IF (ATT.GE.20.96 .AND. ATT.LE.50.) BETA = 0.58417*(ATT-20.96)** * 0.4 + 0.07886*(ATT-20.96) IF (ATT.LT.20.96) BETA = 0. CALL KAISER(NF, W, N, IEO, BETA) WRITE (OTCD2,9981) NF 9981 FORMAT (22H KAISER WINDOW LENGTH=, I4) WRITE (OTCD2,9980) ATT, BETA 9980 FORMAT (14H ATTENUATION=, F14.7, 7H BETA=, F14.7) GO TO 180 C C CHEBYSHEV WINDOW C 170 CALL CHEBY(NF, W, N, IEO, DP, DF, X0, XN) WRITE (OTCD2,9979) NF 9979 FORMAT (25H CHEBYSHEV WINDOW LENGTH=, I4) WRITE (OTCD2,9978) DP, DF 9978 FORMAT (8H RIPPLE=, F14.7, 19H TRANSITION WIDTH=, F14.7) C C WINDOW IDEAL FILTER RESPONSE C CHANGE BACK NF AND N FOR HANNING WINDOW C 180 IF (ITYPE.EQ.5) NF = NF - 2 IF (ITYPE.EQ.5) N = N - 1 DO 190 I=1,N G(I) = G(I)*W(I) 190 CONTINUE C C OPTIONALLY PRINT OUT RESULTS C 200 WRITE (*,'(/,A\)') ' PRINT WINDOW VALUES? (Y/N): ' READ (*,'(A1)') ANSWER IF (ANSWER .EQ. 'N' .OR. ANSWER .EQ. 'N') GO TO 220 WRITE (OTCD2,9975) 9975 FORMAT (14H WINDOW VALUES) DO 210 I=1,N J = N + 1 - I K = NF + 1 - I WRITE (OTCD2,9974) I, W(J), K 9974 FORMAT (10X, 3H W(, I3, 2H)=, E15.8, 4H =W(, I4, 1H)) 210 CONTINUE 220 IF (JTYPE.EQ.1) WRITE (OTCD2,9973) 9973 FORMAT (26H **LOWPASS FILTER DESIGN**) IF (JTYPE.EQ.2) WRITE (OTCD2,9972) 9972 FORMAT (27H **HIGHPASS FILTER DESIGN**) IF (JTYPE.EQ.3) WRITE (OTCD2,9971) 9971 FORMAT (27H **BANDPASS FILTER DESIGN**) IF (JTYPE.EQ.4) WRITE (OTCD2,9970) 9970 FORMAT (27H **BANDSTOP FILTER DESIGN**) IF (JTYPE.EQ.1) WRITE (OTCD2,9969) FC 9969 FORMAT (22H IDEAL LOWPASS CUTOFF=, F14.7) IF (JTYPE.EQ.2) WRITE (OTCD2,9968) FC 9968 FORMAT (23H IDEAL HIGHPASS CUTOFF=, F14.7) IF (JTYPE.EQ.3 .OR. JTYPE.EQ.4) WRITE (OTCD2,9967) FL, FH 9967 FORMAT (26H IDEAL CUTOFF FREQUENCIES=, 2F14.7) IF (JTYPE.EQ.1 .OR. JTYPE.EQ.3) GO TO 240 DO 230 I=2,N G(I) = -G(I) 230 CONTINUE G(1) = 1.0 - G(1) C C WRITE OUT IMPULSE RESPONSE C C C 240 DO 250 I=1,N J = N + 1 - I K = NF + 1 - I WRITE (OTCD2,9966) I, G(J), K 9966 FORMAT (10X, 3H H(, I3, 2H)=, E15.8, 4H =H(, I4, 1H)) 250 CONTINUE CALL FLCHAR(NF, ITYPE, JTYPE, FC, FL, FH, N, IEO, G, OTCD2) WRITE (OTCD2,9965) 9965 FORMAT (1H /) C C OPEN OUTPUT FILE C IOUTD = 1 WRITE (*,'(,A)') ' ENTER NAME OF WINDOW VALUES OUTPUT FILE OR ' WRITE (*,'(,A)') ' (SORRY, NO TILDE-EXPANSION. GIVE PATH RELATIVE + TO YOUR HOME DIRECTORY) ' READ (INCOD,'(A)') FNAME IF (FNAME .EQ. ' ') GO TO 260 OPEN (IOUTD, FILE=FNAME) CALL OUTIMP(G,N,NF,IEO,IOUTD) CALL FLUSH(IOUTD) CLOSE(IOUTD) 260 CLOSE(3) WRITE (*,'(/,A\)') ' ANOTHER DESIGN? (Y/N): ' READ (*,'(A1)') ANSWER IF (ANSWER .EQ. 'Y' .OR. ANSWER .EQ. 'Y') GO TO 10 STOP 'END OF PROGRAM' END C C----------------------------------------------------------------------- C SUBROUTINE: TRIANG C TRIANGULAR WINDOW C----------------------------------------------------------------------- C SUBROUTINE TRIANG(NF, W, N, IEO) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW COEFFICIENTS FOR HALF THE WINDOW C N = HALF WINDOW LENGTH=(NF+1)/2 C IEO = EVEN - ODD INDICATION--IEO=0 FOR NF EVEN C DIMENSION W(1) FN = N DO 10 I=1,N XI = I - 1 IF (IEO.EQ.0) XI = XI + 0.5 W(I) = 1. - XI/FN 10 CONTINUE RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: HAMMIN C GENERALIZED HAMMING WINDOW ROUTINE C WINDOW IS W(N) = ALPHA + BETA * COS( TWOPI*(N-1) / (NF-1) ) C----------------------------------------------------------------------- C SUBROUTINE HAMMIN(NF, W, N, IEO, ALPHA, BETA) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW ARRAY OF SIZE N C N = HALF LENGTH OF FILTER=(NF+1)/2 C IEO = EVEN ODD INDICATOR--IEO=0 IF NF EVEN C ALPHA = CONSTANT OF WINDOW C BETA = CONSTANT OF WINDOW--GENERALLY BETA=1-ALPHA C DIMENSION W(1) PI2 = 8.0*ATAN(1.0) FN = NF - 1 DO 10 I=1,N FI = I - 1 IF (IEO.EQ.0) FI = FI + 0.5 W(I) = ALPHA + BETA*COS((PI2*FI)/FN) 10 CONTINUE RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: KAISER C KAISER WINDOW C----------------------------------------------------------------------- C SUBROUTINE KAISER(NF, W, N, IEO, BETA) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW ARRAY OF SIZE N C N = FILTER HALF LENGTH=(NF+1)/2 C IEO = EVEN ODD INDICATOR--IEO=0 IF NF EVEN C BETA = PARAMETER OF KAISER WINDOW C DIMENSION W(1) REAL INO BES = INO(BETA) XIND = FLOAT(NF-1)*FLOAT(NF-1) DO 10 I=1,N XI = I - 1 IF (IEO.EQ.0) XI = XI + 0.5 XI = 4.*XI*XI W(I) = INO(BETA*SQRT(1.-XI/XIND)) W(I) = W(I)/BES 10 CONTINUE RETURN END C C----------------------------------------------------------------------- C FUNCTION: INO C BESSEL FUNCTION FOR KAISER WINDOW C----------------------------------------------------------------------- C REAL FUNCTION INO(X) Y = X/2. T = 1.E-08 E = 1. DE = 1. DO 10 I=1,25 XI = I DE = DE*Y/XI SDE = DE*DE E = E + SDE IF (E*T-SDE) 10, 10, 20 10 CONTINUE 20 INO = E RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: CHEBC C SUBROUTINE TO GENERATE CHEBYSHEV WINDOW PARAMETERS WHEN C ONE OF THE THREE PARAMETERS NF,DP AND DF IS UNSPECIFIED C----------------------------------------------------------------------- C SUBROUTINE CHEBC(NF, DP, DF, N, X0, XN) C C NF = FILTER LENGTH (IN SAMPLES) C DP = FILTER RIPPLE (ABSOLUTE SCALE) C DF = NORMALIZED TRANSITION WIDTH OF FILTER C N = (NF+1)/2 = FILTER HALF LENGTH C X0 = (3-C0)/(1+C0) WITH C0=COS(PI*DF) = CHEBYSHEV WINDOW CONSTANT C XN = NF-1 C PI = 4.*ATAN(1.0) IF (NF.NE.0) GO TO 10 C C DP,DF SPECIFIED, DETERMINE NF C C1 = COSHIN((1.+DP)/DP) C0 = COS(PI*DF) X = 1. + C1/COSHIN(1./C0) C C INCREMENT BY 1 TO GIVE NF WHICH MEETS OR EXCEEDS SPECS ON DP AND DF C NF = X + 1.0 N = (NF+1)/2 XN = NF - 1 GO TO 30 10 IF (DF.NE.0.0) GO TO 20 C C NF,DP SPECIFIED, DETERMINE DF C XN = NF - 1 C1 = COSHIN((1.+DP)/DP) C2 = COSH(C1/XN) DF = ARCCOS(1./C2)/PI GO TO 30 C C NF,DF SPECIFIED, DETERMINE DP C 20 XN = NF - 1 C0 = COS(PI*DF) C1 = XN*COSHIN(1./C0) DP = 1./(COSH(C1)-1.) 30 X0 = (3.-COS(2.*PI*DF))/(1.+COS(2.*PI*DF)) RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: CHEBY C DOLPH CHEBYSHEV WINDOW DESIGN C----------------------------------------------------------------------- C SUBROUTINE CHEBY(NF, W, N, IEO, DP, DF, X0, XN) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW ARRAY OF SIZE N C N = HALF LENGTH OF FILTER = (NF+1)/2 C IEO = EVEN-ODD INDICATOR--IEO=0 FOR NF EVEN C DP = WINDOW RIPPLE ON AN ABSOLUTE SCALE C DF = NORMALIZED TRANSITION WIDTH OF WINDOW C X0 = WINDOW PARAMETER RELATED TO TRANSITION WIDTH C XN = NF-1 C DIMENSION W(1) DIMENSION PR(1024), PI(1024) PIE = 4.*ATAN(1.0) XN = NF - 1 FNF = NF ALPHA = (X0+1.)/2. BETA = (X0-1.)/2. TWOPI = 2.*PIE C2 = XN/2. DO 40 I=1,NF XI = I - 1 F = XI/FNF X = ALPHA*COS(TWOPI*F) + BETA IF (ABS(X)-1.) 10, 10, 20 10 P = DP*COS(C2*ARCCOS(X)) GO TO 30 20 P = DP*COSH(C2*COSHIN(X)) 30 PI(I) = 0. PR(I) = P C C FOR EVEN LENGTH FILTERS USE A ONE-HALF SAMPLE DELAY C ALSO THE FREQUENCY RESPONSE IS ANTISYMMETRIC IN FREQUENCY C IF (IEO.EQ.1) GO TO 40 PR(I) = P*COS(PIE*F) PI(I) = -P*SIN(PIE*F) IF (I.GT.(NF/2+1)) PR(I) = -PR(I) IF (I.GT.(NF/2+1)) PI(I) = -PI(I) 40 CONTINUE C C USE DFT TO GIVE WINDOW C TWN = TWOPI/FNF DO 60 I=1,N XI = I - 1 SUM = 0. DO 50 J=1,NF XJ = J - 1 SUM = SUM + PR(J)*COS(TWN*XJ*XI) + PI(J)*SIN(TWN*XJ*XI) 50 CONTINUE W(I) = SUM 60 CONTINUE C1 = W(1) DO 70 I=1,N W(I) = W(I)/C1 70 CONTINUE RETURN END C C----------------------------------------------------------------------- C FUNCTION: COSHIN C FUNCTION FOR HYPERBOLIC INVERSE COSINE OF X C----------------------------------------------------------------------- C REAL FUNCTION COSHIN(X) COSHIN = ALOG(X+SQRT(X*X-1.)) RETURN END C C----------------------------------------------------------------------- C FUNCTION: ARCCOS C FUNCTION FOR INVERSE COSINE OF X C----------------------------------------------------------------------- C FUNCTION ARCCOS(X) IF (X) 30, 20, 10 10 A = SQRT(1.-X*X)/X ARCCOS = ATAN(A) RETURN 20 ARCCOS = 2.*ATAN(1.0) RETURN 30 A = SQRT(1.-X*X)/X ARCCOS = ATAN(A) + 4.*ATAN(1.0) RETURN END C C----------------------------------------------------------------------- C FUNCTION: COSH C FUNCTION FOR HYPERBOLIC COSINE OF X C----------------------------------------------------------------------- C REAL FUNCTION COSH(X) COSH = (EXP(X)+EXP(-X))/2. RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: FLCHAR C SUBROUTINE TO DETERMINE FILTER CHARACTERISTICS C----------------------------------------------------------------------- C SUBROUTINE FLCHAR(NF, ITYPE, JTYPE, FC, FL, FH, N, IEO, G, OTCD2) C C NF = FILTER LENGTH IN SAMPLES C ITYPE = WINDOW TYPE C JTYPE = FILTER TYPE C FC = IDEAL CUTOFF OF LP OR HP FILTER C FL = LOWER CUTOFF OF BP OR BS FILTER C FH = UPPER CUTOFF OF BP OR BS FILTER C N = FILTER HALF LENGTH = (NF+1) / 2 C IEO = EVEN ODD INDICATOR C G = FILTER ARRAY OF SIZE N C OTCD2 = OUTPUT CODE FOR LINE PRINTER USED IN WRITE STATEMENTS C DIMENSION G(1) DIMENSION RESP(2048) INTEGER OTCD2 C C NOT FOR FOR TRIANGULAR WINDOW C IF (ITYPE.EQ.2) RETURN C C DFT TO GET FREQ RESP C PI = 4.*ATAN(1.0) C C UP TO 4096 PT DFT C NR = 8*NF IF (NR.GT.2048) NR = 2048 XNR = NR TWN = PI/XNR SUMI = -G(1)/2. IF (IEO.EQ.0) SUMI = 0. DO 20 I=1,NR XI = I - 1 TWNI = TWN*XI SUM = SUMI DO 10 J=1,N XJ = J - 1 IF (IEO.EQ.0) XJ = XJ + .5 SUM = SUM + G(J)*COS(XJ*TWNI) 10 CONTINUE RESP(I) = 2.*SUM 20 CONTINUE C C DISPATCH ON FILTER TYPE C GO TO (30, 40, 50, 60), JTYPE C C LOWPASS C 30 CALL RIPPLE(NR, 1., 0., FC, RESP, F1, F2, DB) WRITE (OTCD2,9999) F2, DB 9999 FORMAT (17H PASSBAND CUTOFF , F6.4, 9H RIPPLE , F8.3, 3H DB) CALL RIPPLE(NR, 0., FC, .5, RESP, F1, F2, DB) WRITE (OTCD2,9998) F1, DB 9998 FORMAT (17H STOPBAND CUTOFF , F6.4, 9H RIPPLE , F8.3, 3H DB) RETURN C C HIGHPASS C 40 CALL RIPPLE(NR, 0., 0., FC, RESP, F1, F2, DB) WRITE (OTCD2,9998) F2, DB CALL RIPPLE(NR, 1., FC, .5, RESP, F1, F2, DB) WRITE (OTCD2,9999) F1, DB RETURN C C BANDPASS C 50 CALL RIPPLE(NR, 0., 0., FL, RESP, F1, F2, DB) WRITE (OTCD2,9998) F2, DB CALL RIPPLE(NR, 1., FL, FH, RESP, F1, F2, DB) WRITE (OTCD2,9997) F1, F2, DB 9997 FORMAT (18H PASSBAND CUTOFFS , F6.4, 2X, F6.4, 8H RIPPLE, F9.3, * 3H DB) CALL RIPPLE(NR, 0., FH, .5, RESP, F1, F2, DB) WRITE (OTCD2,9998) F1, DB RETURN C C STOPBAND C 60 CALL RIPPLE(NR, 1., 0., FL, RESP, F1, F2, DB) WRITE (OTCD2,9999) F2, DB CALL RIPPLE(NR, 0., FL, FH, RESP, F1, F2, DB) WRITE (OTCD2,9996) F1, F2, DB 9996 FORMAT (18H STOPBAND CUTOFFS , F6.4, 2X, F6.4, 8H RIPPLE, F9.3, * 3H DB) CALL RIPPLE(NR, 1., FH, .5, RESP, F1, F2, DB) WRITE (OTCD2,9999) F1, DB RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: RIPPLE C FINDS LARGEST RIPPLE IN BAND AND LOCATES BAND EDGES BASED ON THE C POINT WHERE THE TRANSITION REGION CROSSES THE MEASURED RIPPLE BOUND C----------------------------------------------------------------------- C SUBROUTINE RIPPLE(NR, RIDEAL, FLOW, FHI, RESP, F1, F2, DB) C C NR = SIZE OF RESP C RIDEAL = IDEAL FREQUENCY RESPONSE C FLOW = LOW EDGE OF IDEAL BAND C FHI = HIGH EDGE OF IDEAL BAND C RESP = FREQUENCY RESPONSE OF SIZE NR C F1 = COMPUTED LOWER BAND EDGE C F2 = COMPUTED UPPER BAND EDGE C DB = DEVIATION FROM IDEAL RESPONSE IN DB C DIMENSION RESP(1) XNR = NR C C BAND LIMITS C IFLOW = 2.*XNR*FLOW + 1.5 IFHI = 2.*XNR*FHI + 1.5 IF (IFLOW.EQ.0) IFLOW = 1 IF (IFHI.GE.NR) IFHI = NR - 1 C C FIND MAX AND MIN PEAKS IN BAND C RMIN = RIDEAL RMAX = RIDEAL DO 20 I=IFLOW,IFHI IF (RESP(I).LE.RMAX .OR. RESP(I).LT.RESP(I-1) .OR. * RESP(I).LT.RESP(I+1)) GO TO 10 RMAX = RESP(I) 10 IF (RESP(I).GE.RMIN .OR. RESP(I).GT.RESP(I-1) .OR. * RESP(I).GT.RESP(I+1)) GO TO 20 RMIN = RESP(I) 20 CONTINUE C C PEAK DEVIATION FROM IDEAL C RIPL = AMAX1(RMAX-RIDEAL,RIDEAL-RMIN) C C SEARCH FOR LOWER BAND EDGE C F1 = FLOW IF (FLOW.EQ.0.0) GO TO 50 DO 30 I=IFLOW,IFHI IF (ABS(RESP(I)-RIDEAL).LE.RIPL) GO TO 40 30 CONTINUE 40 XI = I - 1 C C LINEAR INTERPOLATION OF BAND EDGE FREQUENCY TO IMPROVE ACCURACY C X1 = .5*XI/XNR X0 = .5*(XI-1.)/XNR Y1 = ABS(RESP(I)-RIDEAL) Y0 = ABS(RESP(I-1)-RIDEAL) F1 = (X1-X0)/(Y1-Y0)*(RIPL-Y0) + X0 C C SEARCH FOR UPPER BAND EDGE C 50 F2 = FHI IF (FHI.EQ.0.5) GO TO 80 DO 60 I=IFLOW,IFHI J = IFHI + IFLOW - I IF (ABS(RESP(J)-RIDEAL).LE.RIPL) GO TO 70 60 CONTINUE 70 XI = J - 1 C C LINEAR INTERPOLATION OF BAND EDGE FREQUENCY TO IMPROVE ACCURACY C X1 = .5*XI/XNR X0 = .5*(XI+1.)/XNR Y1 = ABS(RESP(J)-RIDEAL) Y0 = ABS(RESP(J+1)-RIDEAL) F2 = (X1-X0)/(Y1-Y0)*(RIPL-Y0) + X0 C C DEVIATION FROM IDEAL IN DB C 80 DB = 20.*ALOG10(RIPL+RIDEAL) RETURN END C C----------------------------------------------------------------------- C C----------------------------------------------------------------------- C SUBROUTINE: OUTIMP C PRINTS OUT ON UNIT IOUTD THE FULL C IMPULSE RESPONSE AS FORMAT 13X,F15.8. THIS MAKES THE IMPULSE C RESPONSE ACCESSIBLE TO OTHER PROGRAMS WITHOUT THE NEED FOR PARSING C THE OUTPUT OF THE PROGRAM. C----------------------------------------------------------------------- C SUBROUTINE OUTIMP(G,N,NF,IEO,IOUTD) DIMENSION G(512) C C REMOVED FILTER ORDER OUTPUT FOR GABRIEL COMPATIBILITY C C WRITE(IOUTD,5) NF C 5 FORMAT(I15) L = N IF(IEO.EQ.0) GO TO 10 L = N-1 10 DO 20 I=1,L J = N + 1 - I WRITE(IOUTD,50) G(J) 20 CONTINUE IF(IEO.EQ.0) GO TO 30 WRITE(IOUTD,50) G(1) 30 DO 40 I=L,1,-1 J = N + 1 - I WRITE(IOUTD,50) G(J) 40 CONTINUE 50 FORMAT(E15.8) CLOSE(UNIT=IOUTD) 60 RETURN END