FIR Linear Phase Filter Design Program
Return to Main Software Page
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 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
C 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(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION H(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DIMENSION EDGE(20),FX(10),WTX(10),DEVIAT(10)
DOUBLE PRECISION PI2,PI
DOUBLE PRECISION AD,DEV,X,Y
DOUBLE PRECISION GEE,D
INTEGER BD1,BD2,BD3,BD4
DATA BD1,BD2,BD3,BD4/1HB,1HA,1HN,1HD/
INPUT=I1MACH(1)
IOUT=I1MACH(2)
PI=4.0*DATAN(1.0D0)
PI2=2.0D00*PI
C
C THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128, 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=128
100 CONTINUE
JTYPE=0
C
C PROGRAM INPUT SECTION
C
READ(INPUT,110) NFILT,JTYPE,NBANDS,LGRID
IF(NFILT.EQ.0)STOP
110 FORMAT(4I5)
IF(NFILT.LE.NFMAX.AND.NFILT.GE.3) GO TO 115
CALL ERROR
STOP
115 IF(NBANDS.LE.0) NBANDS=1
C
C GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED
C OTHERWISE
C
IF(LGRID.LE.0) LGRID=16
JB=2*NBANDS
READ(INPUT,120) (EDGE(J),J=1,JB)
120 FORMAT(4F15.9)
READ(INPUT,120) (FX(J),J=1,NBANDS)
READ(INPUT,120) (WTX(J),J=1,NBANDS)
IF(JTYPE.GT.0.AND.JTYPE.LE.3) GO TO 125
CALL ERROR
STOP
125 NEG=1
IF(JTYPE.EQ.1) NEG=0
NODD=NFILT/2
NODD=NFILT-2*NODD
NFCNS=NFILT/2
IF(NODD.EQ.1.AND.NEG.EQ.0) NFCNS=NFCNS+1
C
C SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID
C IS (FILTER LENGTH + 1)*GRID DENSITY/2
C
GRID(1)=EDGE(1)
DELF=LGRID*NFCNS
DELF=0.5/DELF
IF(NEG.EQ.0) GO TO 135
IF(EDGE(1).LT.DELF) GRID(1)=DELF
135 CONTINUE
J=1
L=1
LBAND=1
140 FUP=EDGE(L+1)
145 TEMP=GRID(J)
C
C CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
C FUNCTION ON THE GRID
C
DES(J)=EFF(TEMP,FX,WTX,LBAND,JTYPE)
WT(J)=WATE(TEMP,FX,WTX,LBAND,JTYPE)
J=J+1
GRID(J)=TEMP+DELF
IF(GRID(J).GT.FUP) GO TO 150
GO TO 145
150 GRID(J-1)=FUP
DES(J-1)=EFF(FUP,FX,WTX,LBAND,JTYPE)
WT(J-1)=WATE(FUP,FX,WTX,LBAND,JTYPE)
LBAND=LBAND+1
L=L+2
IF(LBAND.GT.NBANDS) GO TO 160
GRID(J)=EDGE(L)
GO TO 140
160 NGRID=J-1
IF(NEG.NE.NODD) GO TO 165
IF(GRID(NGRID).GT.(0.5-DELF)) NGRID=NGRID-1
165 CONTINUE
C
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
CALL REMEZ
C
C CALCULATE THE IMPULSE RESPONSE.
C
IF(NEG) 300,300,320
300 IF(NODD.EQ.0) GO TO 310
DO 305 J=1,NM1
NZMJ=NZ-J
305 H(J)=0.5*ALPHA(NZMJ)
H(NFCNS)=ALPHA(1)
GO TO 350
310 H(1)=0.25*ALPHA(NFCNS)
DO 315 J=2,NM1
NZMJ=NZ-J
NF2J=NFCNS+2-J
315 H(J)=0.25*(ALPHA(NZMJ)+ALPHA(NF2J))
H(NFCNS)=0.5*ALPHA(1)+0.25*ALPHA(2)
GO TO 350
320 IF(NODD.EQ.0) GO TO 330
H(1)=0.25*ALPHA(NFCNS)
H(2)=0.25*ALPHA(NM1)
DO 325 J=3,NM1
NZMJ=NZ-J
NF3J=NFCNS+3-J
325 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF3J))
H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(3)
H(NZ)=0.0
GO TO 350
330 H(1)=0.25*ALPHA(NFCNS)
DO 335 J=2,NM1
NZMJ=NZ-J
NF2J=NFCNS+2-J
335 H(J)=0.25*(ALPHA(NZMJ)-ALPHA(NF2J))
H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(2)
C
C PROGRAM OUTPUT SECTION.
C
350 WRITE(IOUT,360)
360 FORMAT(1H1, 70(1H*)//15X,29HFINITE IMPULSE RESPONSE (FIR)/
113X,34HLINEAR PHASE DIGITAL FILTER DESIGN/
217X,24HREMEZ EXCHANGE ALGORITHM/)
IF(JTYPE.EQ.1) WRITE(IOUT,365)
365 FORMAT(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,380)
380 FORMAT(15X,28H***** IMPULSE RESPONSE *****)
DO 381 J=1,NFCNS
K=NFILT+1-J
IF(NEG.EQ.0) WRITE(IOUT,382) J,H(J),K
IF(NEG.EQ.1) WRITE(IOUT,383) J,H(J),K
381 CONTINUE
382 FORMAT(13X,2HH(,I2,4H) = ,E15.8,5H = H(,I3,1H))
383 FORMAT(13X,2HH(,I2,4H) = ,E15.8,6H = -H(,I3,1H))
IF(NEG.EQ.1.AND.NODD.EQ.1) WRITE(IOUT,384) NZ
384 FORMAT(13X,2HH(,I2,8H) = 0.0)
DO 450 K=1,NBANDS,4
KUP=K+3
IF(KUP.GT.NBANDS) KUP=NBANDS
WRITE(IOUT,385) (BD1,BD2,BD3,BD4,J,J=K,KUP)
385 FORMAT(/24X,4(4A1,I3,7X))
WRITE(IOUT,390) (EDGE(2*J-1),J=K,KUP)
390 FORMAT(2X,15HLOWER BAND EDGE,5F14.7)
WRITE(IOUT,395) (EDGE(2*J),J=K,KUP)
395 FORMAT(2X,15HUPPER BAND EDGE,5F14.7)
IF(JTYPE.NE.2) WRITE(IOUT,400) (FX(J),J=K,KUP)
400 FORMAT(2X,13HDESIRED VALUE,2X,5F14.7)
IF(JTYPE.EQ.2) WRITE(IOUT,405) (FX(J),J=K,KUP)
405 FORMAT(2X,13HDESIRED SLOPE,2X,5F14.7)
WRITE(IOUT,410) (WTX(J),J=K,KUP)
410 FORMAT(2X,9HWEIGHTING,6X,5F14.7)
DO 420 J=K,KUP
420 DEVIAT(J)=DEV/WTX(J)
WRITE(IOUT,425) (DEVIAT(J),J=K,KUP)
425 FORMAT(2X,9HDEVIATION,6X,5F14.7)
IF(JTYPE.NE.1) GO TO 450
DO 430 J=K,KUP
430 DEVIAT(J)=20.0*ALOG10(DEVIAT(J)+FX(J))
WRITE(IOUT,435) (DEVIAT(J),J=K,KUP)
435 FORMAT(2X,15HDEVIATION IN DB,5F14.7)
450 CONTINUE
DO 452 J=1,NZ
IX=IEXT(J)
452 GRID(J)=GRID(IX)
WRITE(IOUT,455) (GRID(J),J=1,NZ)
455 FORMAT(/2X,47HEXTREMAL FREQUENCIES--MAXIMA OF THE ERROR CURVE/
1 (2X,5F12.7))
WRITE(IOUT,460)
460 FORMAT(/1X,70(1H*)/1H1)
GO TO 100
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 IS THE
C VALUE OF NORMALIZED FREQUENCY NEEDED FOR EVALUATION.
C-----------------------------------------------------------------------
C
FUNCTION EFF(FREQ,FX,WTX,LBAND,JTYPE)
DIMENSION FX(5),WTX(5)
IF(JTYPE.EQ.2) GO TO 1
EFF=FX(LBAND)
RETURN
1 EFF=FX(LBAND)*FREQ
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)
DIMENSION FX(5),WTX(5)
IF(JTYPE.EQ.2) GO TO 1
WATE=WTX(LBAND)
RETURN
1 IF(FX(LBAND).LT.0.0001) GO TO 2
WATE=WTX(LBAND)/FREQ
RETURN
2 WATE=WTX(LBAND)
RETURN
END
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(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DIMENSION A(66),P(65),Q(65)
DOUBLE PRECISION PI2,DNUM,DDEN,DTEMP,A,P,Q
DOUBLE PRECISION DK,DAK
DOUBLE PRECISION AD,DEV,X,Y
DOUBLE PRECISION GEE,D
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(IOUT,131) DEV
131 FORMAT(1X,12HDEVIATION = ,F12.9)
NU=1
IF(DEV.GT.0.0) NU=-1
DEV=-FLOAT(NU)*DEV
K=NU
DO 140 J=1,NZ
L=IEXT(J)
DTEMP=FLOAT(K)*DEV/WT(L)
Y(J)=DES(L)+DTEMP
140 K=-K
IF(DEV.GT.DEVL) GO TO 150
CALL OUCH
GO TO 400
150 DEVL=DEV
JCHNGE=0
K1=IEXT(1)
KNZ=IEXT(NZ)
KLOW=0
NUT=-NU
J=1
C
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(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DOUBLE PRECISION AD,DEV,X,Y
DOUBLE PRECISION Q
DOUBLE PRECISION PI2
D=1.0
Q=X(K)
DO 3 L=1,M
DO 2 J=L,N,M
IF(J-K)1,2,1
1 D=2.0*D*(Q-X(J))
2 CONTINUE
3 CONTINUE
D=1.0/D
RETURN
END
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(66),AD(66),ALPHA(66),X(66),Y(66)
DIMENSION DES(1045),GRID(1045),WT(1045)
DOUBLE PRECISION P,C,D,XF
DOUBLE PRECISION PI2
DOUBLE PRECISION AD,DEV,X,Y
P=0.0
XF=GRID(K)
XF=DCOS(PI2*XF)
D=0.0
DO 1 J=1,N
C=XF-X(J)
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 **********/
141H0PROBABLE CAUSE IS MACHINE ROUNDING ERROR/
223H0NUMBER OF ITERATIONS =,I4/
339H0IF THE NUMBER OF ITERATIONS EXCEEDS 3,/
462H0THE DESIGN MAY BE CORRECT, BUT SHOULD BE VERIFIED WITH AN FFT)
RETURN
END
Return to Main Software Page