C @(#)OPTFIR.F    1.1    11/2/90
C
C   THIS PROGRAM DESIGNS OPTIMAL EQUIRRIPLE FIR DIGITAL FILTERS
C
      DOUBLE PRECISION ATTN,RIPPLE,SVRIP
      DOUBLE PRECISION PI2,PI,AD,DEV,X,Y
      CHARACTER FNAME*100, ANSWER*1
      COMMON DES,WT,ALPHA,IEXT,NFCNS,NGRID,PI2,AD,DEV,X,Y,GRID
      DIMENSION IEXT(100),AD(100),ALPHA(100),X(100),Y(100),H(100)
      DIMENSION DES(1600),GRID(1600),WT(1600)
      DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
      PI=3.141592653589793
      PI2=6.283185307179586
C MAXIMUM FILTER LENGTH
      NFMAX=256
      WRITE (*,'(,A)') ' PROVISIONAL EQUIRRIPLE FILTER DESIGN '
      WRITE (*,'(,A)') ' USE AT YOUR OWN RISK --------------- '
  100 CONTINUE
      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
         IO = 5
      ELSE
         IO = 3
         OPEN (3, FILE=FNAME)
      ENDIF
      WRITE (*,'(,A,/,20X,A)') ' ENTER FILTER TYPE (1=BANDPASS, 2=DIFFER
     +ENTIATOR,', '3=HILBERT TRANSFORMER, 4=HALF-BAND): '
      READ (IO,*) JTYPE
      WRITE (*,'(,A)') ' ENTER FILTER LENGTH (ENTER 0 FOR ESTIMATE): '
      READ (IO,*) NFILT
      IF(NFILT.GT.NFMAX.OR.(NFILT.LT.3.AND.NFILT.NE.0)) GO TO 515
      WRITE (*,'(,A)') ' ENTER SAMPLING RATE OF FILTER: '
      READ (IO,*) FS
      IF (JTYPE .EQ. 4) THEN
         NFILT = (NFILT + 1) / 2
         NBANDS = 1
         EDGE(1) = 0.0
         WRITE (*,'(/,A\)') ' ENTER PASSBAND EDGE FREQUENCY: '
         READ (IO,*) EDGE(2)
         EDGE(2) = 2.0 * EDGE(2) / FS
         FX(1) = 0.5
         WTX(1) = 1.0
         GO TO 118
      ENDIF
      WRITE (*,'(,A)') ' ENTER NUMBER OF FILTER BANDS: '
      READ (IO,*) NBANDS
      IF(NBANDS.LE.0) NBANDS=1
      DO 115 I=1,NBANDS
      WRITE (*,'(,A,I3,A)') ' ENTER LOWER BAND EDGE FOR BAND', I, ': '
      J=1+2*(I-1)
      READ (IO,*) EDGE(J)
      EDGE(J)=EDGE(J)/FS
      WRITE (*,'(,A,I3,A)') ' ENTER UPPER BAND EDGE FOR BAND', I, ': '
      READ (IO,*) EDGE(J+1)
      EDGE(J+1)=EDGE(J+1)/FS
      WRITE (*,'(,A,I3,A)') ' ENTER DESIRED VALUE FOR BAND', I, ': '
      READ (IO,*) FX(I)
      WRITE (*,'(,A,I3,A)') ' ENTER WEIGHT FACTOR FOR BAND', I, ': '
      READ (IO,*) WTX(I)
  115 CONTINUE
      WRITE (*,'(,A)') ' DO YOU WANT X/SIN(X) PREDISTORTION? (Y/N): '
      READ (IO,'(A1)') ANSWER
      IF (ANSWER .EQ. 'Y' .OR. ANSWER .EQ. 'Y') IDIST=1
  118 WRITE (*,'(,A)') ' ENTER NAME OF COEFFICIENT OUTPUT FILE '
      WRITE (*,'(,A)') ' (SORRY, NO TILDE-EXPANSION.  GIVE PATH RELATIVE
     +VE TO YOUR HOME DIRECTORY): '
      READ (IO,'(A)') FNAME
      OPEN (1, FILE=FNAME)
      IF(NFILT.EQ.0) GO TO 120
      NROX=0
      GO TO 130
  120 CONTINUE
      WRITE (*,'(,A)') ' ENTER DESIRED PASSBAND RIPPLE IN DB: '
      READ (IO,*) RIPPLE
      WRITE (*,'(,A)') ' ENTER DESIRED STOPBAND ATTENUATION IN DB: '
      READ (IO,*) ATTN
      RIPPLE=10.**(RIPPLE/20.)
      ATTN=10.**(-ATTN/20.)
      RIPPLE=(RIPPLE-1)/(RIPPLE+1)
      RIPPLE=DLOG10(RIPPLE)
      ATTN=DLOG10(ATTN)
      SVRIP=RIPPLE
      RIPPLE=(0.005309*(RIPPLE**2)+0.07114*RIPPLE-0.4761)*ATTN
     1-(0.00266*(RIPPLE**2)+0.5941*RIPPLE+0.4278)
      ATTN=11.01217+0.51244*(SVRIP-ATTN)
      DELTAF=EDGE(3)-EDGE(2)
      NFILT=RIPPLE/DELTAF-ATTN*DELTAF+1
      WRITE (*,'(/,A,I4)') ' ESTIMATED FILTER LENGTH = ', NFILT
      IF(NFILT.GT.NFMAX.OR.(NFILT.LT.3.AND.NFILT.NE.0)) GO TO 515
      NROX=1
  130 NEG=1
      WRITE (*,'(/,A,/)') ' EXECUTING ...'
      IF(JTYPE .EQ. 1 .OR. JTYPE .EQ. 4) NEG = 0
      NODD=NFILT/2
      NODD=NFILT-2*NODD
      NFCNS=NFILT/2
      IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS=NFCNS+1
      GRID(1)=EDGE(1)
      DELF=16*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)
      DES(J)=EFF(TEMP,FX,LBAND,JTYPE,IDIST)
      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,LBAND,JTYPE,IDIST)
      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
      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
  200 TEMP=FLOAT(NGRID-1)/FLOAT(NFCNS)
      DO 210 J=1,NFCNS
  210 IEXT(J)=(J-1)*TEMP+1
      IEXT(NFCNS+1)=NGRID
      NM1=NFCNS-1
      NZ=NFCNS+1
      CALL REMEZ(EDGE,NBANDS)
      IF(NEG) 300,300,320
  300 IF(NODD.EQ.0) GO TO 310
      DO 305 J=1,NM1
  305 H(J)=0.5*ALPHA(NZ-J)
      H(NFCNS)=ALPHA(1)
      GO TO 350
  310 H(1)=0.25*ALPHA(NFCNS)
      DO 315 J=2,NM1
  315 H(J)=0.25*(ALPHA(NZ-J)+ALPHA(NFCNS+2-J))
      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
  325 H(J)=0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+3-J))
      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
  335 H(J)=0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))
      H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(2)
  350 WRITE(6,360)
  360 FORMAT(/' FINITE IMPULSE RESPONSE (FIR)'/
     1,' LINEAR PHASE DIGITAL FILTER DESIGN'/
     2,' REMEZ EXCHANGE ALGORITHM'/)
      IF(JTYPE.EQ.1) WRITE(6,365)
  365 FORMAT(' BANDPASS FILTER'/)
      IF(JTYPE.EQ.2) WRITE(6,370)
  370 FORMAT(' DIFFERENTIATOR'/)
      IF(JTYPE.EQ.3) WRITE(6,375)
  375 FORMAT(' HILBERT TRANSFORMER'/)
      IF(JTYPE.EQ.4) WRITE(6,376)
  376 FORMAT(' HALF-BAND FILTER'/)
      WRITE(6,378) NFILT
  378 FORMAT(' FILTER LENGTH = ',I3/)
      IF(NROX.EQ.1) WRITE(6,379)
  379 FORMAT(' FILTER LENGTH DETERMINED BY APPROXIMATION'/)
      WRITE(6,380)
  380 FORMAT(' IMPULSE RESPONSE COEFFICIENTS:'/)
      DO 381 J=1,NFCNS
      K=NFILT+1-J
      IF(NEG.EQ.0) WRITE(6,382)J,H(J),K
      IF(NEG.EQ.1) WRITE(6,383)J,H(J),K
  381 CONTINUE
  382 FORMAT(10X,'H(',I3,') = ', E14.7,' = H(',I4,')')
  383 FORMAT(10X,'H(',I3,') = ', E14.7,' = -H(',I4,')')
      IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(6,384) NZ, 0.0
  384 FORMAT(10X,'H(',I3,') = ', E14.7)
      WRITE (*,*)
C
C       DEBUG
C
C     PAUSE 'PRESS <ENTER> TO CONTINUE ...'
      DO 450 K=1,NBANDS,4
      KUP=K+3
      IF(KUP.GT.NBANDS) KUP=NBANDS
      WRITE(6,385)(J,J=K,KUP)
  385 FORMAT(25X,4('BAND',I3,8X))
      WRITE(6,390)(EDGE(2*J-1),J=K,KUP)
  390 FORMAT(/' LOWER BAND EDGE:',5F15.7)
      WRITE(6,395)(EDGE(2*J),J=K,KUP)
  395 FORMAT(' UPPER BAND EDGE:',5F15.7)
      IF(JTYPE.NE.2) WRITE(6,400) (FX(J),J=K,KUP)
  400 FORMAT(' DESIRED VALUE:',2X,5F15.7)
      IF(JTYPE.EQ.2) WRITE(6,405) (FX(J),J=K,KUP)
  405 FORMAT(' DESIRED SLOPE:',2X,5F15.7)
      WRITE(6,410) (WTX(J),J=K,KUP)
  410 FORMAT(' WEIGHT FACTOR:',2X,5F15.7)
      DO 420 J=K,KUP
  420 DEVIAT(J)=DEV/WTX(J)
      WRITE(6,425) (DEVIAT(J),J=K,KUP)
  425 FORMAT(' DEVIATION:',6X,5F15.7)
      DO 430 J=K,KUP
      IF(FX(J).EQ.1.0) DEVIAT(J)=(1.0+DEVIAT(J))/(1.0-DEVIAT(J))
  430 DEVIAT(J)=20.0*ALOG10(DEVIAT(J))
      WRITE(6,435) (DEVIAT(J),J=K,KUP)
  435 FORMAT(' DEVIATION IN DB:',5F15.7)
  450 CONTINUE
      WRITE(6,455) (GRID(IEXT(J)),J=1,NZ)
  455 FORMAT(/' EXTREMAL FREQUENCIES:'//(2X,5F12.7))
      IF (NODD .EQ. 1) THEN
         IF (NEG .EQ. 1) NFCNS = NFCNS + 1
         NCOEFF = 2 * NFCNS - 1
      ELSE
        NCOEFF = 2 * NFCNS
      ENDIF
      IF (JTYPE .EQ. 4) THEN
C
C   REMOVE THE WRITING OF THE FILTER LENGTH, FOR
C   COMPATIBILITY WITH GABRIEL
C        WRITE(1,500) 2 * NCOEFF - 1
C
  500    FORMAT(I4)
         DO 502 I = 1, NCOEFF - 1
            IF (I .LT. NFCNS) THEN
               WRITE (1,510) H(I)
               WRITE (1,510) 0.0
            ELSEIF (I .EQ. NFCNS) THEN
               WRITE (1,510) H(I)
               WRITE (1,510) 0.5
               WRITE (1,510) H(I)
            ELSE
               WRITE (1,510) 0.0
               WRITE (1,510) H(NCOEFF-I)
            ENDIF
  502    CONTINUE
      ELSE
C
C   REMOVE THE WRITING OF THE FILTER LENGTH, FOR
C   COMPATIBILITY WITH GABRIEL
C        WRITE(1,500) NCOEFF
C
         DO 505 I = 1, NCOEFF
            IF(I .LE. NFCNS) WRITE(1,510) H(I)
            IF(I .GT. NFCNS) THEN
               IF(NEG .EQ. 0) THEN
                  WRITE(1,510) H(NCOEFF-I+1)
               ELSE
                  WRITE(1,510) -H(NCOEFF-I+1)
               ENDIF
            ENDIF
  505    CONTINUE
      ENDIF
  510 FORMAT(1PE14.6)
      WRITE (*,*)
      GO TO 520
C
C   CATCH ERRORS
C
  515 PRINT *, '*** FILTER TOO LONG OR SHORT FOR PROGRAM ***'
C
C   FIND OUT WHETHER TO TERMINATE PROGRAM
C
  520 CALL FLUSH(1)
      CLOSE(1)
      WRITE (*,'(/,A\)') ' ANOTHER DESIGN? (Y/N): '
      READ (*,'(A1)') ANSWER
      IF (ANSWER .EQ. 'Y' .OR. ANSWER .EQ. 'Y') GO TO 100
      STOP 'END OF PROGRAM'
      END
C
      FUNCTION EFF(FREQ,FX,LBAND,JTYPE,IDIST)
      DIMENSION FX(10)
      IF(JTYPE.EQ.2) GO TO 1
      EFF=FX(LBAND)
      IF(IDIST.EQ.0) RETURN
      X=3.141593*FREQ
      IF(X.NE.0) EFF=EFF*X/SIN(X)
      RETURN
    1 EFF=FX(LBAND)*FREQ
      RETURN
      END
C
      FUNCTION WATE(FREQ,FX,WTX,LBAND,JTYPE)
      DIMENSION FX(10),WTX(10)
      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
      SUBROUTINE REMEZ(EDGE,NBANDS)
      DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
      DOUBLE PRECISION D,AD,DEV,X,Y
      COMMON DES,WT,ALPHA,IEXT,NFCNS,NGRID,PI2,AD,DEV,X,Y,GRID
      DIMENSION EDGE(20)
      DIMENSION IEXT(100),AD(100),ALPHA(100),X(100),Y(100)
      DIMENSION DES(1600),GRID(1600),WT(1600)
      DIMENSION A(100),P(100),Q(100)
      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
      DTEMP=GRID(IEXT(J))
      DTEMP=DCOS(DTEMP*PI2)
  110 X(J)=DTEMP
      JET=(NFCNS-1)/15+1
      DO 120 J=1,NZ
      KID=J
  120 AD(J)=D(KID,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=K*AD(J)/WT(L)
      DDEN=DDEN+DTEMP
  130 K=-K
      DEV=DNUM/DDEN
C
C       DEBUG
C
C     WRITE(6,131) DEV
C 131 FORMAT(1X,12HDEVIATION = ,F12.9)
C
      NU=1
      IF(DEV.GT.0.0) NU=-1
      DEV=-NU*DEV
      K=NU
      DO 140 J=1,NZ
      L=IEXT(J)
      DTEMP=K*DEV/WT(L)
      Y(J)=DES(L)+DTEMP
  140 K=-K
      IF(DEV.GE.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
  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=GEE(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=NUT*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 220
      COMP=NUT*ERR
  210 L=L+1
      IF(L.GE.KUP) GO TO 215
      ERR=GEE(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=NUT*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 215
      COMP=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=GEE(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=NUT*ERR-COMP
      IF(DTEMP.GT.0.0) GO TO 230
      IF(JCHNGE.LE.0) GO TO 225
      GO TO 260
  230 COMP=NUT*ERR
  235 L=L-1
      IF(L.LE.KLOW) GO TO 240
      ERR=GEE(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=NUT*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 240
      COMP=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.GE.0) GO TO 215
  255 L=L+1
      IF(L.GE.KUP) GO TO 260
      ERR=GEE(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=NUT*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 255
      COMP=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=GEE(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=NUT*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 310
      COMP=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=GEE(L,NZ)
      ERR=(ERR-DES(L))*WT(L)
      DTEMP=NUT*ERR-COMP
      IF(DTEMP.LE.0.0) GO TO 330
      J=NZZ
      COMP=NUT*ERR
      LUCK=LUCK+10
      GO TO 235
  340 IF(LUCK.EQ.6) GO TO 370
      DO 345 J=1,NFCNS
  345 IEXT(NZZ-J)=IEXT(NZ-J)
      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
  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(EDGE(1).EQ.0.0.AND.EDGE(2*NBANDS).EQ.0.5) 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)*DELF
      XT=DCOS(PI2*FT)
      IF(KKK.EQ.1) GO TO 410
      XT=(XT-BB)/AA
      FT=ACOS(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)=GEE(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)*DDEN
      IF(NM1.LT.1) GO TO 505
      DO 500 K=1,NM1
  500 DTEMP=DTEMP+A(K+1)*DCOS(DNUM*K)
  505 DTEMP=2.0*DTEMP+A(1)
  510 ALPHA(J)=DTEMP
      DO 550 J=2,NFCNS
  550 ALPHA(J)=2.*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)
      Q(1)=Q(1)+ALPHA(NFCNS-1-J)
  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
      DOUBLE PRECISION FUNCTION D(K,N,M)
      DOUBLE PRECISION AD,DEV,X,Y,Q,PI2
      COMMON DES,WT,ALPHA,IEXT,NFCNS,NGRID,PI2,AD,DEV,X,Y,GRID
      DIMENSION IEXT(100),AD(100),ALPHA(100),X(100),Y(100)
      DIMENSION DES(1600),GRID(1600),WT(1600)
      D=1.0
      Q=X(K)
      DO 3 L=1,M
      DO 2 J=L,N,M
      IF(J-K)1,2,1
    1 D=2.0*D*(Q-X(J))
    2 CONTINUE
    3 CONTINUE
      D=1.0/D
      RETURN
      END
C
      FUNCTION GEE(K,N)
      DOUBLE PRECISION P,C,D,XF,PI2,AD,DEV,X,Y
      COMMON DES,WT,ALPHA,IEXT,NFCNS,NGRID,PI2,AD,DEV,X,Y,GRID
      DIMENSION IEXT(100),AD(100),ALPHA(100),X(100),Y(100)
      DIMENSION DES(1600),GRID(1600),WT(1600)
      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)
      GEE=P/D
      RETURN
      END
C
      SUBROUTINE OUCH
C     WRITE (*,'(/A)') ' ********** FAILURE TO CONVERGE **********'
C     WRITE (*,*) 'PROBABLE CAUSE IS MACHINE ROUNDING ERROR'
C     WRITE (*,*) 'THE IMPULSE RESPONSE MAY STILL BE CORRECT'
      RETURN
      END