From:	AITGW::"NRL_ARCHIVE@nrlvax.nrl.navy.mil" 14-OCT-1991 10:41:07.58
To:	ARISIA::EVERHART
CC:	
Subj:	FIR.1

Received: by AITGW.DECnet (utk-mail11 v1.5) ; Mon, 14 Oct 91 10:40:08 EDT
Received: from nrlvx2.nrl.navy.mil by aitgw.ge.com (5.65/GE Gateway 1.4)
	id AA00268; Mon, 14 Oct 91 10:39:20 -0400
Message-Id: <9110141439.AA00268@aitgw.ge.com>
Received: by nrlvax.nrl.navy.mil (MX V2.3) id 16203; Mon, 14 Oct 1991 08:57:03
          EDT
Date: Mon, 14 Oct 1991 08:57:03 EDT
From: NRL_ARCHIVE-Mgr@nrlvax.nrl.navy.mil
Reply-To: NRL_ARCHIVE@nrlvax.nrl.navy.mil
Subject: FIR.1
To: ARISIA::EVERHART
 
X-NEWS: nrlvx1 comp.sources.misc: 682
Relay-Version: VMS News - V6.0-3 14/03/90 VAX/VMS V5.4; site nrlvx1.nrl.navy.mil
Path: nrlvx1.nrl.navy.mil!mimsy!haven.umd.edu!ames!elroy.jpl.nasa.gov!swrinde!cs.utexas.edu!uunet!sparky!kent
Newsgroups: comp.sources.misc
Subject: v22i100:  fir - FIR Linear Phase Filter Design, Part01/01
Message-ID: <1991Sep6.061742.15155@sparky.IMD.Sterling.COM>
From: sture@mdivax1.MDI.COM (John Sture)
Date: 6 Sep 91 06:17:42 GMT
Sender: kent@sparky.IMD.Sterling.COM (Kent Landfield)
Followup-To: comp.sources.d
Organization: Sterling Software, IMD
Approved: kent@sparky.imd.sterling.com
X-Md4-Signature: a66ad237eedc8d4137d8ef633abf6c99
Lines: 1636
 
Submitted-by: John Sture <sture@mdivax1.MDI.COM>
Posting-number: Volume 22, Issue 100
Archive-name: fir/part01
Environment: C++, DOS
 
[ Here's a new approach to portability, distribute 
[ it in multiple languages... :-)    
 
Fir is a C++ implementation of the Parks and McClellan remez exchange 
algorithm for design of fir filters.  Its limited to 52 coefficients 
but can easily be enhanced to extend the number of coefficients to the
published 128.
 
I have checked it against the examples in the book by Rabiner and Gold 
"Digital Filter Design" published by Prentice and Hall.
 
The fir4.cpp program was compiled using Borland's C++ V2.0 Compiler on 
a 386sx platform under DOS 4.01.  The commented lines in the code are 
the original FORTRAN code.
 
John Sture
---
#! /bin/sh
# This is a shell archive.  Remove anything before this line, then feed it
# into a shell via "sh file" or similar.  To overwrite existing files,
# type "sh file -c".
# The tool that generated this appeared in the comp.sources.unix newsgroup;
# send mail to comp-sources-unix@uunet.uu.net if you want that tool.
# Contents:  fir4.cpp
# Wrapped by kent@sparky on Fri Sep  6 01:07:12 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
echo If this archive is complete, you will see the following message:
echo '          "shar: End of archive."'
if test -f 'fir4.cpp' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'fir4.cpp'\"
else
  echo shar: Extracting \"'fir4.cpp'\" \(44237 characters\)
  sed "s/^X//" >'fir4.cpp' <<'END_OF_FILE'
X/*       PROGRAM FIR4
XC	REF: THEORY AND APPLICATION OF DIGITAL SIGNAL PROCESSING
XC	RABINER-GOLD 1975 PRENTICE-HALL PG.194-204
XC-----------------------------------------------------------------------
XC MAIN PROGRAM: FIR LINEAR PHASE FILTER DESIGN PROGRAM
XC               MODIFIED TO RUN ON INTERDATA FORTRAN VI COMPILER
XC
XC AUTHORS: JAMES H. MCCLELLAN
XC          DEPARTMENT OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE
XC          MASSACHUSETTS INSTITUTE OF TECHNOLOGY
XC          CAMBRIDGE, MASS. 02139
XC
XC          THOMAS W. PARKS
XC          DEPARTMENT OF ELECTRICAL ENGINEERING
XC          RICE UNIVERSITY
XC          HOUSTON, TEXAS 77001
XC
XC          LAWRENCE R. RABINER
XC          BELL LABORATORIES
XC          MURRAY HILL, NEW JERSEY 07974
XC
XC MODIFICATIONS: DECLARED PI AND PI2 EXPLICITLY.
XC                MODIFIED INPUT FORMAT.
XC                PETER R. H. MCCONNELL
XC
XC                80.7.17
XC                STRING VARIABLES ADDED.
XC                GENERAL REVISION -- DOUBLE PRECISION ADDED.
XC                80.7.23
XC                MODIFIED TO COMPILE UNDER PERKIN-ELMER FORTRAN 6.
XC
XC                PLOTTING OF IMPULSE, MAGNITUDE AND PHASE RESPONSES ADDED.
XC                ADDED CARD FOR INPUT OF PLOT TITLE (36 CHARACTERS).
XC                PLOTTING PERFORMED BY ROUTINE FIRPLOT
XC
XC         MODIFIED "EFF" AND "WATE" FOR INCLUSION OF EXTERNAL FUNCTIONS.
XC                ADDED ARGUEMENTS "EDGE" AND "NBANDS" TO "EFF" AND "WATE"
XC                FOR GREATER USER FLEXIBILITY.
XC                ADDED OPTION OF INPUT CARD (PER EACH FILTER DESIGN)
XC                FOR USE BY THE EXTERNAL "EFF" AND "WATE"" ROUTINES.
XC
XC                SPECIFICALLY, DESIGNED AN EXTERNAL "EFF" ROUTINE USING
XC                KAISER WINDOW WAITING IN THE PASSBAND OF LOWPASS FILTERS.
XC          ROUTINE USES OPTIONAL INPUT CARD FOR INPUT OF PARAMETER BETA.
XC                ROUTINE IN FILE KSREFF.OBJ.
XC
XC                80.7.31
XC                ADDED PERFORMANCE TESTING ROUTINE FOR L.P. FILTERS.
XC
XC                80.9.9
XC                UPDATED ROUTINES AND DOCUMENTATION FOR GENERAL USAGE.
XC
XC                81.10.14
XC                MODIFIED TO COMPILE UNDER VAX-11 FORTRAN.
XC
XC                GARY W KENWARD
XC
XC                MACDONALD, DETTWILER, AND ASSOCIATES,
XC                RICHMOND, B. C., CANADA, V6X 2Z9
XC
XC		 86.10.10
XC		 MODIFIED TO RUN UNDER CPM - MICROSOFT FORTRAN (F80)
XC		 CORRECTED BUG IN ARCOS CALCULATION AND MADE FRONT-END
XC		 SIMILAR TO ORIGINAL VERSION. TESTED AGAINST PUBLISHED
XC		 RESULTS AND WAS ACCURATE TO ABOUT 6 PLACES FOR FIG.S
XC		 7,9,11. ABOVE VERSION HAD PROBLEMS WITH HILBERT
XC		 TRANSFORMER.
XC		 JOHN STURE
XC
XC		 86.11.16
XC		 MADE COEFFICIENTS PRINTOUT AS ASCII INSTEAD OF FORTRAN
XC		 REAL VALUES ON DATA OUTPUT
XC		 ALSO LOOPS BACK TO GET MORE INPUT (FOR > 1 ENTRY)
XC
XC                91.08.14   JAS
XC		 Ported to C++.
XC		 91.08.24 JAS Finally got it working. Errors included
XC				.for loop incremented wrong variable
XC				.use of abs(integer) versus fabs(floating point)
XC
XC LOGICAL UNITS:
XC  5 -- INPUT
XC  6 -- OUTPUT
XC  7 -- DATA OUTPUT
XC 13 -- DEBUG OUTPUT FOR L.P. FILTER PERFORMANCE TESTING
XC
XC INPUT:
XC  NFILT -- FILTER LENGTH
XC
XC  JTYPE -- TYPE OF FILTER
XC          0 = LOWPASS FILTER (W/ MULTIPLE PASSBAND/STOPBANDS)
XC
XC IF PERFORMANCE TESTING IS DESIRED, JTYPE =0.
XC HOWEVER, THE FILTER HAS TO BE A LOW PASS TYPE (IE: THE UPPER MOST
XC BAND MUST BE THE ONE AND ONLY STOP BAND).
XC
XC          1 = MULTIPLE PASSBAND/STOPBAND FILTER
XC          2 = DIFFERENTIATOR
XC          3 = HILBERT TRANSFORM FILTER
XC
XC  NBANDS -- NUMBER OF BANDS
XC
XC  LGRID -- GRID DENSITY, WILL BE SET TO 16 UNLESS
XC           SPECIFIED OTHERWISE BY A POSITIVE CONSTANT.
XC
XC  JOUT -- IF JOUT = 1 OR JOUT >= 3:
XC                       -- IMPULSE AND MAGNITUDE RESPONSES
XC                          WILL BE PLOTTED ON UNIT 4.
XC                       -- THE IMPULSE RESPONSE OF THE FILTER WILL BE
XC                          OUTPUTTED TO UNIT "COEFF".
XC
XC  EDGE(2*NBANDS) -- BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND
XC                    WITH A MAXIMUM OF 10 BANDS.
XC
XC  FX(NBANDS) -- DESIRED FUNCTION ARRAY (OR DESIRED SLOPE IF A
XC                DIFFERENTIATOR) FOR EACH BAND.
XC
XC  WTX(NBANDS) -- WEIGHT FUNCTION ARRAY IN EACH BAND.  FOR A
XC                 DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY
XC                 PROPORTIONAL TO F.
XC
XC  IF PLOTTING IS TO BE DONE (JOUT = 1 OR >= 3), THE NEXT CARD
XC  IS THE TITLE OF
XC  PLOT, OR IF ONLY ONE FILTER IS BEING DESIGNED, AN EOF WILL SET A BLANK
XC  TITLE (NOTE THAT THIS PRECLUDES THE USE OF THE NEXT INPUT CARD).
XC  THE TITLE CAN BE UP TO 36 CHARACTERS LONG, AND IS REPEATED ON EACH PLOT.
XC
XC  IF AN EXTERNAL "EFF" AND/OR "WATE" IS INCLUDED, THE NEXT
XC  CARD CAN BE USED FOR
XC  UP TO 10 USER COMMON PARAMETERS.
XC
XC  SAMPLE INPUT DATA SETUP:
XC
XC  NFILT,JTYPE,NBANDS,LGRID,JOUT                  5I3
XC  EDGE1,EDGE2,FX,WTX                             4F15.9
XC
XC   32  1  3  0  0
XC       0.0       0.10          0.0           10.0
XC       0.20       0.35          1.0           1.0
XC       0.425      0.5           0.0          10.0
XC
XC     THIS DATA SPECIFIES A LENGTH 32 BANDPASS FILTER WITH
XC     STOPBANDS 0 TO 0.1 AND 0.425 TO 0.5, AND PASSBAND FROM
XC     0.2 TO 0.35 WITH WEIGHTING OF 10 IN THE STOPBANDS AND 1
XC     IN THE PASSBAND.  THE GRID DENSITY DEFAULTS TO 16.
XC     THIS IS THE FILTER IN FIGURE 10.
XC
XC     THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 FULLBAND
XC     DIFFERENTIATOR WITH SLOPE 1 AND WEIGHTING OF 1/F.
XC     THE GRID DENSITY WILL BE SET TO 20.
XC   32  2  1  20  0
XC       0.0        0.5           1.0          1.0
XC
XC-----------------------------------------------------------------------
XC
XC      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
XC*/
X
X#include <stdlib.h>
X//include <string.h>
X#include <iostream.h>
X#include <math.h>
X
X/*    EXTERNAL WATE
X      EXTERNAL EFF   */
X
X// external functions
Xdouble wate ( double*, double*, double*, double*, int, int, int);
Xdouble eff  ( double*, double*, double*, double*, int, int, int);
X
X// local functions
Xdouble d    ( int, int, int);
Xvoid   error( int );
Xvoid   remez( double* , int);
Xvoid   ouch ( void);
Xdouble gee  ( int, int);
X
X// local definitions
X#define false 0;
X#define true 1;
X
X
X
X//  Variables
Xchar rdarg;
Xdouble ad[28+2], alpha[28+2], x[28+2], y[28+2];
X//      DOUBLE PRECISION AD(28),ALPHA(28),X(28),Y(28)
Xdouble des[488+2],grid[488+2], wt[488+2];
X//      DOUBLE PRECISION DES(488),GRID(488),WT(488)
Xdouble edge[20+2], fx[10+2], wtx[10+2],eargs[10+2];
X//      DOUBLE PRECISION EDGE(20),FX(10),WTX(10),EARGS(10)
Xdouble h[18+2], deviat[10+2], devtmp, hout[128+2];
X//      DOUBLE PRECISION H(18),DEVIAT(10),DEVTMP,HOUT(128)
Xdouble pi, pi2, dev, delf, fup, temp, change, xt, dtemp, dnum, dden;
X//      DOUBLE PRECISION PI,PI2,DEV,DELF,FUP,TEMP,CHANGE,XT
Xint nfmax, jtype, nfilt, jout, lgrid, iext[28+2], lptype, jxt;
X//      INTEGER BD1,BD2,BD3,BD4,IEXT(28)
Xint input, output, coeff, nbands, jb, j, k, neg, nodd, nfcns, l;
Xint lband, ngrid, nm1, nz, nzmj, nf2j, nf3j, kup, ix, i, niter;
X//      INTEGER INPUT,OUTPUT,COEFF,NBANDS
X//      INTEGER*1 TIT1(30),TIT2(36),TIT3(26),STIT0(15)
X//      INTEGER*1 STIT1(16),STIT2(14)
X//      INTEGER*1 STIT3(20),TITLI(36)
X// C
X//      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
X//      COMMON /OOPS/NITER,OUTPUT
X//      COMMON /EFFARG/RDARG,EARGS,INPUT
X//      DATA BD1,BD2,BD3,BD4/'b','a','n','d'/
X//      DATA TIT1/  'F','i','n','i','t','e',' ',
X//     -'I','m','p','u','l','s','e',
X//     -' ','R','e','s','p','o','n',
X//     -'s','e',' ','(','F','I','R',')',' '/
X//      DATA TIT2/  'L','i','n','e','a','r',' ',
X//     -'P','h','a','s','e',' ','D',
X//     -'i','g','i','t','a','l',' ','F','i','l',
X//     -'t','e','r',' ','D','e','s',
X//     -'i','g','n',' ',' '/
X//      DATA TIT3/  'R','e','m','e','z',' ','E',
X//     -'x','c','h','a','n','g',
X//     -'e',' ','A','l','g','o','r','i','t','h',
X//     -'m',' ',' '/
X//      DATA STIT0/ 'L','o','w','p','a','s','s',
X//     -' ','f','i','l','t','e',
X//     -'r',' '/
X//      DATA STIT1/ 'B','a','n','d','p','a','s',
X//     -'s',' ','f','i','l','t',
X//     -'e','r',' '/
X//      DATA STIT2/ 'D','i','f','f','e','r','e',
X//     -'n','t','i','a','t','o','r'/
X//      DATA STIT3/ 'H','i','l','b','e','r','t',
X//     -' ','t','r','a','n','s',
X//     -'f','o','r','m','e','r',' '/
X//C
X
X
X
Xmain()
X{
X//      INPUT=6
X//      OUTPUT=7
X//      COEFF=8
X//C
X//      CALL OPEN(INPUT, 'FIRIN   DAT',0)
X//      WRITE(1,2222)
X//2222  FORMAT(' INPUT IN FIRIN.DAT, OUTPUT IN FORT07.DAT, COEFF IN ',
X//     %'FORT08.DAT')
X//C
X      pi  = 3.141592653589793;
X      pi2 = pi*2.0;
X//C
X//      RDARG = .FALSE.
X      rdarg   = false;
X//      NTITL = 36
X//C
X//C  THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128, BUT
X//C  THIS UPPER LIMIT CAN BE CHANGED BY REDIMENSIONING THE
X//C  ARRAYS IEXT, AD, ALPHA, X, Y, H TO BE NFMAX/2 + 2.
X//C  THE ARRAYS DES, GRID, AND WT MUST BE DIMENSIONED
X//C  16(NFMAX/2 + 2).
X//C
X//      NFMAX=52
X      nfmax=52;
Xl100:
X//100   CONTINUE
X//      JTYPE=0
X      jtype=0;
X//      JOUT=0
X      jout=0;
X//C
X//C  PROGRAM INPUT SECTION
X//C
X//      READ(INPUT,110,END=470,ERR=470) NFILT,JTYPE,NBANDS,LGRID,JOUT
X//110   FORMAT(5I3)
X      cout<<"Enter nfilt, jtype, nbands, lgrid, jout\n";
X      cin >> nfilt;  cout<< nfilt  <<"\n";
X      cin >> jtype;  cout<< jtype  <<"\n";
X      cin >> nbands; cout<< nbands <<"\n";
X      cin >> lgrid;  cout<< lgrid  <<"\n";
X      cin >> jout;   cout<< jout   <<"\n";
X
X//      WRITE (OUTPUT,111) NFILT,JTYPE,NBANDS,LGRID,JOUT
X//111   FORMAT (5I3)
X//      WRITE (1,3111) NFILT,JTYPE,NBANDS,LGRID,JOUT
X//3111  FORMAT (5I3)
X//C
X//C     READ(INPUT,*,END=470) NFILT,JTYPE,NBANDS,LGRID,JOUT
X//C
X//      LPTYPE = JTYPE
X      lptype = jtype;
X//      IF (JTYPE.EQ.0) JTYPE = 1
X      if (jtype == 0) error(2);
X//      IF(NFILT.GT.NFMAX.OR.NFILT.LT.3) CALL ERROR
X      if ((nfilt > nfmax) || (nfilt < 3)) error(-1);
X//      IF(NBANDS.LE.0) NBANDS=1
X      if (nbands <= 0) nbands = 1;
X//C
X//C  GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED
X//C  OTHERWISE
X//C
X//      IF(LGRID.LE.0) LGRID=16
X      if (lgrid <= 0 ) lgrid = 16;
X//      JB=2*NBANDS
X      jb=2*nbands;
X//      READ(INPUT,2120) (EDGE(J),J=1,JB)
X      cout<<"Enter "<<jb<<" edges\n";
X      for (j=1; j<=jb; j++) { cin >> edge[j]; cout <<" "<<edge[j]<<"\n";}
X//2120  FORMAT(8F15.9)
X//      WRITE(1,3120) (EDGE(J),J=1,JB)
X      for (j=1; j<=jb; j++) { cout << " edge["<<j<<"]= "<<edge[j]<<"\n";}
X//3120  FORMAT(8F15.9)
X//      READ(INPUT,2121) (FX(J),J=1,NBANDS)
X      cout<<"Enter "<<nbands<<" fx's \n";
X      for (j=1; j<=nbands; j++) {cin >> fx[j]; cout <<" "<<fx[j]<<"\n";}
X//2121  FORMAT(F15.9)
X//      WRITE(1,3121) (FX[J),J=1,NBANDS)
X//3121  FORMAT(8F15.9)
X      for (j=1; j<=nbands; j++) {cout<<" fx["<<j<<"]= "<<fx[j]<<"\n";}
X//      READ(INPUT,2122) (WTX[J),J=1,NBANDS)
X//2122  FORMAT(8F15.9)
X      cout<<"Enter "<<nbands<<" wtx's\n";
X      for (j=1; j<=nbands; j++) {cin >> wtx[j]; cout<<" "<<wtx[j]<<"\n";}
X//      WRITE(1,3122) (WTX[J),J=1,NBANDS)
X//3122  FORMAT(8F15.9)
X      for (j=1; j<=nbands; j++) {cout<<" wtx["<<j<<"[= "<<wtx[j]<<"\n";}
X//C
X//C     MODIFIED READ
X//C
X//C     DO 120 J=1,NBANDS
X//C        J1=2*(J-1)+1
X//C        J2=J1+1
X//C
X//C          READ(INPUT,121,END=470)EDGE(J1),EDGE(J2),FX(J),WTX(J)
X//C121       FORMAT(4F15.9)
X//C	  WRITE(OUTPUT,112)J1,J,EDGE(J1),EDGE(J2),FX(J),WTX(J)
X//C112	  FORMAT(2I3,4F15.9)
X//C
X//C         READ(INPUT,*,END=470)EDGE(J1),EDGE(J2),FX(J),WTX(J)
X//120   CONTINUE
Xl120:
X//      IF(JTYPE.LE.0.OR.JTYPE.GT.3) CALL ERROR
X      if ((jtype<=0) || (jtype > 3)) error(-2);
X//C
X//C IF PLOTTING IS TO BE DONE, READ IN A TITLE FOR THE PLOTS.
X//C
X//C     IF((JOUT.EQ.1).OR.(JOUT.GE.3))
X//
X//C    - READ(INPUT,1120,END=125,ERR=125) TITLI
X//C1120  FORMAT(36A1)
X//C
X//      WRITE(OUTPUT,124)
X//124   FORMAT('1')
X//C
X//      GOTO 126
X      goto l126;
X//125   CONTINUE
Xl125:
X//      CALL ERROR
X      error(-3);
X//C
X//126   CONTINUE
Xl126:
X//      NEG=1
X      neg=1;
X//      IF(JTYPE.EQ.1) NEG=0
X      if (jtype==1) neg=0;
X//      NODD=NFILT/2
X      nodd=nfilt/2;
X//      NODD=NFILT-2*NODD
X      nodd=nfilt-2*nodd;
X//      NFCNS=NFILT/2
X      nfcns=nfilt/2;
X//      IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS=NFCNS+1
X      if((nodd==1) && (neg == 0)) nfcns = nfcns + 1;
X//C
X//C  SET UP THE DENSE GRID.  THE NUMBER OF POINTS IN THE GRID
X//C  IS (FILTER LENGTH + 1)*GRID DENSITY/2
X//C
X//      GRID(1)=EDGE(1)
X      grid[1] = edge[1];
X//      DELF=LGRID*NFCNS
X      delf = (lgrid*nfcns);
X//      DELF=0.5D0/DELF
X	delf = 0.5/delf;
X//      IF(NEG.EQ.0) GO TO 135
X      if (neg==0) goto l135;
X//      IF(EDGE(1).LT.DELF) GRID(1)=DELF
X      if (edge[1]<delf) grid[1]=delf;
X//135   CONTINUE
Xl135:
X//      J=1
X      j=1;
X//      L=1
X      l=1;
X//      LBAND=1
X      lband=1;
X//140   FUP=EDGE[L+1)
Xl140:
X      fup = edge[l+1]; //l+1-1
X//145   TEMP=GRID(J)
Xl145:
X      temp = grid[j];
X//C
X//C  CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
X//C  FUNCTION ON THE GRID
X//C
X//      DES(J)=EFF(TEMP,FX,WTX,EDGE,LBAND,JTYPE,NBANDS)
X     des[j] = eff(&temp, fx, wtx, edge, lband, jtype, nbands);
X//     cout<<"des["<<j<<"]="<<des[j]<<"\n";
X//      WT(J)=WATE(TEMP,FX,WTX,EDGE,LBAND,JTYPE,NBANDS)
X     wt[j] = wate(&temp, fx, wtx, edge, lband, jtype, nbands);
X//     cout<<"wt["<<j<<"]="<<wt[j]<<"\n";
X//      J=J+1
X     j=j+1;
X//      GRID(J)=TEMP+DELF
X     grid[j] = temp + delf;
X//      IF(GRID(J).GT.FUP) GO TO 150
X     if (grid[j]>fup) goto l150;
X//      GO TO 145
X     goto l145;
X//150   GRID(J-1)=FUP
Xl150:
X     grid[j-1] = fup;
X//      DES(J-1)=EFF(FUP,FX,WTX,EDGE,LBAND,JTYPE,NBANDS)
X     des[j-1] = eff(&fup, fx, wtx, edge, lband, jtype, nbands);
X//      WT(J-1)=WATE(FUP,FX,WTX,EDGE,LBAND,JTYPE,NBANDS)
X     wt[j-1] = wate(&fup, fx, wtx, edge, lband, jtype, nbands);
X//      LBAND=LBAND+1
X     lband=lband+1;
X//      L=L+2
X     l=l+2;
X//      IF(LBAND.GT.NBANDS) GO TO 160
X     if (lband>nbands) goto l160;
X//      GRID(J)=EDGE(L)
X     grid[j] = edge[l];
X//      GO TO 140
X     goto l140;
X//160   NGRID=J-1
Xl160:
X     ngrid=j-1;
X//      IF(NEG.NE.NODD) GO TO 165
X     if (neg!=nodd) goto l165;
X//      IF(GRID(NGRID).GT.(0.5D0-DELF)) NGRID=NGRID-1
X     if (grid[ngrid]>(0.5-delf)) ngrid=ngrid-1;
X//165   CONTINUE
Xl165:
X//C
X//C  SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
X//C  TO THE ORIGINAL PROBLEM
X//C
X//      IF(NEG) 170,170,180
X     if (neg <= 0) goto l170;
X     goto l180;
X//170   IF(NODD.EQ.1) GO TO 200
Xl170:
X     if(nodd==1) goto l200;
X//      DO 175 J=1,NGRID
X     for (j=1; j<=ngrid; j++) {
X//      CHANGE=DCOS(PI*GRID(J))
X	     change=cos(pi*grid[j]);
X//      DES(J)=DES(J)/CHANGE
X	     des[j]=des[j]/change;
X//175   WT(J)=WT(J)*CHANGE
Xl175:
X	     wt[j] = wt[j] * change;}
X//      GO TO 200
X     goto l200;
X//180   IF(NODD.EQ.1) GO TO 190
Xl180:
X      if (nodd==1) goto l190;
X//      DO 185 J=1,NGRID
X      for (j=1; j<=ngrid; j++) {
X//      CHANGE=DSIN(PI*GRID(J))
X	      change=sin(pi*grid[j]);
X//      DES(J)=DES(J)/CHANGE
X	      des[j]=des[j]/change;
X//185   WT(J)=WT(J)*CHANGE
X	      wt[j]=wt[j]*change;}
X//      GO TO 200
X      goto l200;
X//190   DO 195 J=1,NGRID
Xl190:
X      for (j=1; j<=ngrid; j++) {
X//      CHANGE=DSIN(PI2*GRID(J))
X	      change=sin(pi2*grid[j]);
X//      DES(J)=DES(J)/CHANGE
X	      des[j]=des[j]/change;
X//195   WT(J)=WT(J)*CHANGE
X	      wt[j]=wt[j]*change; }
X//C
X//C  INITIAL GUESS FOR THE EXTREMAL FREQUENCIES--EQUALLY
X//C  SPACED ALONG THE GRID
X//C
X// 200   TEMP=DBLE(double(NGRID-1))/DBLE(double(NFCNS))
Xl200:
X      temp=(ngrid-1);
X      temp = temp/(nfcns);
X//      DO 210 J=1,NFCNS
X      for(j=1; j<=nfcns; j++) {
X//      XT=J-1
X	      xt=j-1;
X//210   IEXT(J)=XT*TEMP+1
X	      iext[j]=(xt*temp) + 1; }
X//      IEXT(NFCNS+1)=NGRID
X      iext[nfcns+1]=ngrid;   // subtract 1 for test >>>jas 18-aug-1991
X//     NM1=NFCNS-1
X      nm1 =nfcns-1;
X//      NZ=NFCNS+1
X      nz  =nfcns+1;
X//C
X//C  CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION
X//C  PROBLEM
X//C
X//    CALL REMEZ(EDGE,NBANDS)
X      remez(edge, nbands);
X//C
X//C  CALCULATE THE IMPULSE RESPONSE.
X//C
X//      IF(NEG) 300,300,320
X      if (neg<=0) goto l300;
X      goto l320;
X//300   IF(NODD.EQ.0) GO TO 310
Xl300:
X      if (nodd==0) goto l310;
X//      DO 305 J=1,NM1
X      for(j=1; j<=nm1; j++) {
X//      NZMJ=NZ-J
X	      nzmj = nz-j;
X//305   H(J)=0.5D0*ALPHA(NZMJ)
X	      h[j]=0.5 * alpha[nzmj];}
X//      H(NFCNS)=ALPHA(1)
X      h[nfcns]=alpha[1];
X//      GO TO 350
X      goto l350;
X//310   H(1)=0.25D0*ALPHA(NFCNS)
Xl310:
X      h[1]=0.25*alpha[nfcns];
X//      DO 315 J=2,NM1
X      for (j=2; j<=nm1; j++) {
X//      NZMJ=NZ-J
X	      nzmj=nz-j;
X//      NF2J=NFCNS+2-J
X	      nf2j=nfcns+2-j;
X//315   H(J)=0.25D0*(ALPHA(NZMJ)+ALPHA(NF2J))
X	      h[j]=0.25*(alpha[nzmj]+alpha[nf2j]); }
X//      H(NFCNS)=0.5D0*ALPHA(1)+0.25D0*ALPHA(2)
X      h[nfcns]=0.5*alpha[1]+0.25*alpha[2];
X//      GO TO 350
X      goto l350;
X//320   IF(NODD.EQ.0) GO TO 330
Xl320:
X      if(nodd==0) goto l330;
X//      H(1)=0.25D0*ALPHA(NFCNS)
X      h[1]=0.25*alpha[nfcns];
X//      H(2)=0.25D0*ALPHA(NM1)
X      h[2]=0.25*alpha[nm1];
X//      DO 325 J=3,NM1
X      for (j=3; j<=nm1; j++) {
X//325   H(J]=0.25D0*(ALPHA(NZ-J)-ALPHA(NFCNS-3+J))
X	      h[j]=0.25*(alpha[nz-j]-alpha[nfcns+3-j]);  }
X//      H(NFCNS)=0.5D0*ALPHA(1)-0.25D0*ALPHA(3)
X      h[nfcns]=0.5*alpha[1]-0.25*alpha[3];
X//      H(NZ)=0.0
X      h[nz]=0.0;
X//      GO TO 350
X      goto l350;
X//330   H(1)=0.25D0*ALPHA(NFCNS)
Xl330:
X      h[1]=0.25*alpha[nfcns];
X//      DO 335 J=2,NM1
X      for(j=2; j<=nm1; j++) {
X//335   H(J)=0.25D0*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))
X	      h[j]=0.25 * (alpha[nz-j]-alpha[nfcns+2-j]); }
X//      H(NFCNS)=0.5D0*ALPHA(1)-0.25D0*ALPHA(2)
X      h[nfcns]=0.5*alpha[1]-0.25*alpha[2];
X//C
X//C  PROGRAM OUTPUT SECTION.
X//C
X//350   WRITE(OUTPUT,360) TIT1, TIT2, TIT3
Xl350:
X//360   FORMAT('1',70('*')//15X,30A1,/13X,36A1,/17X,26A1)
X//      IF (LPTYPE.EQ.0)
X      if (lptype==1) cout << "Bandpass filter \n";
X//     - WRITE(OUTPUT,364) STIT0
X//364   FORMAT(21X,15A1)
X//      IF(LPTYPE.EQ.1)
X      if (lptype==2) cout << "Differentiator\n";
X//     - WRITE(OUTPUT,365) STIT1
X//365   FORMAT(21X,16A1)
X//      IF(LPTYPE.EQ.2)
X//     - WRITE(OUTPUT,370) STIT2
X//370   FORMAT(21X,14A1)
X//      IF(LPTYPE.EQ.3)
X      if (lptype==3) cout << "Hilbert Transformer \n";
X//     - WRITE(OUTPUT,375) STIT3
X//375   FORMAT(20X,20A1)
X//      IF ((JOUT.EQ.1).OR.(JOUT.GE.3)) WRITE(OUTPUT,377) TITLI
X//377   FORMAT('0',15X,36A1)
X//      WRITE(OUTPUT,380)
X//380   FORMAT('0',18X,'Impulse Response',/' ',18X,'------- --------')
X//      WRITE(OUTPUT,378) NFILT
X//378   FORMAT('0',17X,'Filter length = ',I3/' ')
X      cout << "\n\n------Impulse Response------\n ";
X      cout << "Filter Length = " << nfilt << "\n";
X//      DO 381 J=1,NFCNS
X      for (j=1; j<=nfcns; j++) {
X//      K=NFILT+1-J
X	      k=nfilt+1-j;
X//      IF(NEG.EQ.0) WRITE(OUTPUT,382) J,H(J),K
X	      if (neg==0) cout <<"h("<<j<<") =    "<<h[j]<<"    =h("<<k<<")\n";
X//      IF(NEG.EQ.1) WRITE(OUTPUT,383) J,H(J),K
X	      if (neg==1) cout <<"h("<<j<<") =    "<<h[j]<<"    =-h("<<k<<")\n";
X//      IF (MOD(J,45).EQ.0) WRITE(OUTPUT,3381)
X//3381  FORMAT('1')
X//381   CONTINUE
X	}
X//382   FORMAT(13X,'h(',I2,') = ',E15.8,' = h(',I3,')')
X//383   FORMAT(13X,'h(',I2,') = ',E15.8,' = -h(',I3,')')
X//      IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(OUTPUT,384) NZ
X//384   FORMAT(13X,'h(',I2,') =  0.0')
X      if ((neg==1)&&(nodd==1))  cout <<"h("<<nz<<") = 0.0\n";
X//      DO 450 K=1,NBANDS,4
X      for (k=1; k<=nbands; k +=4) {
X//      KUP=K+3
X	      kup=k+3;
X//      IF(KUP.GT.NBANDS) KUP=NBANDS
X	      if (kup>nbands) kup=nbands;
X//      WRITE(OUTPUT,385) (BD1,BD2,BD3,BD4,J,J=K,KUP)
X//385   FORMAT(/24X,4(4A1,I3,7X))
X	      cout <<"\n band ";
X	      for (j=k; j<=kup; j++) { cout << j<< "      ";}
X	      cout <<"\n";
X//      WRITE(OUTPUT,390) (EDGE(2*J-1),J=K,KUP)
X//390   FORMAT(2X,'Lower band edge',5F14.7)
X	      cout<<"\nLower band edge:\n";
X	      for (j=k;j<=kup; j++){ cout<<edge[2*j-1]<<" ";}
X	      cout<<"\n";
X//      WRITE(OUTPUT,395) (EDGE(2*J),J=K,KUP)
X//395   FORMAT(2X,'Upper band edge',5F14.7)
X	      cout<<"\nUpper band edge:\n";
X	      for (j=k; j<=kup; j++) { cout<<edge[2*j]<<" ";}
X	      cout<<"\n";
X//      IF(LPTYPE.NE.2) WRITE(OUTPUT,400) (FX(J),J=K,KUP)
X//400   FORMAT(2X,'Desired value',2X,5F14.7)
X	      if (lptype!=2) {cout<<"\nDesired values:";
X			for (j=k;j<=kup;j++){cout<<fx[j]<<"     ";}
X			cout <<"\n";
X	      }
X//     IF(LPTYPE.EQ.2) WRITE(OUTPUT,405) (FX(J),J=K,KUP)
X//405   FORMAT(2X,'Desired slope',2X,5F14.7)
X	      if (lptype==2) {cout<<"\nDesired slope:";
X			for (j=k;j<=kup;j++){cout<<fx[j]<<"     ";}
X			cout <<"\n";
X	      }
X//      WRITE(OUTPUT,410) (WTX(J),J=K,KUP)
X//410   FORMAT(2X,'Weighting',6X,5F14.7)
X	       cout<<"\nWeighting:";
X	       for (j=k; j<=kup; j++) {cout<<wtx[j]<<"     ";}
X	       cout<<"\n";
X//    DO 420 J=K,KUP
X	       for (j=k; j<=kup;j++) {
X//420   DEVIAT(J)=DEV/WTX(J)
X			deviat[j]= dev/wtx[j]; }
X//      WRITE(OUTPUT,425) (DEVIAT(J),J=K,KUP)
X//425   FORMAT(2X,'Deviation',6X,5F14.7)
X	       cout<<"\nDeviation:";
X	       for (j=k; j<=kup; j++) {cout<<deviat[j]<<"     ";}
X	       cout<<"\n";
X//      IF(JTYPE.NE.1) GO TO 450
X	       if (jtype!=1) goto l450;
X//      DO 430 J=K,KUP
X	       for (j=k;j<=kup;j++) {
X//      DEVTMP = DEVIAT[J]
X		      devtmp=deviat[j];
X//      DEVIAT(J)= 20.0 * DLOG10(DEVTMP)
X		      deviat[j]=20.0 * log(devtmp)/log(10.0);
X//430   CONTINUE
X	       }
X//      WRITE(OUTPUT,435) (DEVIAT(J),J=K,KUP)
X//435   FORMAT(2X,'Deviation in dB',5F14.7)
X	       cout<<"\nDeviation in dB:";
X	       for (j=k;j<=kup;j++) {cout<<deviat[j]<<"    ";}
X	       cout<<"\n";
X//450   CONTINUE
Xl450:	;}
X//      DO 452 J=1,NZ
X	for (j=1;j<=nz;j++){
X//      IX=IEXT(J)
X		ix=iext[j];
X//452   GRID(J)=GRID(IX)
X		grid[j]=grid[ix];}
X//      WRITE(OUTPUT,455) (GRID(J),J=1,NZ)
X//455   FORMAT(/2X,'Extremal frequencies -- maxima of the error curve:'/
X//     1 (2X,5F12.7))
X	cout<<"\nExtremal frequencies -- maxima of the error curve:\n";
X	for (j=1;j<=nz;j+=2){cout<<grid[j]<<"      "<<grid[j+1]<<"\n";}
X	cout<<"\n";
X//      WRITE(OUTPUT,460)
X//460   FORMAT(/1X,70('*'))
X//C      IF ((JOUT.GE.2).AND.(LPTYPE.EQ.0))
X//C     -  CALL PERFOR(H,NFCNS,NFILT,NODD,EDGE,NBANDS)
X//C      IF ((JOUT.EQ.1).OR.(JOUT.GE.3))
X//C     -  CALL FIRPLOT(NTITL,TITLI,H,NFILT,NFCNS,NODD,NEG)
X//      IF ((JOUT.NE.1).AND.(JOUT.LT.3)) GOTO 470
X//       if ((jout!=1)&&(jout<3)) goto l470;
X//	DO 461 J = 1, NFCNS
X       for (j=1; j<=nfcns; j++){
X//          I = NFILT - J + 1
X		i=nfilt-j+1;
X//          HOUT(J) = H(J)
X		hout[j]=h[j];
X//          HOUT(I) = H(J)
X		hout[i]=h[j];}
X//461     CONTINUE
X//        IF ((NODD.EQ.1).AND.(NEG.EQ.1)) HOUT(NFCNS+1) = 0.0
X	if ((nodd==1)&&(neg==1)) hout[nfcns+1]=0.0;
X//	WRITE(COEFF,471)NFILT,(HOUT(J), J=1,NFILT)
X//471	FORMAT(' ',I4,/,128(F14.11,/))
X	cout<<"\n\nFilter coefficients:\n";
X	cout<<"\nnumber of coefficients:"<<nfilt<<"\n";
X	for (j=1;j<=nfilt; j++){cout<<hout[j]<<"\n";}
X//      GO TO 100
X//	goto l100;
X//470   CONTINUE
Xl470: ;
X//      STOP
X	exit(0);
X//      END
X}
X//C-----------------------------------------------------------------------
X//C
X//C 80.7.23
X//C FUNCTION "EFF" IS NOW IN FILE "EFF.FTN" TO FACILITATE USE OF ALTERNATE/
X//C FUNCTIONS.
X//C GARY W. KENWARD
X//C jas aug-15-1991 returned here for simplicity
X//C-----------------------------------------------------------------------
X//C FUNCTION: EFF
X//C   FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE
X//C   AS A FUNCTION OF FREQUENCY.
X//C   AN ARBITRARY FUNCTION OF FREQUENCY CAN BE
X//C   APPROXIMATED IF THE USER REPLACES THIS FUNCTION
X//C   WITH THE APPROPRIATE CODE TO EVALUATE THE IDEAL
X//C   MAGNITUDE.  NOTE THAT THE PARAMETER FREQ IS THE
X//C   VALUE OF NORMALIZED FREQUENCY NEEDED FOR EVALUATION.
X//C
X//C 81.10.14
X//C MODIFIED TO COMPILE UNDER VAX-11 FORTRAN
X//C
X//C GARY KENWARD
X//C-----------------------------------------------------------------------
X//C
X//      DOUBLE PRECISION FUNCTION EFF(FREQ,FX,WTX,EDGE,LBAND,JTYPE,NBANDS)
X//C
Xdouble eff ( double* freq,double* fx,double* wtx,double* edge, int lband,
X		 int jtype, int nbands)
X{
X//      DOUBLE PRECISION FX(10),WTX(10),EDGE(20),FREQ
Xdouble ret_eff;
X//
X//      IF(JTYPE.EQ.2) GO TO 100
X      if(jtype==2) goto lll100;
X//      EFF=FX(LBAND)
X      ret_eff=fx[lband];
X//      RETURN
X      return(ret_eff);
X// 100  EFF=FX(LBAND)*FREQ
Xlll100:
X      ret_eff=fx[lband]* (*freq);
X//      RETURN
X     return(ret_eff);
X//      END
X}
X//C---------------------------------------------------------------------
X//C---------------------------------------------------------------------
X//C
X//C 80.9.9
X//C FUNCTION "WATE" IS NOW IN FILE "WATE.FTN" TO FACILITATE USE OF ALTERNATE
X//C FUNCTIONS.
X//C GARY W. KENWARD
X//C aug-15-1991 jas . returned here.
X//C---------------------------------------------------------------------
X//C FUNCTION: WATE
X//C   FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION
X//C   OF FREQUENCY.  SIMILAR TO THE FUNCTION EFF, THIS FUNCTION CAN
X//C   BE REPLACED BY A USER-WRITTEN ROUTINE TO CALCULATE ANY
X//C   DESIRED WEIGHTING FUNCTION.
X//C
X//C 80.7.23
X//C ARGUEMENTS "EDGE" AND "NBANDS" ADDED TO GIVE MORE FLEXIBILITY TO USER
X//C FUNCTIONS.
X//C
X//C 80.9.9
X//C ROUTINE MADE EXTERNAL
X//C
X//C 81.10.14
X//C MODIFIED TO COMPILE UNDER VAX-11 FORTRAN
X//C
X//C GARY W. KENWARD
X//C-----------------------------------------------------------------------
X//C
X//      DOUBLE PRECISION FUNCTION WATE(F,FX,WTX,EDGE,LBAND,JTYPE,NBANDS)
Xdouble wate(double* f, double* fx, double* wtx, double* edge, int lband,
X			int jtype, int nbands)
X{
X//C     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
X//      DOUBLE PRECISION FX(10),WTX(10),EDGE(20),F
X//      IF(JTYPE.EQ.2) GO TO 1
Xdouble ret_wate;
X      if(jtype==2) goto v1;
X//      WATE=WTX(LBAND)
X      ret_wate=wtx[lband];
X//      RETURN
X      return(ret_wate);
X//    1 IF(FX(LBAND).LT.0.0001D0) GO TO 2
Xv1:
X      if (fx[lband]<0.0001) goto v2;
X//      WATE=WTX(LBAND)/F
X      ret_wate=wtx[lband]/(*f);
X//      RETURN
X      return(ret_wate);
X//    2 WATE=WTX(LBAND)
Xv2:
X      ret_wate=wtx[lband];
X//      RETURN
X      return(ret_wate);
X//      END
X}
X//C
X//C-----------------------------------------------------------------------
X//C SUBROUTINE: ERROR
X//C /  THIS ROUTINE WRITES AN ERROR MESSAGE IF AN
X//C   ERROR HAS BEEN DETECTED IN THE INPUT DATA.
X//C   80.7.23
X//C   IT STOPS THE PROGRAM BY CALLING EXIT, RETURNING A CODE OF 1.
X//C   GARY W. KENWARD
X//C /
X//C-----------------------------------------------------------------------
X//C
X//      SUBROUTINE ERROR
Xvoid error(int err)
X{
X//      INTEGER NITER,OUTPUT
X//      COMMON /OOPS/NITER,OUTPUT
X//      WRITE(OUTPUT,1)
X//1     FORMAT(' ************ ERROR IN INPUT DATA **********')
X//      CALL EXIT(1)
X//      RETURN
X//      END
X	cout<<" ****** error in input data  ************"<<err<"\n";
X	exit(0);
X}
X//C
X//C-----------------------------------------------------------------------
X//C SUBROUTINE: REMEZ
X//C   THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM
X//C   FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS
X//C   FUNCTION WITH A SUM OF COSINES.  INPUTS TO THE SUBROUTINE
X//C   ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE
X//C   DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE
X//C   GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE
X//C   EXTREMAL FREQUENCIES.  THE PROGRAM MINIMIZES THE CHEBYSHEV
X//C   ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL
X//C   FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES
X//C   THE COEFFICIENTS OF THE BEST APPROXIMATION.
X//C-----------------------------------------------------------------------
X//C
X//      SUBROUTINE REMEZ (EDGE,NBANDS)
Xvoid remez( double* edge, int nbands)
X{
X//C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
X//      INTEGER NITER,OUTPUT
X//      INTEGER IEXT(28)
X//      DOUBLE PRECISION EDGE(20)
X//      DOUBLE PRECISION AD(28),ALPHA(28),X(28),Y(28)
X//      DOUBLE PRECISION DES(488),GRID(488),WT(488)
X//      DOUBLE PRECISION A(66),P(65),Q(65),GTEMP
X//      DOUBLE PRECISION PI2,DEV,DTEMP,DNUM,DDEN,YNZ,COMP,Y1,ERR
X//      DOUBLE PRECISION FSH,CN,DELF,FT,XT,XT1,XE,DAK,DK,AA,BB,DEVL
X//      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
X//      COMMON /OOPS/NITER,OUTPUT
Xdouble a[66+2], p[65+2], q[65+2],
X       gtemp, ynz, comp, y1, err, fsh, cn, aa, bb, ft;
Xdouble xt1, xe, dak, dk, devl, dtemp, dnum, dden;
Xint jb, j, k, l;
Xint nm1, nzmj, nf2j, nf3j, kup, ix, i;
X
Xint    jet, nu, k1, knz, klow, nut, nut1, luck, nzzmj, kn, kkk, nband2, nf1j;
Xint    itrmax, nzz, jchnge, jp1, jm1, nz, jxt;
X//C
X//C  THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25
X//C
X//C      WRITE (OUTPUT,101)
X//C101   FORMAT(' ','REMEZ CALLED')
X//	cout<<" Remez called \n";
X//      ITRMAX=25
X	itrmax=25;
X//      DEVL=-1.0
X	devl=-1.0;
X//      NZ=NFCNS+1
X	nz=nfcns+1;
X//      NZZ=NFCNS+2
X	nzz=nfcns+2;
X//      NITER=0
X	niter=0;
X//100   CONTINUE
Xll100:
X//      IEXT(NZZ)=NGRID+1
X	iext[nzz]=ngrid+1;
X//      NITER=NITER+1
X	niter=niter+1;
X//      IF(NITER.GT.ITRMAX) GO TO 400
X	if(niter>itrmax) goto ll400;
X//      DO 110 J=1,NZ
X	for (j=1; j<=nz; j++){
X//      JXT=IEXT(J)
X		jxt=iext[j];
X//      DTEMP=GRID(JXT)
X		dtemp=grid[jxt];
X//      DTEMP=DCOS(DTEMP*PI2)
X		dtemp=cos(dtemp*pi2);
X//110   X(J)=DTEMP
X		x[j]=dtemp;}
X//      JET=(NFCNS-1)/15+1
X	jet=(nfcns-1)/15+1;
X//      DO 120 J=1,NZ
X	for (j=1; j<=nz; j++){
X//120   AD(J)=D(J,NZ,JET)
X		ad[j]=d(j, nz, jet);}
X//      DNUM=0.0
X	dnum=0.0;
X//      DDEN=0.0
X	dden=0.0;
X//      K=1
X	k=1;
X//      DO 130 J=1,NZ
X	for (j=1; j<=nz; j++){
X//      L=IEXT(J)
X		l=iext[j];
X//      DTEMP=AD(J)*DES(L)
X		dtemp=ad[j]*des[l];
X//      DNUM=DNUM+DTEMP
X		dnum=dnum+dtemp;
X//C
X//C      WRITE (OUTPUT,1222)WT(L),L
X//C1222  FORMAT(' ','WT(L),L ',F15.9,I4)
X//C
X//      DTEMP=DBLE(double(K))*AD(J)/WT(L)
X//		cout<<wt[l]<<"=wt["<<l<<"]\n";
X		dtemp=(k) * ad[j]/wt[l];
X//		cout<<"after divide with wt[l]\n";
X//      DDEN=DDEN+DTEMP
X		dden=dden+dtemp;
X//130   K=-K
X		k=-k;  }
X//C
X//C      WRITE (OUTPUT,1233)DDEN
X//C1233  FORMAT(' ','DDEN ',F15.9)
X//C
X//      DEV=DNUM/DDEN
X//	cout<<dden<<"=dden\n";
X	dev=dnum/dden;
X//C     WRITE(OUTPUT,131) DEV
X//C131   FORMAT(1X,'Deviation = ',F12.9)
X//      NU=1
X	nu=1;
X//      IF(DEV.GT.0.0) NU=-1
X	if (dev>0.0) nu=-1;
X//      DEV=-DBLE(double(NU))*DEV
X	dev=-(nu)*dev;
X//      K=NU
X	k=nu;
X//      DO 140 J=1,NZ
X	for (j=1; j<=nz; j++){
X//      L=IEXT(J)
X		l=iext[j];
X//C
X//C      WRITE (OUTPUT,1223)WT(L),L
X//C1223  FORMAT(' ','2:WT(L),L ',F15.9,I4)
X//C
X//      DTEMP=DBLE(double(K))*DEV/WT(L)
X		dtemp=(k)*dev/wt[l];
X//      Y(J)=DES(L)+DTEMP
X		y[j]=des[l]+dtemp;
X//140   K=-K
X		k=-k;}
X//      IF(DEV.GE.DEVL) GO TO 150
X	if(dev>=devl) goto ll150;
X//      CALL OUCH
X	ouch();
X//      GO TO 400
X	goto ll400;
X//150   DEVL=DEV
Xll150:
X	devl=dev;
X//      JCHNGE=0
X	jchnge=0;
X//      K1=IEXT(1)
X	k1=iext[1];
X//      KNZ=IEXT(NZ)
X	knz=iext[nz];
X//      KLOW=0
X	klow=0;
X//      NUT=-NU
X	nut=-nu;
X//      J=1
X	j=1;
X//C
X//C  SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST
X//C  APPROXIMATION
X//C
X//200   IF(J.EQ.NZZ) YNZ=COMP
Xll200:
X	if(j==nzz) ynz=comp;
X//      IF(J.GE.NZZ) GO TO 300
X	if(j>=nzz) goto ll300;
X//      KUP=IEXT(J+1)
X	kup=iext[j+1];
X//      L=IEXT(J)+1
X	l=iext[j]+1;
X//      NUT=-NUT
X	nut=-nut;
X//      IF(J.EQ.2) Y1=COMP
X	if(j==2)   y1=comp;
X//      COMP=DEV
X	comp=dev;
X//      IF(L.GE.KUP) GO TO 220
X	if(l>=kup)   goto ll220;
X//      ERR=GEE(L,NZ)
X	err=gee(l,nz);
X//      ERR=(ERR-DES(L))*WT(L)
X	err=(err-des[l])*wt[l];
X//      DTEMP=DBLE(double(NUT))*ERR-COMP
X	dtemp=(nut)*err-comp;
X//      IF(DTEMP.LE.0.0) GO TO 220
X	if(dtemp<=0.0) goto ll220;
X//      COMP=DBLE(double(NUT))*ERR
X	comp=(nut)*err;
X//210   L=L+1
Xll210:
X	l=l+1;
X//      IF(L.GE.KUP) GO TO 215
X	if(l>=kup)goto ll215;
X//      ERR=GEE(L,NZ)
X	err=gee(l,nz);
X//      ERR=(ERR-DES(L))*WT(L)
X	err=(err-des[l])*wt[l];
X//      DTEMP=DBLE(double(NUT))*ERR-COMP
X	dtemp=(nut)*err-comp;
X//      IF(DTEMP.LE.0.0) GO TO 215
X	if (dtemp<=0.0) goto ll215;
X//      COMP=DBLE(double(NUT))*ERR
X	comp=(nut)*err;
X//      GO TO 210
X	goto ll210;
X//215   IEXT(J)=L-1
Xll215:
X	iext[j]=l-1;
X//      J=J+1
X	j=j+1;
X//      KLOW=L-1
X	klow=l-1;
X//      JCHNGE=JCHNGE+1
X	jchnge=jchnge+1;
X//      GO TO 200
X	goto ll200;
X//220   L=L-1
Xll220:	l=l-1;
X//225   L=L-1
Xll225:	l=l-1;
X//      IF(L.LE.KLOW) GO TO 250
X	if(l<=klow) goto ll250;
X//      ERR=GEE(L,NZ)
X	err=gee(l,nz);
X//      ERR=(ERR-DES(L))*WT(L)
X	err=(err-des[l])*wt[l];
X//      DTEMP=DBLE(double(NUT))*ERR-COMP
X	dtemp=(nut)*err-comp;
X//      IF(DTEMP.GT.0.0) GO TO 230
X	if(dtemp>0.0) goto ll230;
X//      IF(JCHNGE.LE.0) GO TO 225
X	if(jchnge<=0) goto ll225;
X//      GO TO 260
X	goto ll260;
X//230   COMP=DBLE(double(NUT))*ERR
Xll230:
X	comp=(nut)*err;
X//235   L=L-1
Xll235:
X	l=l-1;
X//      IF(L.LE.KLOW) GO TO 240
X	if(l<=klow) goto ll240;
X//      ERR=GEE(L,NZ)
X	err=gee(l,nz);
X//      ERR=(ERR-DES(L))*WT(L)
X	err=(err-des[l])*wt[l];
X//      DTEMP=DBLE(double(NUT))*ERR-COMP
X	dtemp=nut*err-comp;
X//      IF(DTEMP.LE.0.0) GO TO 240
X	if(dtemp<=0.0) goto ll240;
X//      COMP=DBLE(double(NUT))*ERR
X	comp=(nut)*err;
X//      GO TO 235
X	goto ll235;
X//240   KLOW=IEXT(J)
Xll240:
X	klow=iext[j];
X//      IEXT(J)=L+1
X	iext[j]=l+1;
X//      J=J+1
X	j=j+1;
X//      JCHNGE=JCHNGE+1
X	jchnge=jchnge+1;
X//      GO TO 200
X	goto ll200;
X//250   L=IEXT(J)+1
Xll250:
X	l=iext[j]+1;
X//      IF(JCHNGE.GT.0) GO TO 215
X	if(jchnge>0) goto ll215;
X//255   L=L+1
Xll255:
X	l=l+1;
X//      IF(L.GE.KUP) GO TO 260
X	if(l>=kup) goto ll260;
X//      ERR=GEE(L,NZ)
X	err=gee(l,nz);
X//      ERR=(ERR-DES(L))*WT(L)
X	err=(err-des[l])*wt[l];
X//      DTEMP=DBLE(double(NUT))*ERR-COMP
X	dtemp=(nut)*err-comp;
X//      IF(DTEMP.LE.0.0) GO TO 255
X	if(dtemp<=0.0) goto ll255;
X//      COMP=DBLE(double(NUT))*ERR
X	comp=(nut)*err;
X//      GO TO 210
X	goto ll210;
X//260   KLOW=IEXT(J)
Xll260:
X	klow=iext[j];
X//      J=J+1
X	j=j+1;
X//      GO TO 200
X	goto ll200;
X//300   IF(J.GT.NZZ) GO TO 320
Xll300:
X	if(j>nzz) goto ll320;
X//      IF(K1.GT.IEXT(1)) K1=IEXT(1)
X	if(k1>iext[1]) k1=iext[1];
X//      IF(KNZ.LT.IEXT(NZ)) KNZ=IEXT(NZ)
X	if(knz<iext[nz]) knz=iext[nz];
X//      NUT1=NUT
X	nut1=nut;
X//      NUT=-NU
X	nut=-nu;
X//      L=0
X	l=0;
X//      KUP=K1
X	kup=k1;
X//      COMP=YNZ*(1.00001D0)
X	comp=ynz*(1.00001);
X//      LUCK=1
X	luck=1;
X//310   L=L+1
Xll310:
X	l=l+1;
X//      IF(L.GE.KUP) GO TO 315
X	if(l>=kup) goto ll315;
X//      ERR=GEE(L,NZ)
X	err=gee(l,nz);
X//      ERR=(ERR-DES(L))*WT(L)
X	err=(err-des[l])*wt[l];
X//      DTEMP=DBLE(double(NUT))*ERR-COMP
X	dtemp=(nut)*err-comp;
X//      IF(DTEMP.LE.0.0) GO TO 310
X	if(dtemp<=0.0) goto ll310;
X//      COMP=DBLE(double(NUT))*ERR
X	comp=(nut)*err;
X//      J=NZZ
X	j=nzz;
X//      GO TO 210
X	goto ll210;
X//315   LUCK=6
Xll315:
X	luck=6;
X//      GO TO 325
X	goto ll325;
X//320   IF(LUCK.GT.9) GO TO 350
Xll320:
X	if(luck>9) goto ll350;
X//      IF(COMP.GT.Y1) Y1=COMP
X	if(comp>y1) y1=comp;
X//      K1=IEXT(NZZ)
X	k1=iext[nzz];
X//325   L=NGRID+1
Xll325:
X	l=ngrid+1;
X//      KLOW=KNZ
X	klow=knz;
X//      NUT=-NUT1
X	nut=-nut1;
X//      COMP=Y1*(1.00001D0)
X	comp=y1*(1.00001);
X//330   L=L-1
Xll330:
X	l=l-1;
X//      IF(L.LE.KLOW) GO TO 340
X	if(l<=klow) goto ll340;
X//      ERR=GEE(L,NZ)
X	err=gee(l,nz);
X//      ERR=(ERR-DES(L))*WT(L)
X	err=(err-des[l])*wt[l];
X//      DTEMP=DBLE[double(NUT))*ERR-COMP
X	dtemp=(nut)*err-comp;
X//      IF(DTEMP.LE.0.0) GO TO 330
X	if(dtemp<=0.0) goto ll330;
X//      J=NZZ
X	j=nzz;
X//      COMP=DBLE(double(NUT))*ERR
X	comp=(nut)*err;
X//      LUCK=LUCK+10
X	luck=luck+10;
X//      GO TO 235
X	goto ll235;
X//340   IF(LUCK.EQ.6) GO TO 370
Xll340:
X	if(luck==6) goto ll370;
X//      DO 345 J=1,NFCNS
X	for (j=1;j<=nfcns; j++){
X//      NZZMJ=NZZ-J
X		nzzmj = nzz-j;
X//      NZMJ=NZ-J
X		nzmj=nz-j;
X//345   IEXT(NZZMJ)=IEXT(NZMJ)
X		iext[nzzmj]=iext[nzmj]; }
X//      IEXT(1)=K1
X	iext[1]=k1;
X//      GO TO 100
X	goto ll100;
X//350   KN=IEXT(NZZ)
Xll350:
X	kn=iext[nzz];
X//      DO 360 J=1,NFCNS
X	for (j=1; j<=nfcns; j++){
X//360   IEXT[J)=IEXT(J+1)
X		iext[j]=iext[j+1];}
X//      IEXT(NZ)=KN
X	iext[nz]=kn;
X//      GO TO 100
X	goto ll100;
X//370   IF(JCHNGE.GT.0) GO TO 100
Xll370:
X	if(jchnge>0) goto ll100;
X//C
X//C  CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
X//C  USING THE INVERSE DISCRETE FOURIER TRANSFORM
X//C
X//400   CONTINUE
Xll400:
X//      NM1=NFCNS-1
X	nm1=nfcns-1;
X//      FSH=1.0D-06
X	fsh=1.0e-6;
X//      GTEMP=GRID(1)
X	gtemp=grid[1];
X//      X(NZZ)=-2.0D0
X	x[nzz]=-2.0;
X//      CN=DBLE(double(2*NFCNS-1))
X	cn=(2*nfcns-1);
X//C
X//C      WRITE (OUTPUT,1224)CN
X//C1224  FORMAT(' ','CN ',F15.9)
X//C
X//      DELF=1.0D0/CN
X	delf=1.0/cn;
X//      L=1
X	l=1;
X//      KKK=0
X	kkk=0;
X///C     IF((GRID(1).LT.0.01D0).AND.(GRID(NGRID).GT.0.5D0)) KKK=1
X	if ((grid[1]<0.01)&&(grid[ngrid]>0.5)) kkk=1;
X//C	The following is in text by Rabiner and Gold
X//       NBAND2=NBANDS*2
X//	nband2=nbands*2;
X//       IF((EDGE(1).EQ.0.0D0).AND.(EDGE(NBAND2).EQ.0.5D0)) KKK=1
X//	if((edge[1]==0.0)&&(edge[nband2]==0.5))kkk=1;
X//      IF(NFCNS.LE.3) KKK=1
X	if(nfcns<=3) kkk=1;
X//      IF(KKK.EQ.1) GO TO 405
X	if(kkk==1) goto ll405;
X//      DTEMP=DCOS(PI2*GRID(1))
X	dtemp=cos(pi2*grid[1]);
X//      DNUM=DCOS(PI2*GRID(NGRID))
X	dnum=cos(pi2*grid[ngrid]);
X//C
X//C      WRITE (O,1225)DTEMP,DNUM
X//C1225  FORMAT(' ','DTEMP,DNUM ',2F15.9)
X//C
X//      AA=2.0D0/(DTEMP-DNUM)
X//	cout<<(dtemp-dnum)<<"=dtemp-dnum\n";
X	aa=2.0/(dtemp-dnum);
X//      BB=-(DTEMP+DNUM)/(DTEMP-DNUM)
X	bb=-(dtemp+dnum)/(dtemp-dnum);
X//405   CONTINUE
Xll405:
X//      DO 430 J=1,NFCNS
X	for (j=1; j<=nfcns; j++){
X//      FT=DBLE(double(J-1))*DELF
X		ft=(j-1)*delf;
X//      XT=DCOS(PI2*FT)
X		xt=cos(pi2*ft);
X//      IF(KKK.EQ.1) GO TO 410
X		if(kkk==1) goto ll410;
X//C
X//C      WRITE (OUTPUT,1226)AA
X//C1226  FORMAT(' ','AA ',F15.9)
X//C
X//      XT=(XT-BB)/AA
X//		cout<<aa<<"=aa\n";
X		xt=(xt-bb)/aa;
X//C      WRITE (OUTPUT,1226)XT
X//C1226  FORMAT(' XT:',F15.9)
X//	cout<<xt<<"=xt\n";
X//      XT1=DSQRT(DABS(1.0D0-XT*XT))
X//	xt1=sqrt(abs(1.0-xt*xt));
X//      FT=DATAN2(XT1,XT)
X//	ft=atan2(xt1,xt);
X//      IF (XT.LT.0.0D0) FT = PI-FT
X//	if(xt<0.0) ft=pi-ft;
X//      FT = FT/PI2
X//	ft=ft/pi2;
X		ft=acos(xt)/pi2;
X//410   XE=X(L)
Xll410:
X		xe=x[l];
X//      IF(XT.GT.XE) GO TO 420
X		if(xt>xe) goto ll420;
X//      IF((XE-XT).LT.FSH) GO TO 415
X		if((xe-xt)<fsh) goto ll415;
X//      L=L+1
X		l=l+1;
X//      GO TO 410
X		goto ll410;
X//415   A(J)=Y(L)
Xll415:
X		a[j]=y[l];
X//      GO TO 425
X		goto ll425;
X//420   IF((XT-XE).LT.FSH) GO TO 415
Xll420:
X		if((xt-xe)<fsh) goto ll415;
X//      GRID(1)=FT
X		grid[1]=ft;
X//      A(J)=GEE(1,NZ)
X		a[j]=gee(1,nz);
X//425   CONTINUE
Xll425:
X//      IF(L.GT.1) L=L-1
X		if(l>1) l=l-1;
X//430   CONTINUE
X	}
X//      GRID(1)=GTEMP
X	grid[1]=gtemp;
X//C
X//C      WRITE (OUTPUT,1227)CN
X//C1227  FORMAT(' ','2:CN ',F15.9,I4)
X//C
X//      DDEN=PI2/CN
X	dden=pi2/cn;
X//      DO 510 J=1,NFCNS
X	for (j=1; j<=nfcns; j++){
X//      DTEMP=0.0
X		dtemp=0.0;
X//      DNUM=J-1
X		dnum=j-1;
X//      DNUM=DNUM*DDEN
X		dnum=dnum*dden;
X//      IF(NM1.LT.1) GO TO 505
X		if(nm1<1) goto ll505;
X//      DO 500 K=1,NM1
X	for (k=1; k<=nm1; k++){
X//      DAK=A(K+1)
X		dak=a[k+1];
X//      DK=K
X		dk=k;
X//500   DTEMP=DTEMP+DAK*DCOS(DNUM*DK)
X		dtemp=dtemp+dak*cos(dnum*dk); }
X//505   DTEMP=2.0D0*DTEMP+A(1)
Xll505:
X	dtemp=2.0*dtemp+a[1];
X//510   ALPHA(J)=DTEMP
Xll510:
X	alpha[j]=dtemp;             }
X//      DO 550 J=2,NFCNS
X	for (j=2; j<=nfcns; j++){
X//550   ALPHA(J)=2.0D0*ALPHA(J)/CN
X		alpha[j]=2.0*alpha[j]/cn;}
X//      ALPHA(1)=ALPHA(1)/CN
X	alpha[1]=alpha[1]/cn;
X//      IF(KKK.EQ.1) GO TO 545
X	if(kkk==1) goto ll545;
X//      P(1)=2.0D0*ALPHA(NFCNS)*BB+ALPHA(NM1)
X	p[1]=2.0*alpha[nfcns]*bb+alpha[nm1];
X//      P(2)=2.0D0*AA*ALPHA(NFCNS)
X	p[2]=2.0*aa*alpha[nfcns];
X//      Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)
X	q[1]=alpha[nfcns-2]-alpha[nfcns];
X//      DO 540 J=2,NM1
X	for (j=2;j<=nm1;j++){
X//      IF(J.LT.NM1) GO TO 515
X		if(j<nm1)goto ll515;
X//      AA=0.5D0*AA
X		aa=0.5*aa;
X//      BB=0.5D0*BB
X		bb=0.5*bb;
X//515   CONTINUE
Xll515:
X//      P(J+1)=0.0D0
X		p[j+1]=0.0;
X//      DO 520 K=1,J
X		for (k=1;k<=j;k++){
X//      A(K)=P(K)
X			a[k]=p[k];
X//520   P(K)=2.0D0*BB*A(K)
X			p[k]=2.0*bb*a[k];}
X//      P(2)=P(2)+A(1)*2.0D0*AA
X		p[2]=p[2]+a[1]*2.0*aa;
X//      JM1=J-1
X		jm1=j-1;
X//      DO 525 K=1,JM1
X		for (k=1;k<=jm1;k++){
X//525   P(K)=P(K)+Q(K)+AA*A(K+1)
X			p[k]=p[k]+q[k]+aa*a[k+1]; }
X//      JP1=J+1
X		jp1=j+1;
X//      DO 530 K=3,JP1
X		for (k=3; k<= jp1;k++){
X//530   P(K)=P(K)+AA*A(K-1)
X			p[k]=p[k]+aa*a[k-1]; }
X//      IF(J.EQ.NM1) GO TO 540
X		if(j==nm1) goto ll540;
X//      DO 535 K=1,J
X		for (k=1;k<=j;k++){
X//535   Q(K)=-A(K)
X			q[k]=-a[k]; }
X//      NF1J=NFCNS-1-J
X		nf1j=nfcns-1-j;
X//      Q(1)=Q(1)+ALPHA(NF1J)
X		q[1]=q[1]+alpha[nf1j];
X//540   CONTINUE
Xll540:	;}
X//      DO 543 J=1,NFCNS
X	for (j=1;j<=nfcns;j++){
X//543   ALPHA(J)=P(J)
X		alpha[j]=p[j]; }
X//545   CONTINUE
Xll545:
X//      IF(NFCNS.GT.3) GOTO 9000
X	if(nfcns>3) return;
X//      ALPHA(NFCNS+1)=0.0
X	alpha[nfcns+1]=0.0;
X//      ALPHA(NFCNS+2)=0.0
X	alpha[nfcns+2]=0.0;
X//9000  CONTINUE
X//      RETURN
X//      END
X	return;
X}
X//C
X//C-----------------------------------------------------------------------
X//C FUNCTION: D
X//C   FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION
X//C   COEFFICIENTS FOR USE IN THE FUNCTION GEE.
X//C-----------------------------------------------------------------------
X//C
X//      DOUBLE PRECISION FUNCTION D(K,N,M)
Xdouble d(int k, int n, int m)
X{
X//C     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
X//      INTEGER IEXT(28)
X//C
X//      DOUBLE PRECISION AD(28),ALPHA(28),X(28),Y(28),Q
X//      DOUBLE PRECISION DES(488),GRID(488),WT(488)
X//      DOUBLE PRECISION PI2,DEV
X//      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
X//      D=1.0D0
Xdouble ret_d, q;
Xint l, j;
X//
X     ret_d = 1.0;
X//      Q=X(K)
X     q=x[k];
X//      DO 3 L=1,M
X     for (l=1;l<=m;l++) {
X//      DO 2 J=L,N,M
X	     for (j=l;j<=n;j+=m){
X//      IF(J-K)1,2,1
X	     if ((j-k)!=0) goto l1;
X	     goto l2;
X//1     D=2.0D0*D*(Q-X(J))
Xl1:
X	     ret_d=2.0*ret_d*(q-x[j]);
X//2     CONTINUE
Xl2:	     ;}
X//3     CONTINUE
X	}
X//C     IF (DABS(D).LT.0.00001) D = 0.00001
X	if (fabs(ret_d)<0.00001) ret_d=0.00001;
X//	if (ret_d < 0.0) {
X//	  if (abs(ret_d)<0.00001) ret_d=-0.00001;}
X//	else
X//	  {if (abs(ret_d)<0.00001) ret_d=0.00001;}
X//      D=1.0D0/D
X	ret_d=1.0/ret_d;
X//      RETURN
X//	cout<<"d="<<ret_d<<"\n";
X	return(ret_d);
X//      END
X}
X//C
X//C-----------------------------------------------------------------------
X//C FUNCTION: GEE
X//C   FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE
X//C   LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM
X//C-----------------------------------------------------------------------
X//C
X//      DOUBLE PRECISION FUNCTION GEE(K,N)
Xdouble gee(int k, int n)
X{
X//C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
X//      INTEGER IEXT(28)
X//      DOUBLE PRECISION AD(28),ALPHA(28),X(28),Y(28)
X//      DOUBLE PRECISION DES(488),GRID(488),WT(488),P,XF,PI2,C
X//      DOUBLE PRECISION DEV,D
X//      COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID
X//      P=0.0D0
Xdouble p, xf, d, c;
X	p=0.0;
X//      XF=GRID(K)
X	xf=grid[k];
X//      XF=DCOS(PI2*XF)
X	xf=cos(pi2*xf);
X//      D=0.0D0
X	d = 0.0;
X//      DO 1 J=1,N
X	for (j=1;j<=n;j++){
X//      C=XF-X(J)
X		c=xf-x[j];
X//      C=AD(J)/C
X		c=ad[j]/c;
X//      D=D+C
X		d=d+c;
X//1     P=P+C*Y(J)
X		p=p+c*y[j]; }
X//      GEE=P/D
X	return(p/d);
X//      RETURN
X//      END
X}
X//C
X//C-----------------------------------------------------------------------
X//C SUBROUTINE: OUCH
X//C   WRITES AN ERROR MESSAGE WHEN THE ALGORITHM FAILS TO
X//C   CONVERGE.  THERE SEEM TO BE TWO CONDITIONS UNDER WHICH
X//C   THE ALGORITHM FAILS TO CONVERGE: (1) THE INITIAL
X//C   GUESS FOR THE EXTREMAL FREQUENCIES IS SO POOR THAT
X//C   THE EXCHANGE ITERATION CANNOT GET STARTED, OR
X//C   (2) NEAR THE TERMINATION OF A CORRECT DESIGN,
X//C   THE DEVIATION DECREASES DUE TO ROUNDING ERRORS
X//C   AND THE PROGRAM STOPS.  IN THIS LATTER CASE THE
X//C   FILTER DESIGN IS PROBABLY ACCEPTABLE, BUT SHOULD
X//C   BE CHECKED BY COMPUTING A FREQUENCY RESPONSE.
X//C-----------------------------------------------------------------------
X//C
Xvoid ouch(void){
X//      SUBROUTINE OUCH
X//      INTEGER NITER,OUTPUT
X//      COMMON /OOPS/NITER,OUTPUT
X//      WRITE(OUTPUT,1)NITER
X//1     FORMAT(' ************ FAILURE TO CONVERGE **********'/
X//     -       ' Probable cause is machine rounding error.'/
X//     -       ' Number of interations equals ',I4/
X//     -       ' If the number of iterations exceeds 3, the design may'
X//     -       ' be correct,'/
X//     -       ' but should be verified with an fft!')
X//      RETURN
X//      END
X cout<<"Failed to converge. Rounding error is suspect after"<<niter<<"iterations.\n";
X }
X
END_OF_FILE
  if test 44237 -ne `wc -c <'fir4.cpp'`; then
    echo shar: \"'fir4.cpp'\" unpacked with wrong size!
  fi
  # end of 'fir4.cpp'
fi
echo shar: End of archive.
exit 0
exit 0 # Just in case...
