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 <ENTER
     +> 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 <ENT
     +ER> '
      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
