C C----------------------------------------------------------------------- C MAIN PROGRAM: FIR LINEAR PHASE FILTER DESIGN PROGRAM C C AUTHORS: JAMES H. MCCLELLAN 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 MODIFIED BY: S. A. TRETTER C ELECTRICAL ENGINEERING DEPARTMENT C UNIVERSITY OF MARYLAND C COLLEGE PARK, MD. 20742 C INPUT: C NFILT-- FILTER LENGTH C JTYPE-- TYPE OF FILTER C 1 = MULTIPLE PASSBAND/STOPBAND FILTER C 2 = DIFFERENTIATOR C 3 = HILBERT TRANSFORM FILTER C NBANDS-- NUMBER OF BANDS C LGRID-- GRID DENSITY, WILL BE SET TO 16 UNLESS C SPECIFIED OTHERWISE BY A POSITIVE CONSTANT. C C EDGE(2*NBANDS)-- BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND C WITH A MAXIMUM OF 10 BANDS. C C FX(NBANDS)-- DESIRED FUNCTION ARRAY (OR DESIRED SLOPE IF A C DIFFERENTIATOR) FOR EACH BAND. C C WTX(NBANDS)-- WEIGHT FUNCTION ARRAY IN EACH BAND. FOR A C DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY C PROPORTIONAL TO F. C C SAMPLE INPUT DATA SETUP: C 32,1,3,0 C 0.0,0.1,0.2,0.35,0.425,0.5 C 0.0,1.0,0.0 C 10.0,1.0,10.0 C THIS 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 THIS IS THE FILTER IN FIGURE 10. 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(202),AD(202),ALPHA(202),X(202),Y(202) DIMENSION H(202) DIMENSION DES(3300),GRID(3300),WT(3300) DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10) DIMENSION OMEGA(4001),RESPA(4001) DOUBLE PRECISION PI2,PI DOUBLE PRECISION AD,DEV,X,Y,GRID DOUBLE PRECISION GEE,D,EDGE,SFREQ,BWIDTH,DELF DOUBLE PRECISION DES CHARACTER BD1,BD2,BD3,BD4 CHARACTER*60 LSTFIL,COEFL CHARACTER UDAR,UDWF,KPLT,LINDB LOGICAL FLAG,IXST DATA BD1,BD2,BD3,BD4 /'B','A','N','D'/ INPUT=5 IOUT=6 PI=4.0*DATAN(1.0D0) PI2=2.0D00*PI FLAG=.TRUE. C C THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 400, BUT C THIS UPPER LIMIT CAN BE CHANGED BY REDIMENSIONING THE C ARRAYS IEXT, AD, ALPHA, X, Y, H TO BE NFMAX/2 + 2. C THE ARRAYS DES, GRID, AND WT MUST DIMENSIONED C 16(NFMAX/2 + 2). C NFMAX=400 100 CONTINUE JTYPE=0 C C PROGRAM INPUT SECTION C WRITE(*,800) 800 FORMAT(' ENTER LISTING FILENAME: ',\) READ(*,809)LSTFIL 809 FORMAT(A) OPEN(IOUT,FILE=LSTFIL,STATUS='NEW') WRITE(*,808) 808 FORMAT(' ENTER COEFFICIENT STORAGE FILENAME: ',\) READ(*,810)COEFL 810 FORMAT(A) 8800 WRITE(*,8808) 8808 FORMAT(' LINEAR OR DB AMPLITUDE SCALE FOR PLOTS? (L OR D): ',\) READ(*,810) LINDB IF(LINDB.NE.'L'.AND.LINDB.NE.'l'.AND.LINDB.NE.'D'.AND. 1 LINDB.NE.'d') GOTO 8800 WRITE(*,801) 801 FORMAT(' ENTER SAMPLING FREQUENCY (HZ): ',\) READ(*,*)SFREQ WRITE(*,'(A)')' ENTER START AND STOP FREQUENCIES IN HZ FOR' WRITE(*,'(A,\)')' RESPONSE CALCULATION (FSTART,FSTOP): ' READ(*,*) FSTART,FSTOP FSTART=FSTART/SFREQ FSTOP = FSTOP/SFREQ WRITE(*,802) 802 FORMAT(/' FILTER TYPES AVAILABLE:'/ 1 ' 1 MULTIPLE PASSBAND/STOPBAND FILTER'/ 2 ' 2 DIFFERENTIATOR'/ 3 ' 3 HILBERT TRANSFORM'/) 8033 WRITE(*,803) 803 FORMAT(' ENTER: FILTER LENGTH, TYPE, NO. OF BANDS, GRID DENSITY:' 1 ,\) READ(*,*) NFILT,JTYPE,NBANDS,LGRID IF(NFILT.EQ.0) THEN WRITE(*,*) ' ERROR -- FILTER LENGTH MUST BE GREATER THAN 0' GO TO 8033 WRITE(*,*) ' ' ENDIF 110 FORMAT(4I5) IF(NFILT.LE.NFMAX.AND.NFILT.GE.3) GO TO 115 IF(NFILT.GT.NFMAX) THEN WRITE(*,*) CHAR(7) WRITE(*,*) 'FILTER LENGTH TOO LARGE. MUST BE LESS THAN 400' GO TO 8033 WRITE(*,*) ' ' ENDIF 115 IF(NBANDS.LE.0) THEN WRITE(*,*) ' ERROR -- NUMBER OF BANDS MUST BE GREATER THAN 0' WRITE(*,*) ' ' GO TO 8033 ENDIF IF(JTYPE.LT.1.OR.JTYPE.GT.3) THEN WRITE(*,*) ' ERROR -- TYPE MUST BE 1, 2, OR 3' WRITE(*,*) ' ' GO TO 8033 ENDIF C C GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED C OTHERWISE C IF(LGRID.LE.0) LGRID=16 JB=2*NBANDS 9012 WRITE(*,804) 804 FORMAT( ' ENTER THE BAND EDGES (FREQUENCIES IN HERTZ)') READ(*,*)(EDGE(J),J=1,JB) DO 9010 J=1,JB EDGE(J)=EDGE(J)/SFREQ IF(EDGE(J).LT.0.OR.EDGE(J).GT.0.5) WRITE(*,805) J,EDGE(J) 805 FORMAT(' EDGE(',I4,')=',F10.5,' OUT OF BOUNDS') 9010 IF(EDGE(J).LT.0.OR.EDGE(J).GT.0.5) GOTO 9012 C---------------------------------------------------------------------- C FIND TOTAL WIDTH OF SPECIFIED BANDS BWIDTH=0. DO 9005 KEDG=1,JB-1,2 9005 BWIDTH=EDGE(KEDG+1)-EDGE(KEDG)+BWIDTH C---------------------------------------------------------------------- WRITE(*,9001) 9001 FORMAT(' SPECIAL USER DEFINED AMPLITUDE RESPONSE(Y/N)? ',\) READ(*,810)UDAR WRITE(*,9003) 9003 FORMAT(' SPECIAL USER DEFINED WEIGHTING FUNCTION(Y/N)? ',\) READ(*,810)UDWF IF(UDAR.EQ.'Y'.OR.UDAR.EQ.'y') GO TO 9002 WRITE(*,806) 806 FORMAT(/' ENTER (SEPARATED BY COMMAS):'/ 1' 1. VALUE FOR EACH BAND FOR MULTIPLE PASS/STOP BAND FILTERS'/ 2' 2. SLOPES FOR DIFFERENTIATOR (GAIN = Ki*f -> SLOPE = Ki'/ 3' WHERE Ki = SLOPE OF i th BAND, f IN HERTZ)'/ 4' 3. MAGNITUDE OF DESIRED VALUE FOR HILBERT TRANSFORM') READ(*,*)(FX(J),J=1,NBANDS) IF(JTYPE.EQ.3) THEN DO 8100 J=1,NBANDS 8100 FX(J) = - FX(J) ENDIF 9002 IF(UDWF.EQ.'Y'.OR.UDWF.EQ.'y') GOTO 9004 WRITE(*,807) 807 FORMAT(' ENTER WEIGHT FOR EACH BAND. (FOR A DIFFERENTIATOR'/ 2 ' THE WEIGHT FUNCTION GENERATED BY THE PROGRAM FOR THE i th'/ 3 ' BAND IS WT(i)/f WHERE WT(i) IS THE ENTERED BAND WEIGHT AND'/ 4 ' f IS IN HERTZ.)') READ(*,*)(WTX(J),J=1,NBANDS) 9004 IF(JTYPE.GT.0.AND.JTYPE.LE.3) GO TO 125 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 C SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID C IS NUMBER OF FUNCTIONS*GRID DENSITY ~ FILTER LENGTH*DENSITY/2 C GRID(1)=EDGE(1) DELF=LGRID*NFCNS IF(DELF.GT.3300.) THEN WRITE(*,*) ' FILTER LENGTH*GRID DENSITY too large.' WRITE(*,*) ' Array dimensions exceeded.' STOP ENDIF DELF=BWIDTH/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 C CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT C FUNCTION ON THE GRID C IF(UDAR.NE.'Y'.AND.UDAR.NE.'y')DES(J)=EFF(TEMP,FX,WTX, 1 LBAND,JTYPE,SFREQ) IF(UDAR.EQ.'Y'.OR.UDAR.EQ.'y') DES(J)=EFF1(TEMP,FX,WTX, 1 LBAND,SFREQ,F1,F2,FLAG) IF(UDWF.NE.'Y'.AND.UDWF.NE.'y') WT(J)=WATE(TEMP,FX,WTX, 1 LBAND,JTYPE,SFREQ) IF(UDWF.EQ.'Y'.OR.UDWF.EQ.'y') WT(J)=WATE1(TEMP,FX,WTX, 1 LBAND,SFREQ,F1,F2,FLAG) J=J+1 GRID(J)=TEMP+DELF IF(GRID(J).GT.FUP) GO TO 150 GO TO 145 150 GRID(J-1)=FUP IF(UDAR.EQ.'Y'.OR.UDAR.EQ.'y') DES(J-1)=EFF1(FUP,FX,WTX, 1 LBAND,SFREQ,F1,F2,FLAG) IF(UDAR.NE.'Y'.AND.UDAR.NE.'y') DES(J-1)=EFF(FUP,FX,WTX, 1 LBAND,JTYPE,SFREQ) IF(UDWF.EQ.'Y'.OR.UDWF.EQ.'y') WT(J-1)=WATE1(FUP,FX,WTX, 1 LBAND,SFREQ,F1,F2,FLAG) IF(UDWF.NE.'Y'.AND.UDWF.NE.'y') WT(J-1)=WATE(FUP,FX,WTX, 1 LBAND,JTYPE,SFREQ) 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 C SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT C TO THE ORIGINAL PROBLEM 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 C INITIAL GUESS FOR THE EXTREMAL FREQUENCIES--EQUALLY C SPACED ALONG THE GRID 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 C CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION C PROBLEM C WRITE(*,850) 850 FORMAT(//' STARTING REMEZ ITERATIONS') CALL REMEZ C C CALCULATE THE IMPULSE RESPONSE. C WRITE(*,'(/A)') ' CALCULATING IMPULSE RESPONSE' 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 PROGRAM OUTPUT SECTION. C 350 WRITE(IOUT,360) 360 FORMAT( 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(22X,15HBANDPASS FILTER/) IF(JTYPE.EQ.2) WRITE(IOUT,370) 370 FORMAT(22X,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,'(15X,''SAMPLING FREQUENCY = '',F10.2)') SFREQ WRITE(IOUT,'(20X,''GRID DENSITY = '',I3/)') LGRID 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 9020 J=1,JB 9020 EDGE(J)=EDGE(J)*SFREQ 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(UDAR.NE.'Y'.AND.UDAR.NE.'y') THEN 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) ENDIF IF(UDWF.NE.'Y'.AND.UDWF.NE.'y') THEN 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) ENDIF 450 CONTINUE DO 452 J=1,NZ IX=IEXT(J) 452 GRID(J)=GRID(IX)*SFREQ WRITE(IOUT,455) (GRID(J),J=1,NZ) 455 FORMAT(/2X,47HEXTREMAL FREQUENCIES--MAXIMA OF THE ERROR CURVE/ 1 (2X,5F15.7)) WRITE(IOUT,460) 460 FORMAT(/1X,70(1H*)/) C CALCULATE FREQUENCY RESPONSE WRITE(*,710) 710 FORMAT(/' CALCULATING FREQUENCY RESPONSE ') WRITE(IOUT,830) 830 FORMAT(/15X'*********FREQUENCY RESPONSE*********'/) WRITE(IOUT, 711) 711 FORMAT(/15X,' F ',14X,'ABSOLUTE',15X,' DB') NPOINT=10*NFILT+1 NPOINT=MAX(201,NPOINT) DELW=.5/(NPOINT-1.) NPOINT = 1 + (FSTOP - FSTART)/DELW DO 610 IKA=1,NPOINT OMEGA(IKA)=(IKA-1)*DELW + FSTART SUMAR=0. SUMAC=0. SGN=(-1.)**NEG PI2S=SNGL(PI2) DO 620 NIK=1,NFCNS SUMAR=SUMAR+H(NIK)*(COS(PI2S*OMEGA(IKA)*(NIK-1))+SGN* 1 COS(PI2S*OMEGA(IKA)*(NFILT-NIK))) SUMAC=SUMAC+H(NIK)*(SIN(PI2S*OMEGA(IKA)*(NIK-1))+SGN* 1 SIN(PI2S*OMEGA(IKA)*(NFILT-NIK))) 620 CONTINUE IF(NODD.EQ.1.AND.NEG.EQ.0) SUMAR=SUMAR-H(NFCNS)*COS(PI2S* 1 OMEGA(IKA)*(NFCNS-1)) IF(NODD.EQ.1.AND.NEG.EQ.0) SUMAC=SUMAC-H(NFCNS)*SIN(PI2S* 1 OMEGA(IKA)*(NFCNS-1)) RESPA(IKA)=SQRT(SUMAR**2+SUMAC**2) RESPA1=RESPA(IKA) IF(RESPA1.LE.1.E-15) RESPA(IKA)=1.E-15 RESPA(IKA)=20.*LOG10(RESPA(IKA)) OMEGA(IKA)=OMEGA(IKA)*SFREQ IF(MOD(IKA,2).EQ.0) THEN WRITE(IOUT,630) OMEGA(IKA),RESPA1,RESPA(IKA) ENDIF IF(LINDB.EQ.'L'.OR.LINDB.EQ.'l') RESPA(IKA)=RESPA1 610 CONTINUE 630 FORMAT(F20.3,8X,F15.9,8X,F15.9) WRITE(*,649) 649 FORMAT(/1X,'PLOTTING THE FREQUENCY RESPONSE') CALL PLOT(OMEGA,RESPA,NPOINT,IOUT,RMAX,LINDB) OPEN(2,FILE=COEFL,STATUS='NEW') WRITE(2,500) NFILT 500 FORMAT(' ',I3) WRITE(2,510) (H(J),J=1,NFCNS) 510 FORMAT(' ',E22.8) IF(NODD.EQ.1) GOTO 540 550 DO 520 I=1,NFCNS II=NFCNS-I+1 IF(NEG.EQ.0) WRITE(2,510) H(II) H(II)= -H(II) 520 IF(NEG.EQ.1) WRITE(2,510) H(II) GOTO 530 540 JJ=(NFILT+1)/2 IF(NFCNS.EQ.JJ) GOTO 560 ZERO = 0. WRITE(2,510) ZERO GOTO 550 560 NFCNS=NFCNS-1 GOTO 550 530 CONTINUE WRITE(*,'(/A,\)') ' CREATE (FREQ,RESPONSE) FILE (Y OR N)? ' READ(*,'(A)') KPLT IF(KPLT.NE.'Y'.AND.KPLT.NE.'y') STOP 5530 WRITE(*,'(A,\)') ' ENTER FILENAME: ' READ(*,'(A)') COEFL INQUIRE(FILE=COEFL,EXIST=IXST) IF(IXST) THEN WRITE(*,*) CHAR(7) WRITE(*,'(A)') ' FILE ALREADY EXISTS. USE ANOTHER NAME.' GOTO 5530 ENDIF OPEN(15,FILE=COEFL,STATUS='NEW') DO 5531 IKA=1,NPOINT 5531 WRITE(15,'(1X,F20.3,2X,F20.4)') OMEGA(IKA),RESPA(IKA) END C C SUBROUTINE PLOT(OMEGA,RESPA,NPOINT,IOUT,RMAX,LINDB) C SUBROUTINE PLOT(OMEGA,RESPA,NPOINT,IOUT,RMAX,LINDB) DIMENSION OMEGA(4001),RESPA(4001) CHARACTER LINE(60),BLINE(60),SYMBOL,DASH(60),LINDB DATA BLINE/'I',58*' ','I'/ DATA DASH/'I',58*'-','I'/ RMAX = -1.E20 IF(LINDB.EQ.'D'.OR.LINDB.EQ.'d') THEN WRITE(IOUT,22) 22 FORMAT(/' F ',60X,' DB') ELSE WRITE(IOUT,222) 222 FORMAT(/' F ',60X,' GAIN') ENDIF DO 1 I=1,NPOINT 1 IF(RESPA(I).GE.RMAX) RMAX = RESPA(I) RMIN = -80. IF(LINDB.EQ.'L'.OR.LINDB.EQ.'l') RMIN = 0. DO 2 K=1,NPOINT,2 DO 20 JJ=1,60 20 LINE(JJ)=BLINE(JJ) IF(K.EQ.1.OR.K.EQ.NPOINT.OR.MOD(K-1,10).EQ.0) THEN DO 21 JJ=1,60 21 LINE(JJ)=DASH(JJ) ENDIF IF(RESPA(K).LT.RMIN) THEN DB= RMIN SYMBOL='X' ELSE DB=RESPA(K) SYMBOL='*' ENDIF IND=MAX(NINT(60*(DB-RMIN)/(RMAX-RMIN)),1) LINE(IND)=SYMBOL 2 WRITE(IOUT,3) OMEGA(K),LINE,RESPA(K) 3 FORMAT(F7.2,1X,60A1,1X,F10.4) RETURN END C C----------------------------------------------------------------------- C FUNCTION: EFF C FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE C AS A FUNCTION OF FREQUENCY. C AN ARBITRARY FUNCTION OF FREQUENCY CAN BE C APPROXIMATED IF THE USER REPLACES THIS FUNCTION C WITH THE APPROPRIATE CODE TO EVALUATE THE IDEAL C MAGNITUDE. NOTE THAT THE PARAMETER FREQ = f/SFREQ = C NORMALIZED FREQUENCY. C----------------------------------------------------------------------- C FUNCTION EFF(FREQ,FX,WTX,LBAND,JTYPE,SFREQ) DIMENSION FX(10),WTX(10) DOUBLE PRECISION SFREQ IF(JTYPE.EQ.2) GO TO 1 EFF=FX(LBAND) RETURN 1 EFF=FX(LBAND)*FREQ*SFREQ RETURN END C C----------------------------------------------------------------------- C FUNCTION: WATE C FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION C OF FREQUENCY. SIMILAR TO THE FUNCTION EFF, THIS FUNCTION CAN C BE REPLACED BY A USER-WRITTEN ROUTINE TO CALCULATE ANY C DESIRED WEIGHTING FUNCTION. C----------------------------------------------------------------------- C FUNCTION WATE(FREQ,FX,WTX,LBAND,JTYPE,SFREQ) DIMENSION FX(10),WTX(10) DOUBLE PRECISION SFREQ IF(JTYPE.EQ.1) THEN WATE=WTX(LBAND) ELSE IF(JTYPE.EQ.2) THEN IF (FREQ.LE..00001) THEN WATE=WTX(LBAND)/(.00001*SFREQ) ELSE WATE=WTX(LBAND)/(FREQ*SFREQ) ENDIF ELSE IF(JTYPE.EQ.3) THEN WATE = -WTX(LBAND) ENDIF RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: ERROR C THIS ROUTINE WRITES AN ERROR MESSAGE IF AN C ERROR HAS BEEN DETECTED IN THE INPUT DATA. C----------------------------------------------------------------------- C SUBROUTINE ERROR COMMON /OOPS/NITER,IOUT WRITE(IOUT,1) 1 FORMAT(44H ************ ERROR IN INPUT DATA **********) RETURN END 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 SUBROUTINE REMEZ COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID COMMON /OOPS/NITER,IOUT DIMENSION IEXT(202),AD(202),ALPHA(202),X(202),Y(202) DIMENSION DES(3300),GRID(3300),WT(3300) DIMENSION A(202),P(202),Q(202) DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q DOUBLE PRECISION DK,DAK DOUBLE PRECISION AD,DEV,X,Y,GRID DOUBLE PRECISION GEE,D DOUBLE PRECISION DES C C THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25 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(*,131)DEV 131 FORMAT(1X,12HDEVIATION = ,E20.8) 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 C SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST C APPROXIMATION 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) 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) 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) 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) 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) 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) 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) 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=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) 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 C C----------------------------------------------------------------------- C FUNCTION: D C FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION C COEFFICIENTS FOR USE IN THE FUNCTION GEE. C----------------------------------------------------------------------- C DOUBLE PRECISION FUNCTION D(K,N,M) COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID DIMENSION IEXT(202),AD(202),ALPHA(202),X(202),Y(202) DIMENSION DES(3300),GRID(3300),WT(3300) DOUBLE PRECISION AD,DEV,X,Y,GRID DOUBLE PRECISION Q DOUBLE PRECISION PI2 DOUBLE PRECISION DES D=1.0D0 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.0D0/D RETURN END C C----------------------------------------------------------------------- C FUNCTION: GEE C FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE C LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM C----------------------------------------------------------------------- C DOUBLE PRECISION FUNCTION GEE(K,N) COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID DIMENSION IEXT(202),AD(202),ALPHA(202),X(202),Y(202) DIMENSION DES(3300),GRID(3300),WT(3300) DOUBLE PRECISION P,C,D,XF DOUBLE PRECISION PI2 DOUBLE PRECISION AD,DEV,X,Y,GRID DOUBLE PRECISION DES P=0.0 XF=GRID(K) XF=DCOS(PI2*XF) D=0.0 DO 1 J=1,N C=XF-X(J) C=AD(J)/C D=D+C 1 P=P+C*Y(J) GEE=P/D RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: OUCH C WRITES AN ERROR MESSAGE WHEN THE ALGORITHM FAILS TO C CONVERGE. THERE SEEM TO BE TWO CONDITIONS UNDER WHICH C THE ALGORITHM FAILS TO CONVERGE: (1) THE INITIAL C GUESS FOR THE EXTREMAL FREQUENCIES IS SO POOR THAT C THE EXCHANGE ITERATION CANNOT GET STARTED, OR C (2) NEAR THE TERMINATION OF A CORRECT DESIGN, C THE DEVIATION DECREASES DUE TO ROUNDING ERRORS C AND THE PROGRAM STOPS. IN THIS LATTER CASE THE C FILTER DESIGN IS PROBABLY ACCEPTABLE, BUT SHOULD C BE CHECKED BY COMPUTING A FREQUENCY RESPONSE. C----------------------------------------------------------------------- C SUBROUTINE OUCH COMMON /OOPS/NITER,IOUT WRITE(IOUT,1)NITER 1 FORMAT(44H ************ FAILURE TO CONVERGE **********/ 1 ' PROBABLE CAUSE IS MACHINE ROUNDING ERROR'/ 2 ' NUMBER OF ITERATIONS =',I4/ 3 ' IF THE NUMBER OF ITERATIONS EXCEEDS 3,'/ 4 ' THE DESIGN MAY BE CORRECT. CHECK FREQUENCY RESPONSE.') RETURN END