FIR Windowed Filter Design Program
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: WINDOW DESIGN OF LINEAR PHASE, LOWPASS, HIGHPASS
C BANDPASS, AND BANDSTOP FIR DIGITAL FILTERS
C AUTHOR: LAWRENCE R. RABINER AND CAROL A. MCGONEGAL
C BELL LABORATORIES, MURRAY HILL, NEW JERSEY, 07974
C MODIFIED JAN. 1978 BY DOUG PAUL, MIT LINCOLN LABORATORIES
C TO INCLUDE SUBROUTINES FOR OBTAINING FILTER BAND EDGES AND RIPPLES
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 FL AND FH ARE THE NORMALIZED FILTER CUTOFF FREQUENCIES
C 0 <= FL <= FH <= 0.5
C IWP OPTIONALLY PRINTS OUT THE WINDOW VALUES
C IWP = 0 DO NOT PRINT
C IWP = 1 PRINT
C IMD REQUESTS ADDITIONAL RUNS
C IMD = 1 NEW RUN
C IMD = 0 TERMINATES PROGRAM
C-----------------------------------------------------------------------
C
DIMENSION W(512), G(512)
INTEGER OTCD1, OTCD2
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: ALL OUTPUT IS WRITTEN ON THE STANDARD
C OUTPUT UNIT (OTCD2)
C
INCOD = I1MACH(1)
OTCD1 = I1MACH(4)
OTCD2 = I1MACH(2)
C
C INPUT THE FILTER LENGTH(NF), WINDOW TYPE(ITYPE) AND FILTER TYPE(JTYPE)
C
10 WRITE (OTCD1,9999)
9999 FORMAT (44H SPECIFY FILTER LENGTH(I4), WINDOW TYPE(I2),, 6H FILTE,
* 10HR TYPE(I2))
READ (INCOD,9998) NF, ITYPE, JTYPE
9998 FORMAT (I4, I2, I2)
IF (NF.LE.1024) GO TO 30
20 WRITE (OTCD2,9997) NF
9997 FORMAT (4H NF=, I4, 17H IS OUT OF BOUNDS)
GO TO 10
30 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
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 (OTCD1,9996)
9996 FORMAT (38H SPECIFY IDEAL CUTOFF FREQUENCY(F14.7))
READ (INCOD,9993) FC
IF (FC.GT.0.0 .AND. FC.LT.0.5) GO TO 60
WRITE (OTCD1,9995) FC
9995 FORMAT (13H VALUE OF FC=, 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 (OTCD1,9994)
9994 FORMAT (43H SPECIFY LOWER AND UPPER CUTOFF FREQUENCIES, 7H(2F14.7,
* 1H))
READ (INCOD,9993) FL, FH
9993 FORMAT (2F14.7)
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 (OTCD1,9991)
9991 FORMAT (46H SPECIFY CHEBYSHEV RIPPLE IN DB (F14.7) AND/OR,
*24H TRANSITION WIDTH(F14.7))
READ (INCOD,9993) DPLOG, DF
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)
IF (IEO.EQ.1 .OR. JTYPE.EQ.1 .OR. JTYPE.EQ.3) GO TO 80
WRITE (OTCD1,9990)
9990 FORMAT (48H NF MUST BE ODD INTEGER FOR HP OR BS FILTERS--NF,
* 24H IS 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 (19H HAMMING WINDOW-NF=, 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 (31H GENERALIZED HAMMING WINDOW-NF=, I4)
GO TO 130
C
C HANNING WINDOW
C
150 ALPHA = 0.5
WRITE (OTCD2,9983) NF
9983 FORMAT (19H HANNING WINDOW-NF=, 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 (18H KAISER WINDOW-NF=, I4)
WRITE (OTCD2,9980) ATT, BETA
9980 FORMAT (6H ATT=, 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 (21H CHEBYSHEV WINDOW-NF=, I4)
WRITE (OTCD2,9978) DP, DF
9978 FORMAT (4H DP=, F14.7, 5H DF=, 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 PRINT OUT RESULTS
C
200 WRITE (OTCD1,9977)
9977 FORMAT (36H PRINT OUT WINDOW VALUES(1=YES,0=NO))
READ (INCOD,9976) IWP
9976 FORMAT (I1)
IF (IWP.EQ.0) 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
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 /1H /1H /1H )
WRITE (OTCD2,9964)
9964 FORMAT (1H1)
WRITE (OTCD1,9963)
9963 FORMAT (26H MORE DATA(1=YES,0=NO)(I1))
READ (INCOD,9962) IMD
9962 FORMAT (I1)
IF (IMD.EQ.1) GO TO 10
STOP
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
Return to Main Software Page