A Subroutine for Finite Wordlength FIR Filter Design

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM:  IDEFIX - TEST PROGRAM FOR THE SUBROUTINE IDEFIR
C
C AUTHOR:  ULRICH HEUTE
C INSTITUT FUER NACHRICHTENTECHNIK,UNIV.ERLANGEN-NUERNBERG
C D-8520 ERLANGEN, W-GERMANY
C INPUT:
C     FORMAT...6F8.6 (1 CARD).
C     DEL1....PASSBAND TOLERANCE;
C     DEL2....STOPBAND TOLERANCE;
C     OMGP1...FIRST PASSBAND EDGE FREQUENCY;
C     OMGS1...FIRST STOPBAND EDGE FREQUENCY;
C     OMGP2...SECOND PASSBAND EDGE FREQUENCY;
C     OMGS2...SECOND STOPBAND EDGE FREQUENCY.
C
C     NORMALIZED FREQUENCY OMG = F / SAMPLING RATE.
C     IN LOWPASS AND HIGHPASS CASES, THE LAST TWO PARAMETERS ARE ZERO
C     IN BANDFILTER CASES, THE LAST TWO PARAMETERS HAVE THE VALUES
C               OMGP(S)2 = 0.5 - OMGP(S)1
C     IN HILBERT CASES, DEL2 = DEL1 IS VALID;OMGP1 IS THE HILBERT BAND
C     EDGE FREQUENCY; THE LAST T H R E E PARAMETERS ARE ZERO.
C
C OUTPUT:
C     ALL DATA CONCERNING THE DESIGN ARE DELIVERED TO THE STANDARD
C     OUTPUT UNIT.
C-----------------------------------------------------------------------
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
      COMMON /CHEK/ NL, NG, KA, KB
C
C THE COMMON-DEFINITIONS HAVE TO BE ADAPTED TO THE REQUIREMENTS OF
C THE SUBROUTINE 'DESIGN' CALLED BY THE SUBROUTINE 'IDEFIR'.
C
      DIMENSION H0(256), RH0(128)
C
C SET UP MACHINE CONSTANTS.
C
C IRCHN...CHANNEL ASSIGNED TO INPUT DEVICE;
C IWCHN...CHANNEL ASSIGNED TO OUTPUT DEVICE;
C
      IRCHN = I1MACH(1)
      IWCHN = I1MACH(2)
C
C INPUT OF THE TOLERANCE SCHEME.
C
      READ (IRCHN,9999) DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2
C
C A MAXIMUM OF 20 DESIGN ITERATIONS IS ALLOWED.
C
      ITMAX = 20
C
      DP = DEL1
      OP = OMGP1
      OS = OMGS1
C
      WRITE (IWCHN,9998)
C
      WRITE (IWCHN,9997)
C
      CALL IDEFIR(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP, NF00,
     *    NF, IWCMIN, IWC, IWACT, DEL11, DEL21, EMAX1, EMAX2, H0, RH0,
     *    ITERA, IERR)
C
C OUTPUT OF FILTER TYPE AND DESIGN PARAMETERS.
C
      IF (ITYP.EQ.0) WRITE (IWCHN,9995) DEL1, DEL2, OP, OS
      IF (ITYP.EQ.1) WRITE (IWCHN,9994) DEL1, DEL2, OP, OS
      IF (ITYP.EQ.2) WRITE (IWCHN,9993) DEL1, DEL2, OP, OS, OMGP2, OMGS2
      IF (ITYP.EQ.3) WRITE (IWCHN,9992) DEL1, DEL2, OP, OS, OMGP2, OMGS2
      IF (ITYP.EQ.4) WRITE (IWCHN,9991) DP, OP
C
C TEST FOR POSSIBLE ERROR RETURNS.
C
      GO TO (40, 10, 20, 30), IERR
C
C OUTPUT FOR THE CASE OF INCONSISTENT SPECIFICATIONS.
C
  10  WRITE (IWCHN,9981)
      STOP
C
C OUTPUT, IF THE MAXIMUM ALLOWED FILTER ORDER IS EXCEEDED.
C
  20  WRITE (IWCHN,9980)
      STOP
C
C OUTPUT,IF NO SOLUTION WAS FOUND.
C
  30  WRITE (IWCHN,9982) ITERA
  40  CONTINUE
C
C OUTPUT OF REALIZATION PARAMETERS.
C
      IF (ITYP.NE.0) WRITE (IWCHN,9990) DEL1, DEL2, OMGP1, OMGS1
      WRITE (IWCHN,9989) IWCMIN, NF00
      WRITE (IWCHN,9996) IWC, IWACT, NF
C
C OUTPUT OF THE REDUCED TOLERANCE SCHEME.
C
      WRITE (IWCHN,9988) DEL11, DEL21
C
C OUTPUT OF EXACT AND ROUNDED COEFFICIENTS.
C
      N1 = (NF+1)/2
      WRITE (IWCHN,9985)
      WRITE (IWCHN,9987) (H0(I),I=1,N1)
      WRITE (IWCHN,9984)
      WRITE (IWCHN,9987) (RH0(I),I=1,N1)
C
C OUTPUT OF ACTUAL LOWPASS TOLERANCES.
C
      WRITE (IWCHN,9986) EMAX1, EMAX2
      WRITE (IWCHN,9983) ITERA
C
      STOP
C
9999  FORMAT (6F8.6)
C
9998  FORMAT (1H1)
9997  FORMAT (54H ITERATIVE DESIGN OF LINEAR PHASE FIR FILTERS IN DIREC,
     *    6HT FORM/36H WITH MINIMUM COEFFICIENT WORDLENGTH//)
9996  FORMAT (41H COEFFICIENT WORDLENGTH APPLIED     WC = , I2/6H ACTUA,
     *    48HL WORDLENGTH WITHOUT LEADING ZEROS     WC ACT = ,
     *    I2/31H ACTUAL FILTER LENGTH     NF = , I3/)
9995  FORMAT (24H LOWPASS WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
     *    F8.6/14X, 10HOMEGA P = , F8.6, 11H,OMEGA S = , F8.6/)
9994  FORMAT (25H HIGHPASS WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
     *    F8.6/15X, 10HOMEGA P = , F8.6, 11H,OMEGA S = , F8.6/)
9993  FORMAT (25H BANDPASS WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
     *    F8.6/15X, 10HOMEGAP1 = , F8.6, 11H,OMEGAS1 = , F8.6/15X,
     *    10HOMEGAP2 = , F8.6, 11H,OMEGAS2 = , F8.6/)
9992  FORMAT (25H BANDSTOP WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
     *    F8.6/15X, 10HOMEGAP1 = , F8.6, 11H,OMEGAS1 = , F8.6/15X,
     *    10HOMEGAP2 = , F8.6, 11H,OMEGAS2 = , F8.6/)
9991  FORMAT (34H HILBERT-TRANSFORMER WITH DELTA = , F8.6, 9H,OMEGAP =,
     *    1H , F8.6/)
9990  FORMAT (35H EQUIVALENT LOWPASS WITH DELTA 1 = , F8.6, 9H,DELTA 2 ,
     *    2H= , F8.6/25X, 10HOMEGA P = , F8.6, 11H,OMEGA S = , F8.6/)
9989  FORMAT (31H EST. MIN. WORDLENGTH WC MIN = , I2,      9H BITS (SI,
     *    30HGN+BITS BEHIND THE DUAL POINT)/  22H ESTIMATED LOWPASS LEN,
     *    9HGTH NF = , I3/)
9988  FORMAT (54H REDUCED LOWPASS TOLERANCE SCHEME     DELTA 1 PRIME = ,
     *    F8.6/ 54H                                      DELTA 2 PRIME =
     *    , F8.6/)
9987  FORMAT (4(1X, F13.6))
9986  FORMAT (1H /50H ACTUAL TOLERANCES OF THE EQUIVALENT LOWPASS WITH ,
     *    21HROUNDED COEFFICIENTS /11H DELTA 1 = , F8.6, 11H,DELTA 2 = ,
     *    F8.6/)
9985  FORMAT (54H EXACT FILTER COEFFICIENTS H0 (I), I = 0...((NF-1)/2)
     *    /)
9984  FORMAT (1H /37H ROUNDED FILTER COEFFICIENTS H0 (I)  /)
9983  FORMAT (31H NUMBER OF DESIGN ITERATIONS = , I3)
9982  FORMAT (34H NO SOLUTION FOUND WITHIN ALLOWED , I2, 11H DESIGN ITE,
     *    8HRATIONS-/47H IN THE LAST TRIAL, THE FOLLOWING PARAMETERS WE,
     *    9HRE USED /)
9981  FORMAT (1H ///48H INCONSISTENT SPECIFICATIONS,CHECK PARAMETER ORD,
     *    6HERING )
9980  FORMAT (1H ///48H FILTER REQUIREMENTS EXCEED MAX. ALLOWED FILTER ,
     *    6HORDER )
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: I D E F I R
C
C     DESIGN OF LINEAR-PHASE FIR FILTERS IN DIRECT FORM WITH MINIMUM
C     COEFFICIENT WORDLENGTH. POSSIBLE FILTER TYPES INCLUDE LOWPASS
C     HIGHPASS,SYMMETRICAL BANDPASS AND BANDSTOP FILTERS (WITH CENTER
C     FREQUENCY AT OMEGA = 0.25, WHERE OMEGA = 0.5 IS HALF THE SAMPLING
C     FREQUENCY), AND HILBERT TRANSFORMERS.
C     THE DESIGN IS CARRIED OUT IN THE FOLLOWING WAY - -
C     PROBLEM TRANSFORMATION TO AN EQUIVALENT LOWPASS PROBLEM;
C     FILTER LENGTH,WORDLENGTH, AND ROUNDING ERROR ESTIMATION;
C     LOWPASS DESIGN WITH APPROPRIATELY REDUCED TOLERANCES;
C     COEFFICIENT ROUNDING, TOLERANCE CONTROL; ITERATION IF NECESSARY;
C     COEFFICIENT TRANSFORMATION BACK TO THE DESIRED FILTER TYPE.
C
C  I N P U T
C     I N P U T  P A R A M E T E R S = TOLERANCE SCHEME WITH
C     DEL1...PASSBAND TOLERANCE DELTA 1;
C     DEL2...STOPBAND TOLERANCE DELTA 2;
C     OMGP1..FIRST PASSBAND EDGE FREQUENCY;
C     OMGS1..FIRST STOPBAND EDGE FREQUENCY;
C     OMGP2..SECOND PASSBAND EDGE FREQUENCY;
C     OMGS2..SECOND STOPBAND EDGE FREQUENCY.
C     NORMALIZED FREQUENCY OMG = F / SAMPLING RATE.
C     IN LOWPASS AND HIGHPASS CASES,THE LAST TWO PARAMETERS ARE ZERO.
C     IN BANDFILTER CASES,THE LAST TWO PARAMETERS HAVE THE VALUES
C                OMGP(S)2 = 0.5 - OMGP(S)1
C     IN HILBERT CASES, DEL2 = DEL1 IS VALID;OMGP1 IS THE HILBERT BAND
C     EDGE FREQUENCY; THE LAST T H R E E  PARAMETERS ARE ZERO.
C
C  O U T P U T
C     O U T P U T  P A R A M E T E R S
C     DEL1,DEL2,OMGP1,OMGS1...TOLERANCE SCHEME OF THE EQUIVALENT LOWPASS
C     PROBLEM; THE ORIGINAL TOLERANCE SCHEME IS REPLACED .
C     ITYP...FILTER TYPE INFORMATION- ITYP =0..LOWPASS;
C                                           1..HIGHPASS;
C                                           2..BANDPASS;
C                                           3..BANDSTOP;
C                                           4..HILBERT-TRANSFORMER;
C     NF00...FILTER LENGTH ACCORDING TO TOLERANCES AND ERROR ESTIMATION;
C     NF.....FILTER LENGTH ACTUALLY APPLIED;
C     IWCMIN...ESTIMATED MIN.COEFF.WORDLENGTH (ACCORDING TO MAX.QUANTI-
C              ZATION STEPSIZE + SIGN BIT);
C     IWC......WORDLENGTH USED FOR THE DESIGN (IN THE SAME SENSE);
C     IWACT....ACTUAL WORDLENGTH (SIGN + NECESSARY BITS BEHIND THE DUAL
C              POINT (WITHOUT LEADING ZEROS) FOR THE LARGEST ABS.VALUE);
C     DEL11,DEL21...REDUCED TOLERANCES FOR THE EQUIVALENT LOWPASS DESIGN
C     EMAX1,2..TOLERANCES (I.E. MAX. ERRORS) OF THE EQUIVALENT LOWPASS
C              WITH ROUNDED COEFFICIENTS;
C     H0.......COEFFICIENTS BEFORE ROUNDING;
C     RH0......ROUNDED COEFFICIENTS;
C     ITERA....NUMBER OF DESIGN ITERATIONS;
C     IERR.....ERROR PARAMETER -
C              IERR= 1...NO ERROR;
C                    2...INCONSISTENT SPECIFICATIONS;
C                    3...MAX. ALLOWED FILTER ORDER EXCEEDED;
C                    4...MAX. ALLOWED NUMBER OF ITERATIONS EXCEEDED
C-----------------------------------------------------------------------
C
      SUBROUTINE IDEFIR(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP,
     *    NF00, NF, IWCMIN, IWC, IWACT, DEL11, DEL21, EMAX1, EMAX2, H0,
     *    RH0, ITERA, IERR)
C
C COMMON-DEFINITIONS HAVE TO BE INSERTED HERE ACCORDING TO THE SUB-
C ROUTINE 'DESIGN' CALLED BY THIS SUBROUTINE.
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
      COMMON /CHEK/ NL, NG, KA, KB
C
      DIMENSION H0(256), RH0(128)
      IERR = 1
      SMALL = 2.*R1MACH(3)
C
C THE FOLLOWING PARAMETERS CONCERN THE MAX. ALLOWED NUMBER OF DESIGN
C ITERATIONS (ITMAX), THE MINIMUM USAGE OF THE TOLERANCES (PROZ),THE
C REDUCTION (CR) AND INCREASE (CI) FACTORS FOR THE TOLERANCE CHANGE.
C
      ITMAX = 20
      PROZ = 0.9
      CR = 0.75
      CI = 0.25
C
C TRANSFORMATION TO LOWPASS PROBLEM.
C
      CALL TYPTRA(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP, IERR)
C
C IERR = 2 IN THE CASE OF INCONSISTENT SPECIFICATIONS.
C
      IF (IERR.EQ.2) RETURN
C
      DELOMG = OMGS1 - OMGP1
C
C CHECKING THE FILTER SYMMETRY.
C
      OM = ABS(OMGP1+OMGS1-.5)
      ISYM = 0
      IF ((DEL1.EQ.DEL2) .AND. (OM.LT.SMALL)) ISYM = 1
      ITERA = 0
      ISOLUT = 0
C
C DETERMINATION OF MINIMUM COEFFICIENT WORDLENGTH.
C
      CALL WCMIN(DEL1, DEL2, DELOMG, IWCMIN)
      IWC = IWCMIN
C
C O U T E R  D E S I G N  I T E R A T I O N  L O O P
C
  10  CONTINUE
      DEL11 = DEL1
      DEL21 = DEL2
C
C I N N E R  D E S I G N  I T E R A T I O N  L O O P
C
  20  CONTINUE
C
C STOP AFTER MAXIMUM NUMBER OF ITERATIONS.
C
      IF (ITERA.LT.ITMAX) GO TO 30
      IERR = 4
      GO TO 80
  30  CONTINUE
C
C DETERMINATION OF FILTER LENGTH.
C
      CALL LENGTH(DEL11, DEL21, DELOMG, NF)
      NN = (NF/2)*2
      IF (NN.EQ.NF) NF = NF + 1
C
C IERR = 3 , IF THE MAX. ALLOWED FILTER ORDER IS EXCEEDED.
C
      IF (NF.LE.128) GO TO 40
      IF ((ITYP.NE.2) .AND. (ITYP.NE.3) .AND. (NF.LE.256)) GO TO 40
      IERR = 3
      RETURN
  40  CONTINUE
C
C EXPECTED EFFICIENT ROUNDING ERROR.
C
      IF (ISYM.EQ.1) NF = ((NF+2)/4)*2
      CALL EFFERR(DEL11, NF, IWC, S1)
      CALL EFFERR(DEL21, NF, IWC, S2)
      IF (ISYM.EQ.1) NF = 2*NF - 1
C
C COMPARISON OF ESTIMATED TOTAL ERRORS WITH TOLERANCES.
C
      F1 = DEL11 + S1
      DD1 = F1 - DEL1
      F2 = DEL21 + S2
      DD2 = F2 - DEL2
      IF ((DD1.GT.SMALL).OR.(DD2.GT.SMALL)) GO TO 70
      ITERA = ITERA + 1
      IF (ITERA.EQ.1) NF0 = NF
      IF (ITERA.EQ.1) NF00 = NF0
C
C EQUIVALENT LOWPASS DESIGN FOR REDUCED TOLERANCES.
C
      N1 = (NF+1)/2
      N = N1
      CALL DESIGN(DEL11, DEL21, OMGP1, OMGS1, H0, NF)
C
C COEFFICIENT ROUNDING TO (IWC-1) BITS BEHIND THE DUAL POINT.
C
      Q = 2.**(1-IWC)
      DO 50 I=1,N1
        HR = H0(I)/Q
        RH0(I) = Q*AINT(HR)
  50  CONTINUE
C
C DETERMINATION OF REAL MAX. ERRORS.
C
      CALL MAXERR(RH0, N1, NF, 0., OMGP1, EMAX1)
      CALL MAXERR(RH0, N1, NF, OMGS1, .5, EMAX2)
C
C COMPARISON WITH TOLERANCES.
C
      DD1 = EMAX1 - DEL1
      DD2 = EMAX2 - DEL2
      IF ((DD1.GT.SMALL).OR.(DD2.GT.SMALL)) GO TO 60
      ISOLUT = ISOLUT + 1
C
C CHECKING FOR SUFFICIENT USAGE OF TOLERANCES.
C
      IF ((DD1/DEL1).GE.(PROZ-1.)) DD1 = 0.
      IF ((DD2/DEL2).GE.(PROZ-1.)) DD2 = 0.
      IF ((DD1.EQ.0.) .AND. (DD2.EQ.0.)) GO TO 80
      IF (NF.LE.NF0) GO TO 80
      IF (ISOLUT.EQ.1) NF0 = NF
C
C TOLERANCE CHANGE.
C
  60  CONTINUE
      IF (DD1.LE.0.) DEL11 = DEL11 - CI*DD1
      IF (DD2.LE.0.) DEL21 = DEL21 - CI*DD2
  70  CONTINUE
      IF (DD1.GT.0.) DEL11 = DEL11 - CR*DD1
      IF (DD2.GT.0.) DEL21 = DEL21 - CR*DD2
      IF (ISYM.EQ.1) DEL11 = AMIN1(DEL11,DEL21)
      IF (ISYM.EQ.1) DEL21 = DEL11
C
C CHECK FOR REALIZABILITY.
C
      IF ((DEL11.GT.0.) .AND. (DEL21.GT.0.)) GO TO 20
C
C E N D  O F  I N N E R  D E S I G N  I T E R A T I O N  L O O P .
C
C
C INCREASE OF WORDLENGTH IN CASE OF UNREALIZABILITY.
C
      IWC = IWC + 1
      GO TO 10
C
C E N D  O F  O U T E R  D E S I G N  I T E R A T I O N  L O O P .
C
C
C TRANSFORMATION TO DESIRED FILTER TYPE.
C
  80  CONTINUE
      IF (ITYP.EQ.0) GO TO 90
      CALL BACTRA(ITYP, H0, RH0, NF)
      N1 = (NF+1)/2
  90  CONTINUE
C
C DETERMINATION OF ACTUAL COEFFICIENT WORDLENGTH.
C
      CALL WCACT(RH0, N1, IWC, IWACT)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TYPTRA
C TRANSFORMATION TO EQUIVALENT LOWPASS;SPECIFICATIONS TEST.
C-----------------------------------------------------------------------
C
      SUBROUTINE TYPTRA(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP,
     *    IERR)
C
C     TRANSFORMATION TO EQUIVALENT LOWPASS PROBLEM.
C     IDENTIFICATION OF FILTER TYPES- ITYP = 0..LOWPASS)
C                                     ITYP = 1..HIGHPASS)
C                                     ITYP = 2..BANDPASS)
C                                     ITYP = 3..BANDSTOP)
C                                     ITYP = 4..HILBERT-TRANSFORMER
C
C L O W P A S S
C
      IF ((OMGP2.NE.0.) .OR. (OMGS2.NE.0.)) GO TO 20
      IF (OMGP1.GT.OMGS1) GO TO 10
      ITYP = 0
      GO TO 60
C
C H I G H P A S S
C
  10  CONTINUE
      IF (OMGS1.EQ.0.) GO TO 40
      ITYP = 1
      OMGP1 = .5 - OMGP1
      OMGS1 = .5 - OMGS1
      GO TO 60
  20  CONTINUE
C
C CHECKING FOR INCONSISTENT SPECIFICATIONS.
C SET UP A MACHINE CONSTANT.
C
      EPSMAC = 2.*R1MACH(3)
      DOP = ABS(OMGP1+OMGP2-.5)
      DOS = ABS(OMGS1+OMGS2-.5)
      DO = AMAX1(DOP,DOS)
      IF (DO.GE.EPSMAC) GO TO 50
C
C B A N D P A S S
C
      IF (OMGP1.LT.OMGS1) GO TO 30
      ITYP = 2
      OMGP1 = .5 - 2.*OMGP1
      OMGS1 = .5 - 2.*OMGS1
      GO TO 60
C
C B A N D S T O P
C
  30  CONTINUE
      ITYP = 3
      OMGP1 = 2.*OMGP1
      OMGS1 = 2.*OMGS1
      GO TO 60
C
C H I L B E R T - T R A N S F O R M E R
C
  40  CONTINUE
C
C CHECKING FOR INCONSISTENT SPECIFICATIONS.
C
      IF ((DEL1.NE.DEL2) .AND. (DEL2.NE.0.)) GO TO 50
      ITYP = 4
      OMGS1 = .25 + OMGP1
      OMGP1 = .25 - OMGP1
      DEL1 = .5*DEL1
      DEL2 = DEL1
      GO TO 60
C
C ERROR - INCONSISTENT SPECIFICATIONS.
C
  50  IERR = 2
C
  60  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: BACTRA
C TRANSFORMATION OF THE DESIGNED LOWPASS TO THE DESIRED FILTER TYPE.
C-----------------------------------------------------------------------
C
      SUBROUTINE BACTRA(ITYP, H0, RH0, NF)
C
      DIMENSION H0(256), RH0(128)
C
      N1 = (NF+1)/2
      IF (ITYP-3) 10, 30, 60
  10  CONTINUE
C
C H I G H P A S S  AND B A N D P A S S
C
      N = N1 - 1
      DO 20 I=1,N,2
        J = N1 - I
        H0(J) = -H0(J)
        RH0(J) = -RH0(J)
  20  CONTINUE
      IF (ITYP.EQ.1) GO TO 80
  30  CONTINUE
C
C B A N D P A S S  AND  B A N D S T O P
C
      N11 = N1
      N1 = NF
      NF = 2*NF - 1
      NN = 2*N11
      I = N1
      J = N11
      IF (NN-N1) 40, 40, 50
  40  CONTINUE
C
C LOWPASS FILTER LENGTH EVEN.
C
      H0(I) = 0.
      RH0(I) = 0.
      I = I - 1
  50  CONTINUE
C
C LOWPASS FILTER LENGTH EVEN OR ODD.
C
      H0(I) = H0(J)
      H0(I-1) = 0.
      RH0(I) = RH0(J)
      RH0(I-1) = 0.
      I = I - 2
      J = J - 1
      IF (J-1) 80, 80, 50
  60  CONTINUE
C
C H I L B E R T - T R A N S F O R M E R
C
      I = N1
      H0(I) = 0.
      RH0(I) = 0.
      V = -2.
      I = I - 1
  70  CONTINUE
      V = -V
      H0(I) = V*H0(I)
      RH0(I) = V*RH0(I)
      I = I - 2
      IF (I) 80, 80, 70
  80  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: WCMIN
C MINIMUM COEFFICIENT WORDLENGTH ESTIMATION BY AN EMPIRICAL FORMULA.
C-----------------------------------------------------------------------
C
      SUBROUTINE WCMIN(D1, D2, DOMG, IW)
C
      ALD = ALOG10(D1)
      ALS = ALOG10(D2)
      ALO = ALOG10(DOMG*360.)
      W = 7.22116 - 0.589268*ALD
      W = W + (0.112354*ALD-2.328539)*ALO
      W = W - 3.559051*ALS
      IW = IFIX(W)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: LENGTH
C FILTER LENGTH ESTIMATION BY AN EMPIRICAL FORMULA.
C-----------------------------------------------------------------------
C
      SUBROUTINE LENGTH(D1, D2, DOMG, NF)
C
      DL1 = ALOG10(D1)
      DL2 = ALOG10(D2)
      DQ1 = DL1*DL1
      D = .005309*DQ1 + .07114*DL1 - .4761
      D = D*DL2
      D = D - .00266*DQ1 - .5941*DL1 - .4278
      D = 2. + D/DOMG
      NF = IFIX(D)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: EFFERR
C ITERATIVE CALCULATION OF THE APPROXIMATE STATISTICAL ERROR BOUND S
C FOR A GIVEN TOLERANCE D,FILTER LENGTH NF,COEFF.WORDLENGTH IW.
C-----------------------------------------------------------------------
C
      SUBROUTINE EFFERR(D, NF, IW, S)
C
      QH = 2.**(-IW)
      CALL SSTAT(NF, 0., SQ)
      ALPHA = D/(SQ*QH)
  10  CONTINUE
      CALL SSTAT(NF, ALPHA, SQ)
      ALPHA1 = D/(SQ*QH)
      DALPHA = ALPHA1 - ALPHA
      IF (DALPHA-.05) 30, 30, 20
  20  CONTINUE
      ALPHA = ALPHA1
      GO TO 10
  30  CONTINUE
      S = QH*SQ
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SSTAT
C EVALUATION OF THE APPROXIMATE STATISTICAL ERROR FORMULA.
C-----------------------------------------------------------------------
C
      SUBROUTINE SSTAT(NF, ALPHA, SQ)
C
      PI = 4.*ATAN(1.0)
      PIQ = PI*PI
C
C WEIGHTING FUNCTION M(ALPHA).
C
      AM = 1./(1.+SQRT(ALPHA))
C
C ERROR BOUND.
C
      ANZ = FLOAT(NF-1)
      SQ = (ANZ/PI+.5)*AM
      X = (4.-24.*AM/PIQ)*ANZ
      X = X + 4. - 3.*AM
      X = SQRT(X*AM/3.)
      SQ = SQ + X
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: MAXERR
C SEARCH FOR THE MAX.ABS.ERROR OF A LINEAR-PHASE FIR FILTER.
C-----------------------------------------------------------------------
C
      SUBROUTINE MAXERR(RH0, N1, NF, OMGA, OMGE, EMAX)
C
      DIMENSION RH0(128)
      EMAX = 0.
      OMG = OMGA
      IA = 1
      ISTG = 1
      DOM1 = 1./(4.*FLOAT(NF))
      DOM2 = DOM1/2.
      DOM3 = DOM2/2.
      DFR = .05
      DFRH = .025
      CALL ERROR(RH0, N1, NF, OMG, E)
      BF1 = ABS(E)
      IF (BF1.GT.EMAX) EMAX = BF1
      DELOMG = DOM3
      OMG = OMG + DELOMG
  10  CONTINUE
      BF2 = BF1
  20  CONTINUE
      CALL ERROR(RH0, N1, NF, OMG, E)
      BF1 = ABS(E)
      IF (BF1.GT.EMAX) EMAX = BF1
      IF (IA) 30, 110, 30
  30  CONTINUE
      DELFR = 2.*(BF1-BF2)/(BF1+BF2)
      IF (ISTG) 90, 90, 40
  40  CONTINUE
      IF (DELFR) 80, 80, 50
  50  CONTINUE
      ISTG = 1
      DELOMG = DOM2
      IF (DELFR.LT.DFRH) DELOMG = DOM3
      IF (DELFR.GT.DFR) DELOMG = DOM1
      OMG = OMG + DELOMG
  60  CONTINUE
      IF (OMG-OMGE) 10, 70, 70
  70  CONTINUE
      OMG = OMGE
      IA = 0
      GO TO 10
  80  CONTINUE
      ISTG = -1
      IF (DELOMG.EQ.DOM1) DELOMG = -3.*DOM3
      IF (DELOMG.EQ.DOM2) DELOMG = -DOM3
      IF (DELOMG.EQ.DOM3) DELOMG = DOM1
      OMG = OMG + DELOMG
      IF (DELOMG) 20, 60, 60
  90  CONTINUE
      IF (DELFR) 100, 100, 50
 100  CONTINUE
      DELOMG = DOM1
      OMG = OMG + DELOMG
      GO TO 60
 110  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ERROR
C ERROR E OF A LINEAR-PHASE FIR LOWPASS TRANSFER FUNCTION WITH
C RESPECT TO THE IDEAL LOWPASS.
C-----------------------------------------------------------------------
C
      SUBROUTINE ERROR(RH0, N1, NF, OMG, E)
C
      DIMENSION RH0(128)
C
      PI = 4.*ATAN(1.0)
      PI2 = PI*2.
      NN = 2*N1
      IF (NF-NN) 10, 30, 30
  10  CONTINUE
C
C LENGTH NF ODD.
C
      H = RH0(N1)
      N = N1 - 1
      DO 20 I=1,N
        ANNU = FLOAT(N1-I)
        H = H + 2.*RH0(I)*COS(ANNU*OMG*PI2)
  20  CONTINUE
      GO TO 50
  30  CONTINUE
C
C LENGTH NF EVEN.
C
      H = 0.
      DO 40 I=1,N1
        ANNU = FLOAT(N1-I) + .5
        H = H + RH0(I)*COS(ANNU*OMG*PI2)
  40  CONTINUE
      H = H*2.
  50  CONTINUE
C
C ERROR VALUE E.
C
      E = H
      E1 = H - 1.
      EA = ABS(E)
      EA1 = ABS(E1)
      IF (EA1.LT.EA) E = E1
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: WCACT
C DETERMINATION OF ACTUAL WORDLENGTH.
C-----------------------------------------------------------------------
C
      SUBROUTINE WCACT(RH0, N1, IWC, IWACT)
C
      DIMENSION RH0(128)
      HMAX = ABS(RH0(1))
      DO 10 I=2,N1
        H = ABS(RH0(I))
        HMAX = AMAX1(HMAX,H)
  10  CONTINUE
      L = -2
  20  CONTINUE
      L = L + 1
      IF (HMAX-1.) 30, 40, 40
  30  CONTINUE
      HMAX = HMAX*2.
      GO TO 20
  40  CONTINUE
      IWACT = IWC - L
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESIGN
C AUXILIARY ROUTINE CALLING THE ODD-LENGTH LOWPASS DESIGN ROUTINE
C 'TOMSB2'.
C-----------------------------------------------------------------------
C
      SUBROUTINE DESIGN(DEL11, DEL21, OMGP1, OMGS1, H0, NF)
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
      COMMON /CHEK/ NL, NG, KA, KB
      DIMENSION H0(256)
C
C AUXILIARY VARIABLES FOR 'TOMSB2'.
C
      FP = OMGP1
      FS = OMGS1
      CALL TOMSB2(DEL11, DEL21, FP, FS, H0, NF)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TOMSB2
C     SUBROUTINE ACCORDING TO "PROGRAM TOM" WITH  THE
C     EXCHANGE ALGORITHM FOR DESIGNING LOW PASS FINITE LENGTH IMPULSE
C     RESPONSE DIGITAL FILTERS WITH LINEAR PHASE.
C
C     J.H. MCCLELLAN AND T.W. PARKS---RICE UNIVERSITY---OCTOBER 26, 1971
C
C     THE FILTER LENGTH (NF), THE PASSBAND CUTOFF FREQUENCY (FP),
C     THE STOPBAND FREQUENCY (FS), AND THE WEIGHTING FACTOR (AA) ARE
C     INPUTS TO THE ALGORITHM.  THEN THE STOPBAND DEVIATION (DS) IS
C     MINIMIZED AND THUS THE PASSBAND DEVIATION (DP=AA*DS) IS ALSO
C     MINIMIZED. THE PROGRAM SEKS TO DETERMINE THE BEST LOCATION OF
C     THE PEAKS OF THE ERROR CURVE IN ORDER TO MINIMIZE THE DEVIATIONS.
C
C     REFERENCE--PARKS,T.W. AND J.H. MCCLELLAN,*CHEBYCHEV APPROXIMATION
C     OF NON-RECURSIVE DIGITAL FILTERS WITH LINEAR PHASE*, IEEE TRANS.
C     ON CIRCUIT THEORY, VOL. CT-19, MARCH, 1972.
C
C     SPECIAL PARAMETERS FOR CHANGE OF FS, IF STOPBAND RIPPLE IS
C     SPECIFIED
C     D1D   =  STOPBAND RIPPLE D2
C     STEP  =  INITIAL STEPSIZE FOR CHANGING FS& STEP DECREASES WITH
C              EXP.& IF SOLUTION  FAR FROM OPTIMUM TRY 0.02
C     XAC   =  ACCURACY OF D2
C              1000. MEANS ACCURACY TO 0.1 PERCENT
C              100.  MEANS ACCURACY TO 1 PERCENT ETC.
C     IF XAC = D1D  NO ACTION IS TAKEN AND FS REMAINS UNCHANGED
C
C        DIMENSION OF X, AD, Y, IFR, FOPT = N*/2+3 WHERE N* IS
C        MAXIMUM LENGTH OF N
C        DIMENSION OF AF = 16(N*/2+3)+1
C-----------------------------------------------------------------------
C
      SUBROUTINE TOMSB2(DEL1, DEL2, FP, FS, H0, NF)
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
      COMMON /CHEK/ NL, NG, KA, KB
      DIMENSION H0(256)
C
      IWCHN = I1MACH(2)
      DO 5 I=1,131
        X(I) = 0.
   5  CONTINUE
C
C INITIALIZATION OF THE GRID OF 16N POINTS.  THE CONSTANT LGRID
C CAN BE CHANGED IF A DIFFERENT GRID DENSITY IS DESIRED.
C
      LGRID = 16
      XAC = DEL2
      AA = DEL1/DEL2
      D1D = DEL2
      STEP = 0.02
      ICH = 1
      KITER = 0
      D1C = D1D/XAC
  10  CONTINUE
      PI2 = 8.0*ATAN(1.0)
      KCOUNT = 0
      TW = FS - FP
      NGRID = LGRID*(N-1)
      NGRID = NGRID + LGRID/2
      KGRID = 2
      DELF = 0.5/FLOAT(NGRID)
      AF(1) = 0.0
      DO 20 J=1,NGRID
        AF(J+1) = FLOAT(J)*DELF
  20  CONTINUE
      KFP = FP/DELF - 1.
      KFS = FS/DELF
  30  CONTINUE
      IF (FP.LT.AF(KFP+1)) GO TO 40
      KFP = KFP + 1
      GO TO 30
  40  CONTINUE
      IF (FS.LE.AF(KFS)) GO TO 50
      KFS = KFS + 1
      GO TO 40
  50  CONTINUE
      AF(KFP) = FP
      AF(KFS) = FS
      NGRID = NGRID + 1
      DO 60 J=1,NGRID
        AF(J) = COS(PI2*AF(J))
  60  CONTINUE
C
C INITIAL GUESS FOR THE N+1 OPTIMAL FREQUENCIES
C
      KP = (KFP-1)/2
      KP = 2*KP + 1
      KS = KFS/2
      KS = 2*KS + 1
      XPP = KP - 1
      XSS = NGRID - KS
      XAA = XPP + XSS
      PT = XAA/FLOAT(N-1)
      IF (PT.GT.3.) GO TO 70
      WRITE (IWCHN,9999)
9999  FORMAT (20H ERROR IN INPUT DATA)
      GO TO 560
  70  KR = XPP/PT
      IF (KR.EQ.0) KR = KR + 1
      PT = XPP/FLOAT(KR)
      IFR(1) = 1
      DO 80 J=1,KR
        XTW = FLOAT(J)*PT
        NTW = IFIX(XTW) + 1
        NTW = NTW/2
        IFR(J+1) = 2*NTW + 1
  80  CONTINUE
      KR = KR + 1
      IFR(KR) = KP
      IFR(KR+1) = KS
      NTW = N - KR
      IF (NTW.GE.0) GO TO 90
      WRITE (IWCHN,9999)
      GO TO 560
  90  IF (NTW.EQ.0) GO TO 110
      PT = XSS/FLOAT(NTW)
      DO 100 J=1,NTW
        XTW = FLOAT(J)*PT + FLOAT(KS) - 1.0
        NB = XTW
        NB = NB/2
        K111 = KR + 1 + J
        IFR(K111) = 2*NB - 1
 100  CONTINUE
 110  CONTINUE
C
C MAIN CALCULATION ROUTINE
C CALCULATE RHO ON THE OPTIMAL SET
C THE FIRST 4 PASSES ARE MADE ON A GRID OF 8N POINTS THEN THE
C FINER GRID OF 16N POINTS IS USED.
C
 120  CONTINUE
      N111 = N + 1
      DO 130 J=1,N111
        J111 = IFR(J)
        X(J) = AF(J111)
 130  CONTINUE
      M = (N-1)/15 + 1
      M = (N-1)/20 + 1
      N111 = N + 1
      DO 140 J=1,N111
        AD(J) = D(J,N+1,M)
 140  CONTINUE
      R1 = 0.0
      DO 150 J=1,KR
        R1 = R1 + AD(J)
 150  CONTINUE
      R2 = 0.0
      K = 1
      DO 160 J=1,KR
        R2 = R2 + AD(J)*AA*FLOAT(K)
        K = -K
 160  CONTINUE
      K111 = KR + 1
      N111 = N + 1
      DO 170 J=K111,N111
        R2 = R2 + FLOAT(K)*AD(J)
        K = -K
 170  CONTINUE
      RHO = R1/R2
      NU = 1
      IF (RHO.GT.0.0) NU = -1
      RHO = ABS(RHO)
      K = NU
      XRAY = X(KR+1)
      DO 180 J=1,KR
        AD(J) = AD(J)*(XRAY-X(J))
        Y(J) = 1.0 + AA*FLOAT(K)*RHO
        K = -K
 180  CONTINUE
      K = -K
      K111 = KR + 1
      DO 190 J=K111,N
        X(J) = X(J+1)
        AD(J) = AD(J+1)*(XRAY-X(J))
        Y(J) = FLOAT(K)*RHO
        K = -K
 190  CONTINUE
      Y(N+1) = FLOAT(K)*RHO
      NL = 0
      NG = 0
      TEST = GEE(IFR(KR+1))
      RO = ABS(TEST-RHO)/RHO
      IF (RO.LT.0.1 .OR. KCOUNT.LT.4) GO TO 200
      WRITE (IWCHN,9998)
9998  FORMAT (20H INTERPOLATION ERROR)
      GO TO 560
 200  CONTINUE
C
C SEARCH FOR THE NEW EXTREMAL FREQUENCIES
C
      NTOT = 0
      NPT = IFR(1)
      NG = IFR(2)
      KA = 1
      KB = 1
      KSIGN = NU
      XT = Y(1)
      NFIR = IFR(1)
      NOLD = IFR(N+1)
      IF (KR.EQ.1) GO TO 330
      IF (NPT-1) 210, 260, 210
 210  ZT = GEE(NPT-KGRID)
      IF (FLOAT(KSIGN)*(ZT-XT)) 260, 260, 220
 220  KPT = NPT - KGRID
 230  KPT = KPT - KGRID
      XT = ZT
      IF (KPT-1) 250, 240, 240
 240  ZT = GEE(KPT)
      IF (FLOAT(KSIGN)*(ZT-XT)) 250, 250, 230
 250  IFR(1) = KPT + KGRID
      NTOT = NTOT + 1
      KSIGN = -KSIGN
      GO TO 310
 260  ZT = GEE(NPT+KGRID)
      IF (FLOAT(KSIGN)*(ZT-XT)) 300, 300, 270
 270  KPT = NPT
 280  KPT = KPT + KGRID
      XT = ZT
      ZT = GEE(KPT+KGRID)
      IF (FLOAT(KSIGN)*(ZT-XT)) 290, 290, 280
 290  IFR(1) = KPT
      NTOT = NTOT + 1
 300  KSIGN = -KSIGN
 310  XFIR = ABS(1.0-XT)/AA
      KA = 2
 320  IF (KA.GT.(N+1)) GO TO 420
      NL = NPT
      NG = NGRID + 1
      IF (KA.LE.N) NG = IFR(KA+1)
      NPT = IFR(KA)
      XT = Y(KA)*FLOAT(KB)
      IF (KA.EQ.KR) GO TO 330
      IF (KA.EQ.(KR+1)) GO TO 340
      GO TO 350
 330  IF ((KFP-NPT).LT.KGRID) GO TO 410
      NPT = KFP
      NTOT = 1
      GO TO 410
 340  KB = -1
      IF ((NPT-KFS).LT.KGRID) GO TO 410
      NPT = KFS
      NTOT = 1
      GO TO 410
 350  KPT = NPT
      ZT = GEE(KPT-KGRID)
      IF (FLOAT(KSIGN)*(ZT-XT)) 380, 380, 360
 360  XT = ZT
      KPT = KPT - KGRID
      ZT = GEE(KPT-KGRID)
      IF (FLOAT(KSIGN)*(ZT-XT)) 370, 370, 360
 370  IFR(KA) = KPT
      KSIGN = -KSIGN
      NTOT = NTOT + 1
      KA = KA + 1
      GO TO 320
 380  IF ((KPT+KGRID).GT.NGRID) GO TO 410
      ZT = GEE(KPT+KGRID)
      IF (FLOAT(KSIGN)*(ZT-XT)) 410, 410, 390
 390  XT = ZT
      KPT = KPT + KGRID
      IF ((KGRID+KPT).GT.NGRID) GO TO 400
      ZT = GEE(KPT+KGRID)
      IF (FLOAT(KSIGN)*(ZT-XT)) 400, 400, 390
 400  IFR(KA) = KPT
      KSIGN = -KSIGN
      NTOT = NTOT + 1
      KA = KA + 1
      GO TO 320
 410  IFR(KA) = NPT
      KSIGN = -KSIGN
      KA = KA + 1
      GO TO 320
 420  CONTINUE
      XOLD = ABS(XT)*AA
      NL = 0
      NG = 0
      IF (IFR(1).EQ.1) GO TO 450
      IF (IFR(N+1).EQ.NGRID) GO TO 430
      WRITE (IWCHN,9997)
9997  FORMAT (38H ERROR-NEITHER ENDPOINT IS AN EXTREMUM)
      GO TO 560
 430  IF (NFIR.EQ.1) GO TO 470
      TEST = GEE(1)
      TEST = -FLOAT(NU)*(TEST-1.0) - XOLD
      IF (TEST.LE.0.0) GO TO 470
      NTOT = 1
      DO 440 J=1,N
        N111 = N + 2 - J
        M111 = N + 1 - J
        IFR(N111) = IFR(M111)
 440  CONTINUE
      IFR(1) = 1
      KR = KR + 1
      GO TO 470
 450  IF (IFR(N+1).EQ.NGRID) GO TO 470
      IF (NOLD.EQ.NGRID) GO TO 470
      TEST = GEE(NGRID)
      TEST = FLOAT(KSIGN)*TEST - XFIR
      IF (TEST.LE.0.0) GO TO 470
      NTOT = 1
      DO 460 J=1,N
        IFR(J) = IFR(J+1)
 460  CONTINUE
      IFR(N+1) = NGRID
      KR = KR - 1
 470  CONTINUE
      IF (KR.GT.0) GO TO 480
      WRITE (IWCHN,9996)
9996  FORMAT (6H ERROR)
      GO TO 560
 480  CONTINUE
      KCOUNT = KCOUNT + 1
      IF (KCOUNT.GT.50) GO TO 500
      IF (KCOUNT.EQ.4 .AND. KGRID.EQ.2) GO TO 490
      IF (NTOT.NE.0) GO TO 120
      IF (KGRID.EQ.2) GO TO 490
      GO TO 500
 490  KGRID = 1
      GO TO 120
 500  CONTINUE
      IF (ABS(RHO-D1D).LT.D1C) GO TO 520
      KITER = KITER + 1
      IF (KITER.GT.1000) GO TO 520
      IF (RHO.GT.D1D) GO TO 510
      IF (ICH.EQ.1) STEP = STEP/2.
      FS = FS - STEP
      ICH = -1
      GO TO 10
 510  IF (ICH.EQ.(-1)) STEP = STEP/2.
      FS = FS + STEP
      ICH = 1
      GO TO 10
 520  CONTINUE
      KRP = N - KR
      R1 = AA*RHO
      J = N
 530  CONTINUE
      IF (X(J)) 540, 550, 540
 540  J = J + 1
      GO TO 530
 550  X(J) = -1.
      CALL PREF2(H0, NF)
      GO TO 570
 560  CONTINUE
      STOP
 570  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: D
C THE SUBROUTINE D(,,,) CALCULATES THE LAGRANGIAN INTERPOLATION
C COEFFICIENTS FOR USE IN THE BARYCENTRIC INTERPOLATION FORMULA
C-----------------------------------------------------------------------
C
      FUNCTION D(K, N, M)
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      D = 1.0
      Q = X(K)
      DO 30 L=1,M
        DO 20 J=L,N,M
          IF (J-K) 10, 20, 10
  10      D = 2.0*D*(X(J)-Q)
  20    CONTINUE
  30  CONTINUE
      D = 1.0/D
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: GEE
C THE SUBROUTINE GEE EVALUATES THE FREQUENCY RESPONSE USING THE
C INTERPOLATION COEFFICIENTS WHICH ARE CALCULATED IN THE ROUTINE D.
C-----------------------------------------------------------------------
C
      FUNCTION GEE(K)
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
      COMMON /CHEK/ NL, NG, KA, KB
      IF (K-NL) 10, 40, 10
  10  IF (K-NG) 20, 40, 20
  20  P = 0.0
      XF = AF(K)
      D = 0.0
      DO 30 J=1,N
        C = AD(J)/(XF-X(J))
        D = D + C
        P = P + C*Y(J)
  30  CONTINUE
      GEE = P/D
      RETURN
  40  K111 = KA + KB
      GEE = Y(K111)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: PREF2
C THE SUBROUTINE PREF SAMPLES THE FREQUENCY RESPONSE AT 2**M POINTS
C SO THAT THE FFT (R.C. SINGLETON) CAN BE CALLED TO INVERSE TRANSFORM
C THIS DATA TO OBTAIN THE IMPULSE RESPONSE.
C-----------------------------------------------------------------------
C
      SUBROUTINE PREF2(A, NF)
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
      COMPLEX Z
      DIMENSION A(256), B(256)
      DIMENSION Z(256)
      EQUIVALENCE (B(1),AF(257)), (Z(1),AF(520))
C
C SET UP A MACHINE CONSTANT.
C
      FSH = R1MACH(1)
C
      PI2 = 8.0*ATAN(1.0)
      NX = 2
      MX = 1
  10  IF (NX.GE.NF) GO TO 20
      NX = 2*NX
      MX = MX + 1
      GO TO 10
  20  NN = NX
      MM = MX
      NX = NX/2
      XNN = NN
      DF = 1.0/XNN
      IF (IFR(1).EQ.1) GO TO 30
      P = GOO(1.0)
      GO TO 40
  30  P = Y(1)
  40  A(1) = P
      IF (IFR(N+1).EQ.NGRID) GO TO 50
      P = GOO(-1.0)
      GO TO 60
  50  P = Y(N)
  60  A(NX+1) = P
      L = 1
      N111 = NX - 1
      DO 110 J=1,N111
        AT = DF*FLOAT(J)
        AT = COS(PI2*AT)
  70    AS = X(L)
        IF (AT.GT.AS) GO TO 90
        IF ((AS-AT).LT.FSH) GO TO 80
        L = L + 1
        GO TO 70
  80    A(J+1) = Y(L)
        GO TO 100
  90    IF ((AT-AS).LT.FSH) GO TO 80
        A(J+1) = GOO(AT)
 100    L111 = NN + 1 - J
        A(L111) = A(J+1)
        B(J+1) = 0.0
        IF (L.GT.1) L = L - 1
        M111 = NN + 1 - J
        B(M111) = 0.
 110  CONTINUE
      B(1) = 0.0
      NX1 = NX + 1
      B(NX1) = 0.0
      DO 120 J=1,NN
        Z(J) = CMPLX(A(J),B(J))
 120  CONTINUE
      CALL BITREV(Z, MM)
      CALL FFTFT(Z, MM)
      DO 130 J=1,NX
        B(J) = REAL(Z(J))
 130  CONTINUE
      DO 140 J=1,N
        JJ = N + 1 - J
        A(JJ) = B(J)
 140  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: GOO
C THE SUBROUTINE GOO IS THE SAME AS GEE EXCEPT THAT IT IS CALLED BY
C A FREQUENCY VALUE RATHER THAN BY BY A GRID INDEX VALUE.
C-----------------------------------------------------------------------
C
      FUNCTION GOO(F)
C
      COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
      COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
      XF = F
      P = 0.0
      D = 0.0
      DO 20 J=1,N
        IF (XF.NE.X(J)) GO TO 10
        GOO = Y(J)
        GO TO 30
  10    C = AD(J)/(XF-X(J))
        D = D + C
        P = P + C*Y(J)
  20  CONTINUE
      GOO = P/D
  30  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTFT
C FREQUENCY-TO-TIME FFT ROUTINE; INPLACE ALGORITHM; COMPLEX DATA IN
C ARRAY A; INPUT DATA ASSUMED TO BE IN BITREVERSED ORDER.
C-----------------------------------------------------------------------
C
      SUBROUTINE FFTFT(A, M)
C
      COMPLEX A, U, W, T
      DIMENSION A(1)
      PI = 4.0*ATAN(1.)
      N = 2**M
      DO 30 L=1,M
        LE = 2**L
        LE1 = LE/2
        U = (1.0,0.0)
        W = CMPLX(COS(PI/FLOAT(LE1)),SIN(PI/FLOAT(LE1)))
        DO 20 J=1,LE1
          DO 10 I=J,N,LE
            IP = I + LE1
            T = A(IP)*U
            A(IP) = A(I) - T
            A(I) = A(I) + T
  10      CONTINUE
          U = U*W
  20    CONTINUE
  30  CONTINUE
      DO 40 K=1,N
        A(K) = A(K)/FLOAT(N)
  40  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: BITREV
C REORDERING OF 2**M DATA IN COMPLEX ARRAY A IN BITREVERSED ORDER.
C-----------------------------------------------------------------------
C
      SUBROUTINE BITREV(A, M)
C
      DIMENSION A(1)
      COMPLEX A, T
C
      N = 2**M
      NV2 = N/2
      NM1 = N - 1
      J = 1
      DO 40 I=1,NM1
        IF (I.GE.J) GO TO 10
        T = A(J)
        A(J) = A(I)
        A(I) = T
  10    K = NV2
  20    IF (K.GE.J) GO TO 30
        J = J - K
        K = K/2
        GO TO 20
  30    J = J + K
  40  CONTINUE
      RETURN
      END
Return to Main Software Page