home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / REMEZF.F < prev    next >
Encoding:
Text File  |  1992-12-07  |  19.5 KB  |  716 lines
  1. C
  2. C-----------------------------------------------------------------------
  3. C MAIN PROGRAM: FIR LINEAR PHASE FILTER DESIGN PROGRAM
  4. C
  5. C AUTHORS: JAMES H. MCCLELLAN
  6. C          DEPARTMENT OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE
  7. C          MASSACHUSETTS INSTITUTE OF TECHNOLOGY
  8. C          CAMBRIDGE, MASS. 02139
  9. C
  10. C          THOMAS W. PARKS
  11. C          DEPARTMENT OF ELECTRICAL ENGINEERING
  12. C          RICE UNIVERSITY
  13. C          HOUSTON, TEXAS 77001
  14. C
  15. C          LAWRENCE R. RABINER
  16. C          BELL LABORATORIES
  17. C          MURRAY HILL, NEW JERSEY 07974
  18. C
  19. C INPUT:
  20. C  NFILT-- FILTER LENGTH
  21. C  JTYPE-- TYPE OF FILTER
  22. C          1 = MULTIPLE PASSBAND/STOPBAND FILTER
  23. C          2 = DIFFERENTIATOR
  24. C          3 = HILBERT TRANSFORM FILTER
  25. C  NBANDS-- NUMBER OF BANDS
  26. C  LGRID-- GRID DENSITY, WILL BE SET TO 16 UNLESS
  27. C          SPECIFIED OTHERWISE BY A POSITIVE CONSTANT.
  28. C
  29. C  EDGE(2*NBANDS)-- BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND
  30. C                   WITH A MAXIMUM OF 10 BANDS.
  31. C
  32. C  FX(NBANDS)-- DESIRED FUNCTION ARRAY (OR DESIRED SLOPE IF A
  33. C               DIFFERENTIATOR) FOR EACH BAND.
  34. C
  35. C  WTX(NBANDS)-- WEIGHT FUNCTION ARRAY IN EACH BAND.  FOR A
  36. C                DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY
  37. C                PROPORTIONAL TO F.
  38. C
  39. C  SAMPLE INPUT DATA SETUP:
  40. C       32,1,3,0
  41. C       0.0,0.1,0.2,0.35
  42. C    0.425,0.5
  43. C       0.0,1.0,0.0
  44. C       10.0,1.0,10.0
  45. C     THIS DATA SPECIFIES A LENGTH 32 BANDPASS FILTER WITH
  46. C     STOPBANDS 0 TO 0.1 AND 0.425 TO 0.5, AND PASSBAND FROM
  47. C     0.2 TO 0.35 WITH WEIGHTING OF 10 IN THE STOPBANDS AND 1
  48. C     IN THE PASSBAND.  THE GRID DENSITY DEFAULTS TO 16.
  49. C     THIS IS THE FILTER IN FIGURE 10.
  50. C
  51. C     THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 FULLBAND
  52. C     DIFFERENTIATOR WITH SLOPE 1 AND WEIGHTING OF 1/F.
  53. C     THE GRID DENSITY WILL BE SET TO 20.
  54. C       32,2,1,20
  55. C       0,0.5
  56. C       1.0
  57. C       1.0
  58. C
  59. C-----------------------------------------------------------------------
  60. C
  61.  
  62.       SUBROUTINE REMEZF ( NFILT, JTYPE, NBANDS, LGRID, EDGE, FX, WTX, H)
  63.  
  64.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  65.       DOUBLE PRECISION EDGE, FX, WTX, H
  66.  
  67.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  68.       COMMON /OOPS/NITER
  69.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  70.       DIMENSION H(66)
  71.       DIMENSION DES(1045),GRID(1045),WT(1045)
  72.       DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
  73.       DOUBLE PRECISION PI2,PI
  74.       DOUBLE PRECISION AD,DEV,X,Y
  75.       DOUBLE PRECISION GEED,DD
  76.       INTEGER BD1,BD2,BD3,BD4
  77.       DATA BD1,BD2,BD3,BD4/1HB,1HA,1HN,1HD/
  78.       PI=4.0*DATAN(1.0D0)
  79.       PI2=2.0D00*PI
  80. C
  81. C  THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128, BUT
  82. C  THIS UPPER LIMIT CAN BE CHANGED BY REDIMENSIONING THE
  83. C  ARRAYS IEXT, AD, ALPHA, X, Y, H TO BE NFMAX/2 + 2.
  84. C  THE ARRAYS DES, GRID, AND WT MUST DIMENSIONED
  85. C  16(NFMAX/2 + 2).
  86. C
  87.       NFMAX=128
  88.   100 CONTINUE
  89. C *      JTYPE=0
  90. C
  91. C  PROGRAM INPUT SECTION
  92. C
  93. C *      READ(*,110) NFILT,JTYPE,NBANDS,LGRID
  94.       IF(NFILT.EQ.0)STOP
  95.   110 FORMAT(4I5)
  96.       IF(NFILT.LE.NFMAX.AND.NFILT.GE.3) GO TO 115
  97.       CALL ERR51
  98.       STOP
  99.   115 IF(NBANDS.LE.0) NBANDS=1
  100. C
  101. C  GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED
  102. C  OTHERWISE
  103. C
  104.       IF(LGRID.LE.0) LGRID=16
  105.       JB=2*NBANDS
  106. C *      READ(*,120) (EDGE(J),J=1,JB)
  107. C *  120 FORMAT(4F15.9)
  108. C *      READ(*,120) (FX(J),J=1,NBANDS)
  109. C *      READ(*,120) (WTX(J),J=1,NBANDS)
  110.       IF(JTYPE.GT.0.AND.JTYPE.LE.3) GO TO 125
  111.       CALL ERR51
  112.       STOP
  113.   125 NEG=1
  114.       IF(JTYPE.EQ.1) NEG=0
  115.       NODD=NFILT/2
  116.       NODD=NFILT-2*NODD
  117.       NFCNS=NFILT/2
  118.       IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS=NFCNS+1
  119. C
  120. C  SET UP THE DENSE GRID.  THE NUMBER OF POINTS IN THE GRID
  121. C  IS (FILTER LENGTH + 1)*GRID DENSITY/2
  122. C
  123.       GRID(1)=EDGE(1)
  124.       DELF=LGRID*NFCNS
  125.       DELF=0.5/DELF
  126.       IF(NEG.EQ.0) GO TO 135
  127.       IF(EDGE(1).LT.DELF) GRID(1)=DELF
  128.   135 CONTINUE
  129.       J=1
  130.       L=1
  131.       LBAND=1
  132.   140 FUP=EDGE(L+1)
  133.   145 TEMP=GRID(J)
  134. C
  135. C  CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
  136. C  FUNCTION ON THE GRID
  137. C
  138.       DES(J)=EFF(TEMP,FX,WTX,LBAND,JTYPE)
  139.       WT(J)=WATE(TEMP,FX,WTX,LBAND,JTYPE)
  140.       J=J+1
  141.       GRID(J)=TEMP+DELF
  142.       IF(GRID(J).GT.FUP) GO TO 150
  143.       GO TO 145
  144.   150 GRID(J-1)=FUP
  145.       DES(J-1)=EFF(FUP,FX,WTX,LBAND,JTYPE)
  146.       WT(J-1)=WATE(FUP,FX,WTX,LBAND,JTYPE)
  147.       LBAND=LBAND+1
  148.       L=L+2
  149.       IF(LBAND.GT.NBANDS) GO TO 160
  150.       GRID(J)=EDGE(L)
  151.       GO TO 140
  152.   160 NGRID=J-1
  153.       IF(NEG.NE.NODD) GO TO 165
  154.       IF(GRID(NGRID).GT.(0.5-DELF)) NGRID=NGRID-1
  155.   165 CONTINUE
  156. C
  157. C  SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
  158. C  TO THE ORIGINAL PROBLEM
  159. C
  160.       IF(NEG) 170,170,180
  161.   170 IF(NODD.EQ.1) GO TO 200
  162.       DO 175 J=1,NGRID
  163.       CHANGE=DCOS(PI*GRID(J))
  164.       DES(J)=DES(J)/CHANGE
  165.   175 WT(J)=WT(J)*CHANGE
  166.       GO TO 200
  167.   180 IF(NODD.EQ.1) GO TO 190
  168.       DO 185 J=1,NGRID
  169.       CHANGE=DSIN(PI*GRID(J))
  170.       DES(J)=DES(J)/CHANGE
  171.   185 WT(J)=WT(J)*CHANGE
  172.       GO TO 200
  173.   190 DO 195 J=1,NGRID
  174.       CHANGE=DSIN(PI2*GRID(J))
  175.       DES(J)=DES(J)/CHANGE
  176.   195 WT(J)=WT(J)*CHANGE
  177. C
  178. C  INITIAL GUESS FOR THE EXTREMAL FREQUENCIES--EQUALLY
  179. C  SPACED ALONG THE GRID
  180. C
  181.   200 TEMP=FLOAT(NGRID-1)/FLOAT(NFCNS)
  182.       DO 210 J=1,NFCNS
  183.       XT=J-1
  184.   210 IEXT(J)=XT*TEMP+1.0
  185.       IEXT(NFCNS+1)=NGRID
  186.       NM1=NFCNS-1
  187.       NZ=NFCNS+1
  188. C
  189. C  CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION
  190. C  PROBLEM
  191. C
  192.       CALL REMEZ
  193. C
  194. C  CALCULATE THE IMPULSE RESPONSE.
  195. C
  196.       IF(NEG) 300,300,320
  197.   300 IF(NODD.EQ.0) GO TO 310
  198.       DO 305 J=1,NM1
  199.       NZMJ=NZ-J
  200.   305 H(J)=0.5*ALPHA(NZMJ)
  201.       H(NFCNS)=ALPHA(1)
  202.       GO TO 350
  203.   310 H(1)=0.25*ALPHA(NFCNS)
  204.       DO 315 J=2,NM1
  205.       NZMJ=NZ-J
  206.       NF2J=NFCNS+2-J
  207.   315 H(J)=0.25*(ALPHA(NZMJ)+ALPHA(NF2J))
  208.       H(NFCNS)=0.5*ALPHA(1)+0.25*ALPHA(2)
  209.       GO TO 350
  210.   320 IF(NODD.EQ.0) GO TO 330
  211.       H(1)=0.25*ALPHA(NFCNS)
  212.       H(2)=0.25*ALPHA(NM1)
  213.       DO 325 J=3,NM1
  214.       NZMJ=NZ-J
  215.       NF3J=NFCNS+3-J
  216.   325 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF3J))
  217.       H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(3)
  218.       H(NZ)=0.0
  219.       GO TO 350
  220.   330 H(1)=0.25*ALPHA(NFCNS)
  221.       DO 335 J=2,NM1
  222.       NZMJ=NZ-J
  223.       NF2J=NFCNS+2-J
  224.   335 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF2J))
  225.       H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(2)
  226. C
  227. C  PROGRAM OUTPUT SECTION.
  228. C
  229.  
  230.   350 GO TO 470
  231.  
  232. C *  350 WRITE(*,360)
  233. C *  360 FORMAT(1H1, 70(1H*)//15X,29HFINITE IMPULSE RESPONSE (FIR)/
  234. C *     113X,34HLINEAR PHASE DIGITAL FILTER DESIGN/
  235. C *     217X,24HREMEZ EXCHANGE ALGORITHM/)
  236. C *      IF(JTYPE.EQ.1) WRITE(*,365)
  237. C *  365 FORMAT(22X,15HBANDPASS FILTER/)
  238. C *      IF(JTYPE.EQ.2) WRITE(*,370)
  239. C *  370 FORMAT(22X,14HDIFFERENTIATOR/)
  240. C *      IF(JTYPE.EQ.3) WRITE(*,375)
  241. C *  375 FORMAT(20X,19HHILBERT TRANSFORMER/)
  242. C *      WRITE(*,378) NFILT
  243. C *  378 FORMAT(20X,16HFILTER LENGTH = ,I3/)
  244. C *      WRITE(*,380)
  245. C *  380 FORMAT(15X,28H***** IMPULSE RESPONSE *****)
  246. C *      DO 381 J=1,NFCNS
  247. C *      K=NFILT+1-J
  248. C *      IF(NEG.EQ.0) WRITE(*,382) J,H(J),K
  249. C *      IF(NEG.EQ.1) WRITE(*,383) J,H(J),K
  250. C *  381 CONTINUE
  251. C *  382 FORMAT(13X,2HH(,I2,4H) = ,E15.8,5H = H(,I3,1H))
  252. C *  383 FORMAT(13X,2HH(,I2,4H) = ,E15.8,6H = -H(,I3,1H))
  253. C *      IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(*,384) NZ
  254. C *  384 FORMAT(13X,2HH(,I2,8H) =  0.0)
  255. C *      DO 450 K=1,NBANDS,4
  256. C *      KUP=K+3
  257. C *      IF(KUP.GT.NBANDS) KUP=NBANDS
  258. C *      WRITE(*,385) (BD1,BD2,BD3,BD4,J,J=K,KUP)
  259. C *  385 FORMAT(/24X,4(4A1,I3,7X))
  260. C *      WRITE(*,390) (EDGE(2*J-1),J=K,KUP)
  261. C *  390 FORMAT(2X,15HLOWER BAND EDGE,5F14.7)
  262. C *      WRITE(*,395) (EDGE(2*J),J=K,KUP)
  263. C *  395 FORMAT(2X,15HUPPER BAND EDGE,5F14.7)
  264. C *      IF(JTYPE.NE.2) WRITE(*,400) (FX(J),J=K,KUP)
  265. C *  400 FORMAT(2X,13HDESIRED VALUE,2X,5F14.7)
  266. C *      IF(JTYPE.EQ.2) WRITE(*,405) (FX(J),J=K,KUP)
  267. C *  405 FORMAT(2X,13HDESIRED SLOPE,2X,5F14.7)
  268. C *      WRITE(*,410) (WTX(J),J=K,KUP)
  269. C *  410 FORMAT(2X,9HWEIGHTING,6X,5F14.7)
  270. C *      DO 420 J=K,KUP
  271. C *  420 DEVIAT(J)=DEV/WTX(J)
  272. C *      WRITE(*,425) (DEVIAT(J),J=K,KUP)
  273. C *  425 FORMAT(2X,9HDEVIATION,6X,5F14.7)
  274. C *      IF(JTYPE.NE.1) GO TO 450
  275. C *      DO 430 J=K,KUP
  276. C *  430 DEVIAT(J)=20.0*DLOG10(DEVIAT(J)+FX(J))
  277. C *      WRITE(*,435) (DEVIAT(J),J=K,KUP)
  278. C *  435 FORMAT(2X,15HDEVIATION IN DB,5F14.7)
  279. C *  450 CONTINUE
  280. C *      DO 452 J=1,NZ
  281. C *      IX=IEXT(J)
  282. C *  452 GRID(J)=GRID(IX)
  283. C *      WRITE(*,455) (GRID(J),J=1,NZ)
  284. C *  455 FORMAT(/2X,47HEXTREMAL FREQUENCIES--MAXIMA OF THE ERROR CURVE/
  285. C *     1 (2X,5F12.7))
  286. C *      WRITE(*,460)
  287. C *  460 FORMAT(/1X,70(1H*)/1H1)
  288. C *        GO TO 100
  289.  
  290.   470 RETURN
  291.  
  292.       END
  293. C
  294. C-----------------------------------------------------------------------
  295. C FUNCTION: EFF
  296. C   FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE
  297. C   AS A FUNCTION OF FREQUENCY.
  298. C   AN ARBITRARY FUNCTION OF FREQUENCY CAN BE
  299. C   APPROXIMATED IF THE USER REPLACES THIS FUNCTION
  300. C   WITH THE APPROPRIATE CODE TO EVALUATE THE IDEAL
  301. C   MAGNITUDE.  NOTE THAT THE PARAMETER FREQ IS THE
  302. C   VALUE OF NORMALIZED FREQUENCY NEEDED FOR EVALUATION.
  303. C-----------------------------------------------------------------------
  304. C
  305.       DOUBLE PRECISION FUNCTION EFF(FREQ,FX,WTX,LBAND,JTYPE)
  306.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  307.       DOUBLE PRECISION FX,WTX
  308.       DIMENSION FX(5),WTX(5)
  309.       IF(JTYPE.EQ.2) GO TO 1
  310.       EFF=FX(LBAND)
  311.       RETURN
  312.     1 EFF=FX(LBAND)*FREQ
  313.       RETURN
  314.       END
  315. C
  316. C-----------------------------------------------------------------------
  317. C FUNCTION: WATE
  318. C   FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION
  319. C   OF FREQUENCY.  SIMILAR TO THE FUNCTION EFF, THIS FUNCTION CAN
  320. C   BE REPLACED BY A USER-WRITTEN ROUTINE TO CALCULATE ANY
  321. C   DESIRED WEIGHTING FUNCTION.
  322. C-----------------------------------------------------------------------
  323. C
  324.       DOUBLE PRECISION FUNCTION WATE(FREQ,FX,WTX,LBAND,JTYPE)
  325.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  326.       DOUBLE PRECISION FX,WTX
  327.       DIMENSION FX(5),WTX(5)
  328.       IF(JTYPE.EQ.2) GO TO 1
  329.       WATE=WTX(LBAND)
  330.       RETURN
  331.     1 IF(FX(LBAND).LT.0.0001) GO TO 2
  332.       WATE=WTX(LBAND)/FREQ
  333.       RETURN
  334.     2 WATE=WTX(LBAND)
  335.       RETURN
  336.       END
  337. C
  338. C-----------------------------------------------------------------------
  339. C SUBROUTINE: ERR51
  340. C   THIS ROUTINE WRITES AN ERROR MESSAGE IF AN
  341. C   ERROR HAS BEEN DETECTED IN THE INPUT DATA.
  342. C-----------------------------------------------------------------------
  343. C
  344.       SUBROUTINE ERR51
  345.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  346.       COMMON /OOPS/NITER
  347. C *      WRITE(*,1)
  348. C *    1 FORMAT(44H ************ ERROR IN INPUT DATA **********)
  349.       RETURN
  350.       END
  351. C
  352. C-----------------------------------------------------------------------
  353. C SUBROUTINE: REMEZ
  354. C   THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
  355. C   FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS
  356. C   FUNCTION WITH A SUM OF COSINES.  INPUTS TO THE SUBROUTINE
  357. C   ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
  358. C   DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
  359. C   GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE
  360. C   EXTREMAL FREQUENCIES.  THE PROGRAM MINIMIZES THE CHEBYSHEV
  361. C   ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
  362. C   FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
  363. C   THE COEFFICIENTS OF THE BEST APPROXIMATION.
  364. C-----------------------------------------------------------------------
  365. C
  366.       SUBROUTINE REMEZ
  367.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  368.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  369.       COMMON /OOPS/NITER
  370.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  371.       DIMENSION DES(1045),GRID(1045),WT(1045)
  372.       DIMENSION A(66),P(65),Q(65)
  373.       DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
  374.       DOUBLE PRECISION DK,DAK
  375.       DOUBLE PRECISION AD,DEV,X,Y
  376.       DOUBLE PRECISION GEED,DD
  377. C
  378. C  THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
  379. C
  380.       ITRMAX=25
  381.       DEVL=-1.0
  382.       NZ=NFCNS+1
  383.       NZZ=NFCNS+2
  384.       NITER=0
  385.   100 CONTINUE
  386.       IEXT(NZZ)=NGRID+1
  387.       NITER=NITER+1
  388.       IF(NITER.GT.ITRMAX) GO TO 400
  389.       DO 110 J=1,NZ
  390.       JXT=IEXT(J)
  391.       DTEMP=GRID(JXT)
  392.       DTEMP=DCOS(DTEMP*PI2)
  393.   110 X(J)=DTEMP
  394.       JET=(NFCNS-1)/15+1
  395.       DO 120 J=1,NZ
  396.   120 AD(J)=DD(J,NZ,JET)
  397.       DNUM=0.0
  398.       DDEN=0.0
  399.       K=1
  400.       DO 130 J=1,NZ
  401.       L=IEXT(J)
  402.       DTEMP=AD(J)*DES(L)
  403.       DNUM=DNUM+DTEMP
  404.       DTEMP=FLOAT(K)*AD(J)/WT(L)
  405.       DDEN=DDEN+DTEMP
  406.   130 K=-K
  407.       DEV=DNUM/DDEN
  408. C *      WRITE(*,131) DEV
  409. C *  131 FORMAT(1X,12HDEVIATION = ,F12.9)
  410.       NU=1
  411.       IF(DEV.GT.0.0) NU=-1
  412.       DEV=-FLOAT(NU)*DEV
  413.       K=NU
  414.       DO 140 J=1,NZ
  415.       L=IEXT(J)
  416.       DTEMP=FLOAT(K)*DEV/WT(L)
  417.       Y(J)=DES(L)+DTEMP
  418.   140 K=-K
  419.       IF(DEV.GT.DEVL) GO TO 150
  420.       CALL OUCH
  421.       GO TO 400
  422.   150 DEVL=DEV
  423.       JCHNGE=0
  424.       K1=IEXT(1)
  425.       KNZ=IEXT(NZ)
  426.       KLOW=0
  427.       NUT=-NU
  428.       J=1
  429. C
  430. C  SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST
  431. C  APPROXIMATION
  432. C
  433.   200 IF(J.EQ.NZZ) YNZ=COMP
  434.       IF(J.GE.NZZ) GO TO 300
  435.       KUP=IEXT(J+1)
  436.       L=IEXT(J)+1
  437.       NUT=-NUT
  438.       IF(J.EQ.2) Y1=COMP
  439.       COMP=DEV
  440.       IF(L.GE.KUP) GO TO 220
  441.       ERR=GEED(L,NZ)
  442.       ERR=(ERR-DES(L))*WT(L)
  443.       DTEMP=FLOAT(NUT)*ERR-COMP
  444.       IF(DTEMP.LE.0.0) GO TO 220
  445.       COMP=FLOAT(NUT)*ERR
  446.   210 L=L+1
  447.       IF(L.GE.KUP) GO TO 215
  448.       ERR=GEED(L,NZ)
  449.       ERR=(ERR-DES(L))*WT(L)
  450.       DTEMP=FLOAT(NUT)*ERR-COMP
  451.       IF(DTEMP.LE.0.0) GO TO 215
  452.       COMP=FLOAT(NUT)*ERR
  453.       GO TO 210
  454.   215 IEXT(J)=L-1
  455.       J=J+1
  456.       KLOW=L-1
  457.       JCHNGE=JCHNGE+1
  458.       GO TO 200
  459.   220 L=L-1
  460.   225 L=L-1
  461.       IF(L.LE.KLOW) GO TO 250
  462.       ERR=GEED(L,NZ)
  463.       ERR=(ERR-DES(L))*WT(L)
  464.       DTEMP=FLOAT(NUT)*ERR-COMP
  465.       IF(DTEMP.GT.0.0) GO TO 230
  466.       IF(JCHNGE.LE.0) GO TO 225
  467.       GO TO 260
  468.   230 COMP=FLOAT(NUT)*ERR
  469.   235 L=L-1
  470.       IF(L.LE.KLOW) GO TO 240
  471.       ERR=GEED(L,NZ)
  472.       ERR=(ERR-DES(L))*WT(L)
  473.       DTEMP=FLOAT(NUT)*ERR-COMP
  474.       IF(DTEMP.LE.0.0) GO TO 240
  475.       COMP=FLOAT(NUT)*ERR
  476.       GO TO 235
  477.   240 KLOW=IEXT(J)
  478.       IEXT(J)=L+1
  479.       J=J+1
  480.       JCHNGE=JCHNGE+1
  481.       GO TO 200
  482.   250 L=IEXT(J)+1
  483.       IF(JCHNGE.GT.0) GO TO 215
  484.   255 L=L+1
  485.       IF(L.GE.KUP) GO TO 260
  486.       ERR=GEED(L,NZ)
  487.       ERR=(ERR-DES(L))*WT(L)
  488.       DTEMP=FLOAT(NUT)*ERR-COMP
  489.       IF(DTEMP.LE.0.0) GO TO 255
  490.       COMP=FLOAT(NUT)*ERR
  491.       GO TO 210
  492.   260 KLOW=IEXT(J)
  493.       J=J+1
  494.       GO TO 200
  495.   300 IF(J.GT.NZZ) GO TO 320
  496.       IF(K1.GT.IEXT(1)) K1=IEXT(1)
  497.       IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
  498.       NUT1=NUT
  499.       NUT=-NU
  500.       L=0
  501.       KUP=K1
  502.       COMP=YNZ*(1.00001)
  503.       LUCK=1
  504.   310 L=L+1
  505.       IF(L.GE.KUP) GO TO 315
  506.       ERR=GEED(L,NZ)
  507.       ERR=(ERR-DES(L))*WT(L)
  508.       DTEMP=FLOAT(NUT)*ERR-COMP
  509.       IF(DTEMP.LE.0.0) GO TO 310
  510.       COMP=FLOAT(NUT)*ERR
  511.       J=NZZ
  512.       GO TO 210
  513.   315 LUCK=6
  514.       GO TO 325
  515.   320 IF(LUCK.GT.9) GO TO 350
  516.       IF(COMP.GT.Y1) Y1=COMP
  517.       K1=IEXT(NZZ)
  518.   325 L=NGRID+1
  519.       KLOW=KNZ
  520.       NUT=-NUT1
  521.       COMP=Y1*(1.00001)
  522.   330 L=L-1
  523.       IF(L.LE.KLOW) GO TO 340
  524.       ERR=GEED(L,NZ)
  525.       ERR=(ERR-DES(L))*WT(L)
  526.       DTEMP=FLOAT(NUT)*ERR-COMP
  527.       IF(DTEMP.LE.0.0) GO TO 330
  528.       J=NZZ
  529.       COMP=FLOAT(NUT)*ERR
  530.       LUCK=LUCK+10
  531.       GO TO 235
  532.   340 IF(LUCK.EQ.6) GO TO 370
  533.       DO 345 J=1,NFCNS
  534.       NZZMJ=NZZ-J
  535.       NZMJ=NZ-J
  536.   345 IEXT(NZZMJ)=IEXT(NZMJ)
  537.       IEXT(1)=K1
  538.       GO TO 100
  539.   350 KN=IEXT(NZZ)
  540.       DO 360 J=1,NFCNS
  541.   360 IEXT(J)=IEXT(J+1)
  542.       IEXT(NZ)=KN
  543.       GO TO 100
  544.   370 IF(JCHNGE.GT.0) GO TO 100
  545. C
  546. C  CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
  547. C  USING THE INVERSE DISCRETE FOURIER TRANSFORM
  548. C
  549.   400 CONTINUE
  550.       NM1=NFCNS-1
  551.       FSH=1.0E-06
  552.       GTEMP=GRID(1)
  553.       X(NZZ)=-2.0
  554.       CN=2*NFCNS-1
  555.       DELF=1.0/CN
  556.       L=1
  557.       KKK=0
  558.       IF(GRID(1).LT.0.01.AND.GRID(NGRID).GT.0.49) KKK=1
  559.       IF(NFCNS.LE.3) KKK=1
  560.       IF(KKK.EQ.1) GO TO 405
  561.       DTEMP=DCOS(PI2*GRID(1))
  562.       DNUM=DCOS(PI2*GRID(NGRID))
  563.       AA=2.0/(DTEMP-DNUM)
  564.       BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
  565.   405 CONTINUE
  566.       DO 430 J=1,NFCNS
  567.       FT=J-1
  568.       FT=FT*DELF
  569.       XT=DCOS(PI2*FT)
  570.       IF(KKK.EQ.1) GO TO 410
  571.       XT=(XT-BB)/AA
  572.       XT1=SQRT(1.0-XT*XT)
  573.       FT=ATAN2(XT1,XT)/PI2
  574.   410 XE=X(L)
  575.       IF(XT.GT.XE) GO TO 420
  576.       IF((XE-XT).LT.FSH) GO TO 415
  577.       L=L+1
  578.       GO TO 410
  579.   415 A(J)=Y(L)
  580.       GO TO 425
  581.   420 IF((XT-XE).LT.FSH) GO TO 415
  582.       GRID(1)=FT
  583.       A(J)=GEED(1,NZ)
  584.   425 CONTINUE
  585.       IF(L.GT.1) L=L-1
  586.   430 CONTINUE
  587.       GRID(1)=GTEMP
  588.       DDEN=PI2/CN
  589.       DO 510 J=1,NFCNS
  590.       DTEMP=0.0
  591.       DNUM=J-1
  592.       DNUM=DNUM*DDEN
  593.       IF(NM1.LT.1) GO TO 505
  594.       DO 500 K=1,NM1
  595.       DAK=A(K+1)
  596.       DK=K
  597.   500 DTEMP=DTEMP+DAK*DCOS(DNUM*DK)
  598.   505 DTEMP=2.0*DTEMP+A(1)
  599.   510 ALPHA(J)=DTEMP
  600.       DO 550 J=2,NFCNS
  601.   550 ALPHA(J)=2.0*ALPHA(J)/CN
  602.       ALPHA(1)=ALPHA(1)/CN
  603.       IF(KKK.EQ.1) GO TO 545
  604.       P(1)=2.0*ALPHA(NFCNS)*BB+ALPHA(NM1)
  605.       P(2)=2.0*AA*ALPHA(NFCNS)
  606.       Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
  607.       DO 540 J=2,NM1
  608.       IF(J.LT.NM1) GO TO 515
  609.       AA=0.5*AA
  610.       BB=0.5*BB
  611.   515 CONTINUE
  612.       P(J+1)=0.0
  613.       DO 520 K=1,J
  614.       A(K)=P(K)
  615.   520 P(K)=2.0*BB*A(K)
  616.       P(2)=P(2)+A(1)*2.0*AA
  617.       JM1=J-1
  618.       DO 525 K=1,JM1
  619.   525 P(K)=P(K)+Q(K)+AA*A(K+1)
  620.       JP1=J+1
  621.       DO 530 K=3,JP1
  622.   530 P(K)=P(K)+AA*A(K-1)
  623.       IF(J.EQ.NM1) GO TO 540
  624.       DO 535 K=1,J
  625.   535 Q(K)=-A(K)
  626.       NF1J=NFCNS-1-J
  627.       Q(1)=Q(1)+ALPHA(NF1J)
  628.   540 CONTINUE
  629.       DO 543 J=1,NFCNS
  630.   543 ALPHA(J)=P(J)
  631.   545 CONTINUE
  632.       IF(NFCNS.GT.3) RETURN
  633.       ALPHA(NFCNS+1)=0.0
  634.       ALPHA(NFCNS+2)=0.0
  635.       RETURN
  636.       END
  637. C
  638. C-----------------------------------------------------------------------
  639. C FUNCTION: DD
  640. C   FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION
  641. C   COEFFICIENTS FOR USE IN THE FUNCTION GEED.
  642. C-----------------------------------------------------------------------
  643. C
  644.       DOUBLE PRECISION FUNCTION DD(K,N,M)
  645.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  646.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  647.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  648.       DIMENSION DES(1045),GRID(1045),WT(1045)
  649.       DOUBLE PRECISION AD,DEV,X,Y
  650.       DOUBLE PRECISION Q
  651.       DOUBLE PRECISION PI2
  652.       DD=1.0
  653.       Q=X(K)
  654.       DO 3 L=1,M
  655.       DO 2 J=L,N,M
  656.       IF(J-K)1,2,1
  657.     1 DD=2.0*DD*(Q-X(J))
  658.     2 CONTINUE
  659.     3 CONTINUE
  660.       DD=1.0/DD
  661.       RETURN
  662.       END
  663. C
  664. C-----------------------------------------------------------------------
  665. C FUNCTION: GEED
  666. C   FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE
  667. C   LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM
  668. C-----------------------------------------------------------------------
  669. C
  670.       DOUBLE PRECISION FUNCTION GEED(K,N)
  671.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  672.       COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
  673.       DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66)
  674.       DIMENSION DES(1045),GRID(1045),WT(1045)
  675.       DOUBLE PRECISION P,C,D,XF
  676.       DOUBLE PRECISION PI2
  677.       DOUBLE PRECISION AD,DEV,X,Y
  678.       P=0.0
  679.       XF=GRID(K)
  680.       XF=DCOS(PI2*XF)
  681.       D=0.0
  682.       DO 1 J=1,N
  683.       C=XF-X(J)
  684.       C=AD(J)/C
  685.       D=D+C
  686.     1 P=P+C*Y(J)
  687.       GEED=P/D
  688.       RETURN
  689.       END
  690. C
  691. C-----------------------------------------------------------------------
  692. C SUBROUTINE: OUCH
  693. C   WRITES AN ERROR MESSAGE WHEN THE ALGORITHM FAILS TO
  694. C   CONVERGE.  THERE SEEM TO BE TWO CONDITIONS UNDER WHICH
  695. C   THE ALGORITHM FAILS TO CONVERGE: (1) THE INITIAL
  696. C   GUESS FOR THE EXTREMAL FREQUENCIES IS SO POOR THAT
  697. C   THE EXCHANGE ITERATION CANNOT GET STARTED, OR
  698. C   (2) NEAR THE TERMINATION OF A CORRECT DESIGN,
  699. C   THE DEVIATION DECREASES DUE TO ROUNDING ERRORS
  700. C   AND THE PROGRAM STOPS.  IN THIS LATTER CASE THE
  701. C   FILTER DESIGN IS PROBABLY ACCEPTABLE, BUT SHOULD
  702. C   BE CHECKED BY COMPUTING A FREQUENCY RESPONSE.
  703. C-----------------------------------------------------------------------
  704. C
  705.       SUBROUTINE OUCH
  706.       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  707.       COMMON /OOPS/NITER
  708. C *      WRITE(*,1)NITER
  709. C *    1 FORMAT(44H ************ FAILURE TO CONVERGE **********/
  710. C *     141H0PROBABLE CAUSE IS MACHINE ROUNDING ERROR/
  711. C *     223H0NUMBER OF ITERATIONS =,I4/
  712. C *     339H0IF THE NUMBER OF ITERATIONS EXCEEDS 3,/
  713. C *     462H0THE DESIGN MAY BE CORRECT, BUT SHOULD BE VERIFIED WITH AN FFT)
  714.       RETURN
  715.       END
  716.