C
C-----------------------------------------------------------------------
C MAIN PROGRAM: FIR LINEAR PHASE FILTER DESIGN PROGRAM
C
C AUTHORS: JAMES H. MCCLELLAN
C          DEPARTMENT OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE
C          MASSACHUSETTS INSTITUTE OF TECHNOLOGY
C          CAMBRIDGE, MASS. 02139
C
C          THOMAS W. PARKS
C          DEPARTMENT OF ELECTRICAL ENGINEERING
C          RICE UNIVERSITY
C          HOUSTON, TEXAS 77001
C
C          LAWRENCE R. RABINER
C          BELL LABORATORIES
C          MURRAY HILL, NEW JERSEY 07974
C
C INPUT:
C  NFILT-- FILTER LENGTH
C  JTYPE-- TYPE OF FILTER
C          1 = MULTIPLE PASSBAND/STOPBAND FILTER
C          2 = DIFFERENTIATOR
C          3 = HILBERT TRANSFORM FILTER
C  NBANDS-- NUMBER OF BANDS
C  LGRID-- GRID DENSITY, WILL BE SET TO 16 UNLESS
C          SPECIFIED OTHERWISE BY A POSITIVE CONSTANT.
C
C  EDGE(2*NBANDS)-- BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND
C                   WITH A MAXIMUM OF 10 BANDS.
C
C  FX(NBANDS)-- DESIRED FUNCTION ARRAY (OR DESIRED SLOPE IF A
C               DIFFERENTIATOR) FOR EACH BAND.
C
C  WTX(NBANDS)-- WEIGHT FUNCTION ARRAY IN EACH BAND.  FOR A
C                DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY
C                PROPORTIONAL TO F.
C
C  SAMPLE INPUT DATA SETUP:
C       32,1,3,0
C       0.0,0.1,0.2,0.35
C	0.425,0.5
C       0.0,1.0,0.0
C       10.0,1.0,10.0
C     THIS DATA SPECIFIES A LENGTH 32 BANDPASS FILTER WITH
C     STOPBANDS 0 TO 0.1 AND 0.425 TO 0.5, AND PASSBAND FROM
C     0.2 TO 0.35 WITH WEIGHTING OF 10 IN THE STOPBANDS AND 1
C     IN THE PASSBAND.  THE GRID DENSITY DEFAULTS TO 16.
C     THIS IS THE FILTER IN FIGURE 10.
C
C     THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 FULLBAND
C     DIFFERENTIATOR WITH SLOPE 1 AND WEIGHTING OF 1/F.
C     THE GRID DENSITY WILL BE SET TO 20.
C       32,2,1,20
C       0,0.5
C       1.0
C       1.0
C
C-----------------------------------------------------------------------
C

      SUBROUTINE REMEZF ( NFILT, JTYPE, NBANDS, LGRID, EDGE, FX, WTX, H)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION EDGE, FX, WTX, H

      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
      COMMON /OOPS/NITER
      DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
      DIMENSION H(66)
      DIMENSION DES(1045),GRID(1045),WT(1045)
      DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
      DOUBLE PRECISION PI2,PI
      DOUBLE PRECISION AD,DEV,X,Y
      DOUBLE PRECISION GEED,DD
      INTEGER BD1,BD2,BD3,BD4
      DATA BD1,BD2,BD3,BD4/1HB,1HA,1HN,1HD/
      PI=4.0*DATAN(1.0D0)
      PI2=2.0D00*PI
C
C  THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128, BUT
C  THIS UPPER LIMIT CAN BE CHANGED BY REDIMENSIONING THE
C  ARRAYS IEXT, AD, ALPHA, X, Y, H TO BE NFMAX/2 + 2.
C  THE ARRAYS DES, GRID, AND WT MUST DIMENSIONED
C  16(NFMAX/2 + 2).
C
      NFMAX=128
  100 CONTINUE
C *      JTYPE=0
C
C  PROGRAM INPUT SECTION
C
C *      READ(*,110) NFILT,JTYPE,NBANDS,LGRID
      IF(NFILT.EQ.0)STOP
  110 FORMAT(4I5)
      IF(NFILT.LE.NFMAX.AND.NFILT.GE.3) GO TO 115
      CALL ERR51
      STOP
  115 IF(NBANDS.LE.0) NBANDS=1
C
C  GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED
C  OTHERWISE
C
      IF(LGRID.LE.0) LGRID=16
      JB=2*NBANDS
C *      READ(*,120) (EDGE(J),J=1,JB)
C *  120 FORMAT(4F15.9)
C *      READ(*,120) (FX(J),J=1,NBANDS)
C *      READ(*,120) (WTX(J),J=1,NBANDS)
      IF(JTYPE.GT.0.AND.JTYPE.LE.3) GO TO 125
      CALL ERR51
      STOP
  125 NEG=1
      IF(JTYPE.EQ.1) NEG=0
      NODD=NFILT/2
      NODD=NFILT-2*NODD
      NFCNS=NFILT/2
      IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS=NFCNS+1
C
C  SET UP THE DENSE GRID.  THE NUMBER OF POINTS IN THE GRID
C  IS (FILTER LENGTH + 1)*GRID DENSITY/2
C
      GRID(1)=EDGE(1)
      DELF=LGRID*NFCNS
      DELF=0.5/DELF
      IF(NEG.EQ.0) GO TO 135
      IF(EDGE(1).LT.DELF) GRID(1)=DELF
  135 CONTINUE
      J=1
      L=1
      LBAND=1
  140 FUP=EDGE(L+1)
  145 TEMP=GRID(J)
C
C  CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
C  FUNCTION ON THE GRID
C
      DES(J)=EFF(TEMP,FX,WTX,LBAND,JTYPE)
      WT(J)=WATE(TEMP,FX,WTX,LBAND,JTYPE)
      J=J+1
      GRID(J)=TEMP+DELF
      IF(GRID(J).GT.FUP) GO TO 150
      GO TO 145
  150 GRID(J-1)=FUP
      DES(J-1)=EFF(FUP,FX,WTX,LBAND,JTYPE)
      WT(J-1)=WATE(FUP,FX,WTX,LBAND,JTYPE)
      LBAND=LBAND+1
      L=L+2
      IF(LBAND.GT.NBANDS) GO TO 160
      GRID(J)=EDGE(L)
      GO TO 140
  160 NGRID=J-1
      IF(NEG.NE.NODD) GO TO 165
      IF(GRID(NGRID).GT.(0.5-DELF)) NGRID=NGRID-1
  165 CONTINUE
C
C  SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
C  TO THE ORIGINAL PROBLEM
C
      IF(NEG) 170,170,180
  170 IF(NODD.EQ.1) GO TO 200
      DO 175 J=1,NGRID
      CHANGE=DCOS(PI*GRID(J))
      DES(J)=DES(J)/CHANGE
  175 WT(J)=WT(J)*CHANGE
      GO TO 200
  180 IF(NODD.EQ.1) GO TO 190
      DO 185 J=1,NGRID
      CHANGE=DSIN(PI*GRID(J))
      DES(J)=DES(J)/CHANGE
  185 WT(J)=WT(J)*CHANGE
      GO TO 200
  190 DO 195 J=1,NGRID
      CHANGE=DSIN(PI2*GRID(J))
      DES(J)=DES(J)/CHANGE
  195 WT(J)=WT(J)*CHANGE
C
C  INITIAL GUESS FOR THE EXTREMAL FREQUENCIES--EQUALLY
C  SPACED ALONG THE GRID
C
  200 TEMP=FLOAT(NGRID-1)/FLOAT(NFCNS)
      DO 210 J=1,NFCNS
      XT=J-1
  210 IEXT(J)=XT*TEMP+1.0
      IEXT(NFCNS+1)=NGRID
      NM1=NFCNS-1
      NZ=NFCNS+1
C
C  CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION
C  PROBLEM
C
      CALL REMEZ
C
C  CALCULATE THE IMPULSE RESPONSE.
C
      IF(NEG) 300,300,320
  300 IF(NODD.EQ.0) GO TO 310
      DO 305 J=1,NM1
      NZMJ=NZ-J
  305 H(J)=0.5*ALPHA(NZMJ)
      H(NFCNS)=ALPHA(1)
      GO TO 350
  310 H(1)=0.25*ALPHA(NFCNS)
      DO 315 J=2,NM1
      NZMJ=NZ-J
      NF2J=NFCNS+2-J
  315 H(J)=0.25*(ALPHA(NZMJ)+ALPHA(NF2J))
      H(NFCNS)=0.5*ALPHA(1)+0.25*ALPHA(2)
      GO TO 350
  320 IF(NODD.EQ.0) GO TO 330
      H(1)=0.25*ALPHA(NFCNS)
      H(2)=0.25*ALPHA(NM1)
      DO 325 J=3,NM1
      NZMJ=NZ-J
      NF3J=NFCNS+3-J
  325 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF3J))
      H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(3)
      H(NZ)=0.0
      GO TO 350
  330 H(1)=0.25*ALPHA(NFCNS)
      DO 335 J=2,NM1
      NZMJ=NZ-J
      NF2J=NFCNS+2-J
  335 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF2J))
      H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(2)
C
C  PROGRAM OUTPUT SECTION.
C

  350 GO TO 470

C *  350 WRITE(*,360)
C *  360 FORMAT(1H1, 70(1H*)//15X,29HFINITE IMPULSE RESPONSE (FIR)/
C *     113X,34HLINEAR PHASE DIGITAL FILTER DESIGN/
C *     217X,24HREMEZ EXCHANGE ALGORITHM/)
C *      IF(JTYPE.EQ.1) WRITE(*,365)
C *  365 FORMAT(22X,15HBANDPASS FILTER/)
C *      IF(JTYPE.EQ.2) WRITE(*,370)
C *  370 FORMAT(22X,14HDIFFERENTIATOR/)
C *      IF(JTYPE.EQ.3) WRITE(*,375)
C *  375 FORMAT(20X,19HHILBERT TRANSFORMER/)
C *      WRITE(*,378) NFILT
C *  378 FORMAT(20X,16HFILTER LENGTH = ,I3/)
C *      WRITE(*,380)
C *  380 FORMAT(15X,28H***** IMPULSE RESPONSE *****)
C *      DO 381 J=1,NFCNS
C *      K=NFILT+1-J
C *      IF(NEG.EQ.0) WRITE(*,382) J,H(J),K
C *      IF(NEG.EQ.1) WRITE(*,383) J,H(J),K
C *  381 CONTINUE
C *  382 FORMAT(13X,2HH(,I2,4H) = ,E15.8,5H = H(,I3,1H))
C *  383 FORMAT(13X,2HH(,I2,4H) = ,E15.8,6H = -H(,I3,1H))
C *      IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(*,384) NZ
C *  384 FORMAT(13X,2HH(,I2,8H) =  0.0)
C *      DO 450 K=1,NBANDS,4
C *      KUP=K+3
C *      IF(KUP.GT.NBANDS) KUP=NBANDS
C *      WRITE(*,385) (BD1,BD2,BD3,BD4,J,J=K,KUP)
C *  385 FORMAT(/24X,4(4A1,I3,7X))
C *      WRITE(*,390) (EDGE(2*J-1),J=K,KUP)
C *  390 FORMAT(2X,15HLOWER BAND EDGE,5F14.7)
C *      WRITE(*,395) (EDGE(2*J),J=K,KUP)
C *  395 FORMAT(2X,15HUPPER BAND EDGE,5F14.7)
C *      IF(JTYPE.NE.2) WRITE(*,400) (FX(J),J=K,KUP)
C *  400 FORMAT(2X,13HDESIRED VALUE,2X,5F14.7)
C *      IF(JTYPE.EQ.2) WRITE(*,405) (FX(J),J=K,KUP)
C *  405 FORMAT(2X,13HDESIRED SLOPE,2X,5F14.7)
C *      WRITE(*,410) (WTX(J),J=K,KUP)
C *  410 FORMAT(2X,9HWEIGHTING,6X,5F14.7)
C *      DO 420 J=K,KUP
C *  420 DEVIAT(J)=DEV/WTX(J)
C *      WRITE(*,425) (DEVIAT(J),J=K,KUP)
C *  425 FORMAT(2X,9HDEVIATION,6X,5F14.7)
C *      IF(JTYPE.NE.1) GO TO 450
C *      DO 430 J=K,KUP
C *  430 DEVIAT(J)=20.0*DLOG10(DEVIAT(J)+FX(J))
C *      WRITE(*,435) (DEVIAT(J),J=K,KUP)
C *  435 FORMAT(2X,15HDEVIATION IN DB,5F14.7)
C *  450 CONTINUE
C *      DO 452 J=1,NZ
C *      IX=IEXT(J)
C *  452 GRID(J)=GRID(IX)
C *      WRITE(*,455) (GRID(J),J=1,NZ)
C *  455 FORMAT(/2X,47HEXTREMAL FREQUENCIES--MAXIMA OF THE ERROR CURVE/
C *     1 (2X,5F12.7))
C *      WRITE(*,460)
C *  460 FORMAT(/1X,70(1H*)/1H1)
C *        GO TO 100

  470 RETURN

      END
C
C-----------------------------------------------------------------------
C FUNCTION: EFF
C   FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE
C   AS A FUNCTION OF FREQUENCY.
C   AN ARBITRARY FUNCTION OF FREQUENCY CAN BE
C   APPROXIMATED IF THE USER REPLACES THIS FUNCTION
C   WITH THE APPROPRIATE CODE TO EVALUATE THE IDEAL
C   MAGNITUDE.  NOTE THAT THE PARAMETER FREQ IS THE
C   VALUE OF NORMALIZED FREQUENCY NEEDED FOR EVALUATION.
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION FUNCTION EFF(FREQ,FX,WTX,LBAND,JTYPE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION FX,WTX
      DIMENSION FX(5),WTX(5)
      IF(JTYPE.EQ.2) GO TO 1
      EFF=FX(LBAND)
      RETURN
    1 EFF=FX(LBAND)*FREQ
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: WATE
C   FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION
C   OF FREQUENCY.  SIMILAR TO THE FUNCTION EFF, THIS FUNCTION CAN
C   BE REPLACED BY A USER-WRITTEN ROUTINE TO CALCULATE ANY
C   DESIRED WEIGHTING FUNCTION.
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION FUNCTION WATE(FREQ,FX,WTX,LBAND,JTYPE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION FX,WTX
      DIMENSION FX(5),WTX(5)
      IF(JTYPE.EQ.2) GO TO 1
      WATE=WTX(LBAND)
      RETURN
    1 IF(FX(LBAND).LT.0.0001) GO TO 2
      WATE=WTX(LBAND)/FREQ
      RETURN
    2 WATE=WTX(LBAND)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ERR51
C   THIS ROUTINE WRITES AN ERROR MESSAGE IF AN
C   ERROR HAS BEEN DETECTED IN THE INPUT DATA.
C-----------------------------------------------------------------------
C
      SUBROUTINE ERR51
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /OOPS/NITER
C *      WRITE(*,1)
C *    1 FORMAT(44H ************ ERROR IN INPUT DATA **********)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: REMEZ
C   THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
C   FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS
C   FUNCTION WITH A SUM OF COSINES.  INPUTS TO THE SUBROUTINE
C   ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
C   DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
C   GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE
C   EXTREMAL FREQUENCIES.  THE PROGRAM MINIMIZES THE CHEBYSHEV
C   ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
C   FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
C   THE COEFFICIENTS OF THE BEST APPROXIMATION.
C-----------------------------------------------------------------------
C
      SUBROUTINE REMEZ
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
      COMMON /OOPS/NITER
      DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
      DIMENSION DES(1045),GRID(1045),WT(1045)
      DIMENSION A(66),P(65),Q(65)
      DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
      DOUBLE PRECISION DK,DAK
      DOUBLE PRECISION AD,DEV,X,Y
      DOUBLE PRECISION GEED,DD
C
C  THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
C
      ITRMAX=25
      DEVL=-1.0
      NZ=NFCNS+1
      NZZ=NFCNS+2
      NITER=0
  100 CONTINUE
      IEXT(NZZ)=NGRID+1
      NITER=NITER+1
      IF(NITER.GT.ITRMAX) GO TO 400
      DO 110 J=1,NZ
      JXT=IEXT(J)
      DTEMP=GRID(JXT)
      DTEMP=DCOS(DTEMP*PI2)
  110 X(J)=DTEMP
      JET=(NFCNS-1)/15+1
      DO 120 J=1,NZ
  120 AD(J)=DD(J,NZ,JET)
      DNUM=0.0
      DDEN=0.0
      K=1
      DO 130 J=1,NZ
      L=IEXT(J)
      DTEMP=AD(J)*DES(L)
      DNUM=DNUM+DTEMP
      DTEMP=FLOAT(K)*AD(J)/WT(L)
      DDEN=DDEN+DTEMP
  130 K=-K
      DEV=DNUM/DDEN
C *      WRITE(*,131) DEV
C *  131 FORMAT(1X,12HDEVIATION = ,F12.9)
      NU=1
      IF(DEV.GT.0.0) NU=-1
      DEV=-FLOAT(NU)*DEV
      K=NU
      DO 140 J=1,NZ
      L=IEXT(J)
      DTEMP=FLOAT(K)*DEV/WT(L)
      Y(J)=DES(L)+DTEMP
  140 K=-K
      IF(DEV.GT.DEVL) GO TO 150
      CALL OUCH
      GO TO 400
  150 DEVL=DEV
      JCHNGE=0
      K1=IEXT(1)
      KNZ=IEXT(NZ)
      KLOW=0
      NUT=-NU
      J=1
C
C  SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST
C  APPROXIMATION
C
  200 IF(J.EQ.NZZ) YNZ=COMP
      IF(J.GE.NZZ) GO TO 300
      KUP=IEXT(J+1)
      L=IEXT(J)+1
      NUT=-NUT
      IF(J.EQ.2) Y1=COMP
      COMP=DEV
      IF(L.GE.KUP) GO TO 220
      ERR=GEED(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=FLOAT(NUT)*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 220
      COMP=FLOAT(NUT)*ERR
  210 L=L+1
      IF(L.GE.KUP) GO TO 215
      ERR=GEED(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=FLOAT(NUT)*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 215
      COMP=FLOAT(NUT)*ERR
      GO TO 210
  215 IEXT(J)=L-1
      J=J+1
      KLOW=L-1
      JCHNGE=JCHNGE+1
      GO TO 200
  220 L=L-1
  225 L=L-1
      IF(L.LE.KLOW) GO TO 250
      ERR=GEED(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=FLOAT(NUT)*ERR-COMP
      IF(DTEMP.GT.0.0) GO TO 230
      IF(JCHNGE.LE.0) GO TO 225
      GO TO 260
  230 COMP=FLOAT(NUT)*ERR
  235 L=L-1
      IF(L.LE.KLOW) GO TO 240
      ERR=GEED(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=FLOAT(NUT)*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 240
      COMP=FLOAT(NUT)*ERR
      GO TO 235
  240 KLOW=IEXT(J)
      IEXT(J)=L+1
      J=J+1
      JCHNGE=JCHNGE+1
      GO TO 200
  250 L=IEXT(J)+1
      IF(JCHNGE.GT.0) GO TO 215
  255 L=L+1
      IF(L.GE.KUP) GO TO 260
      ERR=GEED(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=FLOAT(NUT)*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 255
      COMP=FLOAT(NUT)*ERR
      GO TO 210
  260 KLOW=IEXT(J)
      J=J+1
      GO TO 200
  300 IF(J.GT.NZZ) GO TO 320
      IF(K1.GT.IEXT(1)) K1=IEXT(1)
      IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
      NUT1=NUT
      NUT=-NU
      L=0
      KUP=K1
      COMP=YNZ*(1.00001)
      LUCK=1
  310 L=L+1
      IF(L.GE.KUP) GO TO 315
      ERR=GEED(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=FLOAT(NUT)*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 310
      COMP=FLOAT(NUT)*ERR
      J=NZZ
      GO TO 210
  315 LUCK=6
      GO TO 325
  320 IF(LUCK.GT.9) GO TO 350
      IF(COMP.GT.Y1) Y1=COMP
      K1=IEXT(NZZ)
  325 L=NGRID+1
      KLOW=KNZ
      NUT=-NUT1
      COMP=Y1*(1.00001)
  330 L=L-1
      IF(L.LE.KLOW) GO TO 340
      ERR=GEED(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=FLOAT(NUT)*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 330
      J=NZZ
      COMP=FLOAT(NUT)*ERR
      LUCK=LUCK+10
      GO TO 235
  340 IF(LUCK.EQ.6) GO TO 370
      DO 345 J=1,NFCNS
      NZZMJ=NZZ-J
      NZMJ=NZ-J
  345 IEXT(NZZMJ)=IEXT(NZMJ)
      IEXT(1)=K1
      GO TO 100
  350 KN=IEXT(NZZ)
      DO 360 J=1,NFCNS
  360 IEXT(J)=IEXT(J+1)
      IEXT(NZ)=KN
      GO TO 100
  370 IF(JCHNGE.GT.0) GO TO 100
C
C  CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
C  USING THE INVERSE DISCRETE FOURIER TRANSFORM
C
  400 CONTINUE
      NM1=NFCNS-1
      FSH=1.0E-06
      GTEMP=GRID(1)
      X(NZZ)=-2.0
      CN=2*NFCNS-1
      DELF=1.0/CN
      L=1
      KKK=0
      IF(GRID(1).LT.0.01.AND.GRID(NGRID).GT.0.49) KKK=1
      IF(NFCNS.LE.3) KKK=1
      IF(KKK.EQ.1) GO TO 405
      DTEMP=DCOS(PI2*GRID(1))
      DNUM=DCOS(PI2*GRID(NGRID))
      AA=2.0/(DTEMP-DNUM)
      BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
  405 CONTINUE
      DO 430 J=1,NFCNS
      FT=J-1
      FT=FT*DELF
      XT=DCOS(PI2*FT)
      IF(KKK.EQ.1) GO TO 410
      XT=(XT-BB)/AA
      XT1=SQRT(1.0-XT*XT)
      FT=ATAN2(XT1,XT)/PI2
  410 XE=X(L)
      IF(XT.GT.XE) GO TO 420
      IF((XE-XT).LT.FSH) GO TO 415
      L=L+1
      GO TO 410
  415 A(J)=Y(L)
      GO TO 425
  420 IF((XT-XE).LT.FSH) GO TO 415
      GRID(1)=FT
      A(J)=GEED(1,NZ)
  425 CONTINUE
      IF(L.GT.1) L=L-1
  430 CONTINUE
      GRID(1)=GTEMP
      DDEN=PI2/CN
      DO 510 J=1,NFCNS
      DTEMP=0.0
      DNUM=J-1
      DNUM=DNUM*DDEN
      IF(NM1.LT.1) GO TO 505
      DO 500 K=1,NM1
      DAK=A(K+1)
      DK=K
  500 DTEMP=DTEMP+DAK*DCOS(DNUM*DK)
  505 DTEMP=2.0*DTEMP+A(1)
  510 ALPHA(J)=DTEMP
      DO 550 J=2,NFCNS
  550 ALPHA(J)=2.0*ALPHA(J)/CN
      ALPHA(1)=ALPHA(1)/CN
      IF(KKK.EQ.1) GO TO 545
      P(1)=2.0*ALPHA(NFCNS)*BB+ALPHA(NM1)
      P(2)=2.0*AA*ALPHA(NFCNS)
      Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
      DO 540 J=2,NM1
      IF(J.LT.NM1) GO TO 515
      AA=0.5*AA
      BB=0.5*BB
  515 CONTINUE
      P(J+1)=0.0
      DO 520 K=1,J
      A(K)=P(K)
  520 P(K)=2.0*BB*A(K)
      P(2)=P(2)+A(1)*2.0*AA
      JM1=J-1
      DO 525 K=1,JM1
  525 P(K)=P(K)+Q(K)+AA*A(K+1)
      JP1=J+1
      DO 530 K=3,JP1
  530 P(K)=P(K)+AA*A(K-1)
      IF(J.EQ.NM1) GO TO 540
      DO 535 K=1,J
  535 Q(K)=-A(K)
      NF1J=NFCNS-1-J
      Q(1)=Q(1)+ALPHA(NF1J)
  540 CONTINUE
      DO 543 J=1,NFCNS
  543 ALPHA(J)=P(J)
  545 CONTINUE
      IF(NFCNS.GT.3) RETURN
      ALPHA(NFCNS+1)=0.0
      ALPHA(NFCNS+2)=0.0
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: DD
C   FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION
C   COEFFICIENTS FOR USE IN THE FUNCTION GEED.
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION FUNCTION DD(K,N,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
      DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
      DIMENSION DES(1045),GRID(1045),WT(1045)
      DOUBLE PRECISION AD,DEV,X,Y
      DOUBLE PRECISION Q
      DOUBLE PRECISION PI2
      DD=1.0
      Q=X(K)
      DO 3 L=1,M
      DO 2 J=L,N,M
      IF(J-K)1,2,1
    1 DD=2.0*DD*(Q-X(J))
    2 CONTINUE
    3 CONTINUE
      DD=1.0/DD
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION: GEED
C   FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE
C   LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION FUNCTION GEED(K,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
      DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
      DIMENSION DES(1045),GRID(1045),WT(1045)
      DOUBLE PRECISION P,C,D,XF
      DOUBLE PRECISION PI2
      DOUBLE PRECISION AD,DEV,X,Y
      P=0.0
      XF=GRID(K)
      XF=DCOS(PI2*XF)
      D=0.0
      DO 1 J=1,N
      C=XF-X(J)
      C=AD(J)/C
      D=D+C
    1 P=P+C*Y(J)
      GEED=P/D
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUCH
C   WRITES AN ERROR MESSAGE WHEN THE ALGORITHM FAILS TO
C   CONVERGE.  THERE SEEM TO BE TWO CONDITIONS UNDER WHICH
C   THE ALGORITHM FAILS TO CONVERGE: (1) THE INITIAL
C   GUESS FOR THE EXTREMAL FREQUENCIES IS SO POOR THAT
C   THE EXCHANGE ITERATION CANNOT GET STARTED, OR
C   (2) NEAR THE TERMINATION OF A CORRECT DESIGN,
C   THE DEVIATION DECREASES DUE TO ROUNDING ERRORS
C   AND THE PROGRAM STOPS.  IN THIS LATTER CASE THE
C   FILTER DESIGN IS PROBABLY ACCEPTABLE, BUT SHOULD
C   BE CHECKED BY COMPUTING A FREQUENCY RESPONSE.
C-----------------------------------------------------------------------
C
      SUBROUTINE OUCH
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /OOPS/NITER
C *      WRITE(*,1)NITER
C *    1 FORMAT(44H ************ FAILURE TO CONVERGE **********/
C *     141H0PROBABLE CAUSE IS MACHINE ROUNDING ERROR/
C *     223H0NUMBER OF ITERATIONS =,I4/
C *     339H0IF THE NUMBER OF ITERATIONS EXCEEDS 3,/
C *     462H0THE DESIGN MAY BE CORRECT, BUT SHOULD BE VERIFIED WITH AN FFT)
      RETURN
      END
