C @(#)REMEZ.F
C MAIN PROGRAM: FIR LINEAR PHASE FILTER DESIGN PROGRAM
C
C AUTHORS: JAMES H. MCCLELLAN
C
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:        NFILT IS THE FILTER LENGTH
C
C               JTYPE IS THE TYPE OF FILTER
C                   JTYPE = 1     MULTIPLE PASSBAND/STOPBAND FILTER
C                   JTYPE = 2     DIFFERENTIATOR
C                   JTYPE = 3     HILBERT TRANSFORM FILTER
C
C               NBANDS IS THE NUMBER OF BANDS
C
C               LGRID IS THE GRID DENSITY, WILL BE SET TO 16 UNLESS
C                   SPECIFIED OTHERWISE BY A POSITIVE CONSTANT.
C
C               EDGE(2*NBANDS) IS THE BANDEDGE ARRAY, LOWER AND UPPER
C                   EDGES FOR EACH BAND WITH A MAXIMUM OF 10 BANDS.
C
C               FX(NBANDS) IS THE DESIRED FUNCTION ARRAY (OR DESIRED 
C                   SLOPE IF A DIFFERENTIATOR) FOR EACH BAND.
C
C               WTX(NBANDS) IS THE WEIGHT FUNCTION ARRAY IN EACH BAND.  
C                   FOR A DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY
C                   PROPORTIONAL TO F.
C
C
C    SAMPLE INPUT DATA SETUP:
C
C    THE FOLLOWING INPUT 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               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
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
        COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
        COMMON /OOPS/ NITER,IOUT
        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 GEE,D
        INTEGER BD1,BD2,BD3,BD4
        DATA BD1,BD2,BD3,BD4/1HB,1HA,1HN,1HD/
        INPUT=6
        IOUT=6
        PI=4.0*DATAN(1.0D0)
        PI2=2.0D0*PI

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



        FUNCTION EFF(FREQ,FX,WTX,LBAND,JTYPE)
C
C       DIMENSION FX(10), WTX(10)
        DIMENSION FX(5), WTX(5)
        IF (JTYPE.EQ.2) GO TO 1
        EFF=FX(LBAND)
        RETURN
1       EFF=FX(LBAND)*FREQ
        RETURN
        END


        FUNCTION WATE(FREQ,FX,WTX,LBAND,JTYPE)
C
C       DIMENSION FX(10),WTX(10)
        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


        SUBROUTINE ERROR
        COMMON /OOPS/NITER,IOUT
        WRITE (IOUT,1)
1       FORMAT(25X,44H ********* ERROR IN INPUT DATA ************/)
        RETURN
        END

        SUBROUTINE REMEZ
C
        COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
        COMMON /OOPS/NITER,IOUT
        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
C       DOUBLE PRECISION GEE,D
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)=D(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
        WRITE (IOUT,131) DEV
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
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,1)
        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=GEE(L,NZ,2)
        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=GEE(L,NZ,3)
        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=GEE(L,NZ,4)
        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=GEE(L,NZ,5)
        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=GEE(L,NZ,6)
        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=GEE(L,NZ,7)
        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
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)=GEE(1,NZ,8)
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

        DOUBLE PRECISION FUNCTION D(K,N,M)
C
        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
        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



        DOUBLE PRECISION FUNCTION GEE(K,N,KCALL)
C
        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)
        IF (C.EQ.0.0) PRINT 5,KCALL,K,J,N,XF,X(J),C,AD(J)
        C=AD(J)/C
        D=D+C
1       P=P+C*Y(J)
        GEE=P/D
        RETURN
5       FORMAT ('KCALL,K,J;N;XF,X(J),C,AD(J)',I3,I3,I3,I3,D15.6,D15.6,
     1D15.6,D15.6)
        END


        SUBROUTINE OUCH
        COMMON /OOPS/NITER,IOUT
        WRITE (IOUT,1) NITER
1       FORMAT(2X,31H***** FAILURE TO CONVERGE *****/
     140HPROBABLE CAUSE IS MACHINE ROUNDING ERROR/
     223HNUMBER OF ITERATIONS = ,I4/
     337HIF THE NUMBER OF ITERATIONS EXCEEDS 3/
     460HTHE DESIGN MAY BE CORRECT BUT SHOULD BE VERIFIED WITH AN FFT)
        RETURN
        END