C REMEZ C C----------------------------------------------------------------------- C SUBROUTINE: REMEZ C THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM C FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS C FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE C ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE C DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE C GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE C EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV C ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL C FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES C THE COEFFICIENTS OF THE BEST APPROXIMATION. C C REFERENCE: THEORY AND APPLICATION OF DIGITAL SIGNAL PROCESSING, C BY RABINER AND GOLD, PRENTICE-HALL, 1975, PP. 136-138. C----------------------------------------------------------------------- C C ARGUMENTS: C C IEXT ARRAY OF EXTREMAL FREQUENCIES. C AD ARRAY OF ??? C ALPHA ARRAY OF COEFFICIENTS OF THE OPTIMAL APPROX. C C IF THE ERROR DECREASES DURING THE ITERATIONS, NITER IS C RETURNED NEGATED TO INDICATE THIS. THIS MAY HAPPEN IF C (1) THE INITIAL VALUES GIVEN IN IEXT ARE TOO FAR OFF C (2) MACHINE ROUNDOFF ERROR C CASE (2) FAILURE IS GENERALLY AT A GOOD APPROXIMATION C WHILE CASE (1) FAILURES MEAN THE ALGORITHM COULD NOT C GET STARTED. C THIS CHECK FOR NON-CONVERGENCE IS OPERATIVE ONLY IF C ITRMAX EXCEEDS ONE. C C LET THE NUMBER OF COSINES IN THE APPROXIMATION BE NALPHA. C THEN THE REQUIRED DIMENSIONS OF THE ARRAYS IS THE FOLLOWING: C ARRAYS IEXT, AD, ALPHA, X, Y, H AT LEAST NALPHA + 2. C ARRAYS DES, GRID, AND WT MUST BE LGRID*(NALPHA + 2), C WHERE LGRID IS THE NUMBER OF GRID POINTS PER COSINE. C C AN EXAMPLE SET OF DIMENSIONS FOR APPROXIMATION BY 66 COSINES C WITH 16 GRID POINTS PER COSINE IS AS FOLLOWS: C C DIMENSION IEXT(66),AD(66),ALPHA(66),X(66),Y(66) C DIMENSION DES(1045),GRID(1045),WT(1045) C DIMENSION A(66),P(65),Q(65) C SUBROUTINE REMEZ(NFCNS,ITRMAX,NGRID,DES,WT,IEXT,ALPHA,DEV, & NITER,GRID,AD,X,Y,A,P,Q) DIMENSION IEXT(1),AD(1),ALPHA(1),X(1),Y(1) DIMENSION DES(1),GRID(1),WT(1) DIMENSION A(1),P(1),Q(1) DEVL=-1.0 NZ=NFCNS+1 NZZ=NFCNS+2 NITER=0 C PI2=8.0D00*DATAN(1.0D00) PI2=8.0*ATAN(1.0) C C THE FOLLOWING PARAMETER ALLOWS CONTINUATION THROUGH NON-CONVERGENCE C INDICATED BY THE MAXIMUM DEVIATION NOT STRICTLY INCREASING EACH C ITERATION. C TOLER = 0.000001 C 100 CONTINUE IEXT(NZZ)=NGRID+1 NITER=NITER+1 IF(NITER.GT.ITRMAX) GO TO 400 C C FILL ARRAY X WITH THE COSINE OF THE EXTREMAL FREQUENCIES C (NEEDED FOR INTERPOLATION AND DEVIATION CALCULATION) C DO 110 J=1,NZ JXT=IEXT(J) DTEMP=GRID(JXT) DTEMP=COS(DTEMP*PI2) 110 X(J)=DTEMP C C PLACE THE LAGRANGE INTERPOLATION COEFFICIENTS FOR THESE C FREQUENCIES INTO ARRAY AD C JET=(NFCNS-1)/15+1 DO 120 J=1,NZ 120 AD(J)=D(J,NZ,JET,X) C C CALCULATE THE DEVIATION FOR THE CURRENT EXTREMAL FREQUENCIES C VIA EQUATION (3.121) OF RABINER AND GOLD. (OBTAINED BY APPLYING C CRAMER'S RULE TO THE LINEAR SYSTEM OF EQUATIONS (3.120)). C DNUM=0.0 DDEN=0.0 K=1 DO 130 J=1,NZ L=IEXT(J) DTEMP=AD(J)*DES(L) DNUM=DNUM+DTEMP DTEMP=FLOAT(K)*AD(J)/WT(L) DDEN=DDEN+DTEMP 130 K=-K DEV=DNUM/DDEN C C FILL ARRAY Y WITH THE VALUES THROUGH WHICH THE CURRENT C APPROXIMATION WILL INTERPOLATE. C 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 C C CHECK FOR FAILURE TO CONVERGE C IF(DEV.GT.DEVL-TOLER) GO TO 150 NITER = -NITER GO TO 400 150 DEVL=DEV C C SEARCH FOR THE EXTREMAL FREQUENCIES OF THE NEXT C APPROXIMATION WHICH ARE THE LOCATIONS OF MAXIMUM C ERROR GIVEN THE CURRENT APPROXIMATION. C 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,X,Y,GRID,AD) 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,X,Y,GRID,AD) 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,X,Y,GRID,AD) 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,X,Y,GRID,AD) 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,X,Y,GRID,AD) 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,X,Y,GRID,AD) 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,X,Y,GRID,AD) ERR=(ERR-DES(L))*WT(L) DTEMP=FLOAT(NUT)*ERR-COMP IF(DTEMP.LE.0.0) GO TO 330 J=NZZ COMP=FLOAT(NUT)*ERR LUCK=LUCK+10 GO TO 235 340 IF(LUCK.EQ.6) GO TO 370 DO 345 J=1,NFCNS NZZMJ=NZZ-J NZMJ=NZ-J 345 IEXT(NZZMJ)=IEXT(NZMJ) IEXT(1)=K1 GO TO 100 350 KN=IEXT(NZZ) DO 360 J=1,NFCNS 360 IEXT(J)=IEXT(J+1) IEXT(NZ)=KN GO TO 100 370 IF(JCHNGE.GT.0) GO TO 100 C C CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION C USING THE INVERSE DISCRETE FOURIER TRANSFORM C 400 CONTINUE NM1=NFCNS-1 FSH=1.0E-06 GTEMP=GRID(1) X(NZZ)=-2.0 CN=2*NFCNS-1 DELF=1.0/CN L=1 KKK=0 IF(GRID(1).LT.0.01.AND.GRID(NGRID).GT.0.49) KKK=1 IF(NFCNS.LE.3) KKK=1 IF(KKK.EQ.1) GO TO 405 DTEMP=COS(PI2*GRID(1)) DNUM=COS(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=COS(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,X,Y,GRID,AD) 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*COS(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 C C THE FOLLOWING HAS BEEN DISABLED C C IF(NFCNS.GT.3) RETURN C ALPHA(NFCNS+1)=0.0 C ALPHA(NFCNS+2)=0.0 C RETURN END C D C----------------------------------------------------------------------- C FUNCTION: D C FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION C COEFFICIENTS FOR USE IN THE FUNCTION GEE. C C RETURNED IS C C N*M 1 C PRODUCT ----------- C I=1 X(K) - X(I) C IK C C C IF THE X ARRAY CONTAINS THE COSINE OF THE INTERPOLATION FREQUENCIES, C X(I) = COS(PI2*F(I)), I=1,...,NALPHA, THEN D RETURNS THE C INTERPOLATION COEFFICIENT FOR A POLE AT X(K). C C THE FACTORS N AND M CONTROL THE ORDERING OF THE PRODUCT TERMS. C----------------------------------------------------------------------- C C FUNCTION D(K,N,M,X) DIMENSION X(1) 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 C GEE C----------------------------------------------------------------------- C FUNCTION: GEE C FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE C LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM C----------------------------------------------------------------------- C FUNCTION GEE(K,N,X,Y,GRID,AD) DIMENSION X(1),Y(1),GRID(1),AD(1) P=0.0 XF=GRID(K) XF=COS(6.283185307*XF) D=0.0 DO 1 J=1,N C=XF-X(J) IF(C.EQ.0.)CALL ERRF4(-10) C=AD(J)/C D=D+C 1 P=P+C*Y(J) GEE=P/D RETURN END C