<< >>


() - , , , , .

. .

1. 24- 0,08, 0,16 1,0.

: , , . ,

. A.3.1. 24-

.

1. 24, 1, 2, 0, 16

2. 0, 0.08, 0.16, 0.5

3. 1, 0

4. 1, 1

1 ( ) , , , ; 2 ; 3 ; 4 . . .3.1 ( Honeywell 6000), . .3.2 . . .3.1 , (-56,4 ).

. .3.2. ( ) , . .3.1.

2. 32- 0,1 0,425, 0,2 0,35 , 10 1 .

1. 32, 1, 3, 0, 16

2. 0, 0.1, 0.2, 0.35, 0.425, 0.5

3. 0, 1, 0

4. 10, 1, 10

. .3.3 ; . 3.80.

. A.3.3. 32-

.

3. 32- 0,5.

1. 32, 2, 1, 0, 16

2. 0, 0.5

3. 1

4. 1

. .3.4 , . 3.68.

. A.3.4. 32- .

4. 20- 0,05 0,5.

1. 20, 3, 1, 0, 16

2. 0.05, 0.5

3. 1

4. 1

. A.3.5. 20-

.

. . 3.5 , . .3.6.

. A.3.6. , . .3.5.

 

-

 

Ѡ PROGRAM FOP. THE DESIGN OF LINEAR PHASE FINITE IMPULSE

Ѡ RESPONSE (FIR) FILTERS USING THE REMEZ EXCHANGE ALGORITHM

Ѡ JIM MCCLELLAN, RICE UNIVERSITY, APRIL 13, 1973

C THREE TYPES OF FILTERS ARE INCLUDEDBANDPASS FILTFRS

C DIFFERENTIATORS, AND HILBERT TRANSFORM FILTERS C

C THE INPUT DATA CONSISTS OF 5 CARDS C

C CARD 1FILTER LENGTH, TYPE OF FILTER. 1-NULTIPLE

C PASSBAND/STOPBAND, 2-DIFFERENTIATOR, 3-HILBERT TRANSFORM

C FILTER. NUMBER OF BANDS, CARD PUNCH DESIRED, AND GRID

C DENSITY.

C

C CARD 2BANDEDGES, LOWER AND UPPER EDGES FOR EACH BAND

C WITH A MAXIMUM OF 10 BANDS.

C CARD 3DESIRED FUNCTION (OR DESIRED SLOPE IF A

Ѡ DIFFERENTIATOR) FOR EACH BAND.

C

C CARD 4WEIGHT FUNCTION IN EACH BAND. FOR A

C DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY

C PROPORTIONAL TO F.

C

C THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 BANDPASS

C FILTER WITH STOPBANDS 0 TO 0.1 AND 0.425 TO 0.5, AND

C PASSBAND FROM 0.2 TO 0.35 WITH WEIGHTING OF 10 IN THE

C STOPBANDS AND 1 IN THE PASSBAND. THE IMPULSE RESPONSE

C WILL BE PUNCHED AND THE GRID DENSITY IS 32. C

Ѡ SAMPLE INPUT DATA SETUP

C 32,1,3,1,32

C 0,0.1,0.2,0.35,0.425,0.5

Ѡ 0,1,0

Ѡ 10,1,10

Ѡ THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 WIDEBAND

Ѡ DIFFERENTIATOR WITH SLOPE 1 AND WEIGHTING OF 1/F. THE

Ѡ IMPULSE RESPONSE WILL NOT BE PUNCHED AND THE GPID

Ѡ DENSITY IS ASSUMED TO BE 16.

Ѡ 32,2,1,0,0

Ѡ 0,0.5

Ѡ 1.0

Ѡ 1.0

 

COMMO NPI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID

DIMENSION I EXT (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

PI2=6.2831853 07179586

PI=3.141592653589793

Ѡ THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128, BUT

Ѡ THIS UPPER LIMIT CAN BE CHANGED BY REDIMENSIONING THE

Ѡ ARRAYS I EXT, AD, ALPHA, X, Y, H TO BE NFMAX/2 ♦ 2.

Ѡ THE ARRAYS DES, GRID, AND WT MUST DIMENSIONED

Ѡ 16(NFMAX/2 ♦ 2) .

NFMAX=128

100 CONTINUE

JTYPE=0

Ѡ PROGRAM INPUT SECTION

READ *,NFILT,JTYPE,NBANDS,JPUNCH,LGRID IF(NFILT.GT.NFMAX.OR,NFILT.LT.3) CALL ERROR

IF(NBANDS.LE.0) NBANDS=1

Ѡ GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED

Ѡ OTHERWISE

IF(LGRID.LE.O) LGRID=16

JB=2*NBANDS

READ , (EDGE(J) ,J=1,JB)

READ ,(FX(J) ,J = 1, NBANDS)

READ ♦ , (WTX(J),J=1,NBANDS)

IF(JTYPE.EQ.0) CALL ERROR

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

Ѡ SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID

Ѡ IS (FILTER LENGTH ♦ 1)♦GRID DENSITY/2

GRID (1) =EDGE (1)

DELF=LGRID* NFCNS

DELF=0.5/DELF

IF(NEG.EQ.O) GO TO 135

IF(EDGE(1).LT.DELF) GRID(1)= DELF 135

135 CONTINUE

J=1

L= 1

LBAND=1 14 0

140 FUP=EDGE(L*1) 145

145 TEMP=GRID(J)

Ѡ CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT

Ѡ FUNCTION ON THE GRID

DES(J) =EFF(TEMP,FX,WTX,LBAND,JTYPE) WT(J)=WATE(TEMP,FX,WTX,LBAND,JTYPE)

GRID(J)=TEMP*DELF

IF (GRID (J) .GT.FUP)

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

160 NGRID=J-1

IF(NEG.NE.NODD) GO TO 165

I F (GRID (NGRID) .GT. (0.5-DELF) ) NGRI D=NGRID-1

Ѡ SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT

Ѡ TO THE ORIGINAL PROBLEM

IF(NEG) 170,170,180

170 IF(NODD.EQ.1)

GOTO 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(PIGRID(J)) DES (J) *DES (J) /CHANGE

185 WT (J) =WT (J) CHANGE

GO TO 200

190 DO 195 J*1,NGRID

CHANGE=DSIN(PI2GRID(J))

DES (J) =DES (J) /CHANGE

195 WT(J)=WT(J)CHANGE

Ѡ INITIAL GUESS FOR THE EXTREMAL FREQUENCIESEQUALLY

Ѡ SPACED ALONG THE GRID

200 TEMP=FLOAT(NGRID-1)/FLOAT(NFCNS)

DO 210 J=1,NFCNS

210 IEXT(J) = (J-1) TEMP*1

IEXT(NFCNS*1)=NGRID

NM1=NFCNS-1

NZ-NFCNS*1

Ѡ CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION

Ѡ PROBLEM

C ALL REMEZ(EDGE,NBANDS)

Ѡ CALCULATE THE IMPULSE RESPONSE.

IF(NEG) 300,300,320

300 IF(NODD.EQ.O) GO TO 310

DO 305 J=1,NM1

305 H (J) =0.5*ALPHA(NZ-J)

H(NFCNS)-ALPHA(1) GO TO 350

310 H(1)=0.25*ALPHA (NFCNS)

DO 315 J=2,NM1

315 H(J)=0.25*(ALPHA(NZ-J)♦ALPHA(NFCNS*2-J))

H(NFCNS) =0.5ALPHA(1) ♦0.25ALPHA(2)

GO TO 350

320 IF(NODD.EQ.0) GO TO 330

H(1)~Q.25ALPHA(NFCNS)

H(2)=0.2SALPHA (NM1)

DO 325 J=3,NM1

325 H(J)=0.25(ALPHA(NZ-J) -ALPHA(NFCNS*3-J)) H(NFCNS)=0.5ALPHA(1)-0.25ALPHA(3)

H (NZ) =0.0

GO TO 350

330 H(1)=0.25ALPHA(NFCNS)

DO 335 J=2,NM1

335 (J)=0.25*(ALPHA(NZ-J)-ALPHA(NFCNS+2-J))

H(NFCNS)=0.5*ALPHA(1)-0.25*ALPHA(2)

PROGRAM OUTPUT SECTION.

350 PRINT 360

360 FORMAT(1H1, 7 0(1H*)//25X,'FINITE IMPULSE RESPONSE (FIR)' 125X,"LINEAR PHASE DIGITAL FILTER DESIGN*/

225X,'REMEZ EXCHANGE ALGORITHM'/) IF(JTYPE.EQ.1) PRINT 365

365 FORMAT(25X,'BANDPASS FILTER'/)

IF(JTYPE.EQ.2) PRINT 370

370 FORMAT(25X,'DIFFERENTIATOR'/)

IF (JTYPE.EQ.3) PRINT 375

375 FORMAT(25X,'HILEERT TRANSFORMER'/)

PRINT 378 ,NFILT

378 FORMAT(15X,*FILTER LENGTH = ',13/)

PRINT 350

380 FORMAT(15X,* ***** IMPULSE RESPONSE ***')

DO 381 J=1,NFCNS

K=NFILT+1-J

IF(NEG.EQ.O) PRINT 382,J,H(J),K

IF(NEG.EQ.I) PRINT 383,J,H(J),K

385 CONTINUE

382 FORMAT(20X,'H(',13,*) = '.E15.8,* = H(*,I4,*)*)

383 FORMAT (20X,'H(*,I3,*) = *,E15.8," = -H(*,I4,')*) IF(NEG.EQ.1.AND.NODD.EQ.1) PRINT 384,NZ

384 FORMAT(20X,'H(*,13,') = 0.0') DO 450 K=1,NBANDS,4

KUP=K+3

IF(KUP.GT.NBANDS) KUP=NBANDS PRINT 385,(J,J=K,KUP)

385 FORMAT(/24X,4('BAND*,I 3,8X) )

PRINT 390, (EDGE (2*J-1) ,J=K,KUP)

390 FORMAT(2X,'LOWER BAND EDGE',5F15.9)

PRINT 395,(EDGE(2*J),J=K,KUP)

395 FORMAT(2X,'UPPER BAND EDGE',5F15.9)

IF(JTYPE.NE.2) PRINT

400, (FX (J) , J=K, KUP)

FORMAT(2X,'DESIRED VALUE*,2X,5F15.9)

IF(JTYPE.EQ.2) PRINT 405,(FX(J),J=K,KUP)

405 FORMAT(2X,*DESIRED SLOPE',2X,5F15.9)

PRINT

410, (WTX(J) , J=K,KUP)

410 FORMAT(2X,* WEIGHTING*,6X,5F15.9)

DO 420 J=K,KUP

420 DEVIAT (J) =DEV/WTX (J)

PRINT 425,(DEVIAT(J),J=K,KUP)

425 FORMAT(2X,'DEVIATION',6X,5F15.9)

IF(JTYPE.NE.1)

GO TO 450 DO

430 J=K,KUP 430 DEVIAT(J)=20.0*ALOG10(DEVIAT(J) )

PRINT 435, (DEVIAT(J),J=K, KUP)

435 FORMAT(2X,'DEVIATION IN DB',5F15.9)

450 CONTINUE

PRINT 455,(GRID (IEXT(J)),J=1,NZ)

455 FORMAT(/2X,'EXTREMAL FREQUENCIES'/(2X,5F12.7))

PRINT 460

460 FORMAT(/1X,70 (1H*)/1 Hi)

IF (JPUNCH. NE. 0) PUNCH * , (H (J) , J= 1,NFCNS)

IF(NFILT.NE.O) GO TO 100

RETURN

END

 

FUNCTION EFF(TEMP,FX,WTX,LBAND,JTYPE)

FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE

AS A FUNCTION OF FREQUENCY.

DIMENSION FX(5),WTX(5)

IF(JTYPE.EQ.2) GO TO

EFF=FX (LBAND)

RETURN 1 EFF=FX (LBAND) TEMP

RETURN

END

 

 

FUNCTION WATE (TEMP,FX,WTX,LbAND,JTYPE)

FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION

OF FREQUENCY.

DIMENSION FX(5),WTX(5)

IF(JTYPE.EQ.2) GO TO 1

WATE=WTX (LBAND)

RETURN

IF(FX(LBAND).LT.0.0001) GO TO 2

WATE=WTX (LBAND) /TEMP

RETURN

WATE=WTX(LBAND)

RETURN

END

 

SUBROUTINE ERROR PRINT 1

1 FORMAT(' ************ ERROR IN INPUT DATA **********')

STOP

END

 

 

SUBROUTINE REMEZ(EDGE,NBANDS)

THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM FOR THE WEIGHTED CHEBYCHEV APPROXIMATION OF A CONTINUOUS

FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE

ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE

DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE

GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE

EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYCHEV

ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL

FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES THE COEFFICIENTS OF THE EEST APPROXIMATION.

COMMON PI2,AD,DEV,X,Y,GRID,DES,WT,ALPHA,IEXT,NFCNS,NGRID DIMENSION EDGE(20)

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 P12,DNUM,DDEN,DTEMP , A , P , Q

DOUBLE PRECISION AD,DEVrX,Y

THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25

ITRMAX=25

DEVL=-1.0

NZ=NFCNS*1

NZZ=NFCNS*2

NITER=0

100 CONTINUE

lEXT(NZZ) =NGRID1

NITER=NITER*1

IF(NITER.GT.ITRMAX) GO TO 400

DO 110 J=1,NZ DTEMP=GRID(IEXT(J))

DTEMP= DCOS (DTEMP*PI 2)

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=K*AD (J) /WT (L)

DDEN= DDEN♦DTEMP

130 K=-K

DEV=DNUM/DDEN NU=1

IF(DEV.GT.O.O) Nl)=-1

DEV=-NU*DEV

K=NU

DO 140 J=1,NZ L=IEXT(J)

DTEMP=K*DEV/WT(L) Y (J) =DES (L) ♦DTEMP

140 K=-K

IF(DEV.GE,DEVI) 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

SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST

APPROXIMATION

200 IF(J.EQ.NZZ) YNZ=COMP

IF(J.GE.NZZ) GO TO 300

KUP'IEXT (J+l)

L-IEXT (J) ♦ 1

NUT=-NUT

IF (J. EQ. 2) Y 1=C0MP C

OMP=DEV

IF(L.GE.KUP) GO TO 220

ERR=GEE(L,NZ)

ERR=(ERR-DES(L))*WT(L)

DTEMP=NUT* ERR-COMP

IF(DTEMP.LE.0.0) GO TO 220 COM P= NUT* ERR

210 L=L*1

IF(L.GE.KUP) GO TO 215

ERR=GEE (L, NZ)

ERR=(ERR-DES(L))*WT(L)

DTEMP=NUT*ERR-COMP

IF(DTEMP.LE.0.0) GO TO 215

COMP=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= NUT* ERR-COMP

IF(DTEMP.GT.0.0) GO TO 230

IF(JCHNGE.LE.0) GO TO 225

GO TO 2 60 230

COMP=NUT*ERR

235 L=L-1

IF(L.LE.KLOW) GO TO 240

ERR=GEE(L,NZ)

ERR=(ERR-DES(L))*WT(L)

DTEMP= NUT* ERR-COMP

IF(DTEMP.LE.0.0) GO TO 240

COMP=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-NUT*EPR-COMP

IF(DTEMP.LE.0.0) GO TO 255

COMP=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 (N2) )

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=NUT*ERR-COMP

IF(DTEMP.LE.0.0) GO TO 310

COMP=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=NUT*ERR-COMP

IF (DTEMP.LF.0.0) GO TO 330

J=NZZ

COMP=NUT* ERR

LUCK=LUCK*10 GO TO 235

340 IF(LUCK.EQ.6) GO TO 370

DO 345 J=1,NFCNS

345 I EXT (NZZ-J) = I EXT (NZ-J)

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

 

CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION

USING THE INVERSE DISCRETE FOURIER TRANSFORM

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(EDGE(1).EQ.0.0.AND.EDGE(2*NBANDS).EQ.0.5) KKK=1

IF(NFCNS.LE.3) =1

IF(.EQ.1) QO 405

DTEMP=DCOS(PI2GFID(1))

DNUM=DCOS(PI2GRID(NGRID))

AA=2.0/(DTEMP-DNOM)

BB=-(DTEMP*DNUM)/(DTEMP-DNUM)

405 CONTINUE

DO 430 Js 1,NFCNS

FT=(J-1)DELE

XT=DCOS(PI2 FT)

IF(KKK.EQ.1) GO TO 410

XT= (XT-PB)/

FT=ARCOS(XT)/PI 2

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

4 20 IF ((XT-XE) .LT.FSH) GO TO 415

GRID(1) =FT A (J) =GEF (1,N2)

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)DDEN

IF(NM1.LT.1) GO TO 505

DO 500 K=1,NM1

500 DTEMP=DTEMP*A(K*1)DCOS(DNUMK)

505 DTEMP=2.0DTEMP*A(1)

510 ALPHA(J)sDTEMP

DO 550 J2,NFCNS

550 ALPHA(J)=2ALPHA(J)/CN

ALPHA(1)=ALPHA(1)/CN IF(KKK.EQ.1) GO TO 545 P(1)=2.0ALPHA(NFCNS)BB*ALPHA(NM1)

P(2)=2.0AAALPHA(NFCNS)

Q(1)=ALPHA(NFCNS-2)-ALPHA(NFCNS)

DO 540 J=2,NM1

IF(J.LT.NMI) GO TO 515

AA=0.5AA BB=0.5BB

515 CONTINUE

P (J*1)=0.0

DO 520 K=1,

J A(K)=P(K) 520

P(K) =2.0BBA(K)

P (2) =P (2) *A(1) 2.0AA JM1=J-1

DO 525 K=1,JM1

525 P(K)=P(K) *Q(K) *AAA(K*1)

JP1=J*1

DO 530 K=3,JP1

530 P(K)=P(K) *AAA (K-1)

IF(J.EQ.NMI) GO TO 540

DO 535 K=1,J

535 Q(K)=-A(K)

Q(1)=Q(1) *ALPhA (NFCNS-1-J)

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

 

DOUBLE PRECISION FUNCTION D(K,N,M)

FUNCTION TO CALCULATE THE LAGRANGE INTERPOIATION

COEFFICIENTS FOR USE IN THE FUNCTION GEE.

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.0D*(Q-X(J))

2)      CONTINUE

3)      CONTINUE D=1.0/D

RETURN

END

 

DOUBLE PRECISION FUNCTION GEE(K,N)

FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE

LAGRANGE INTERPOIATION FORMULA IN THE BARYCENTRIC FORM

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

P=P*CY(J)

GEE=P/D RETURN

END

SUBROUTINE OUCH PRINT 1

1 FORMAT('**************FAILURE TO CONVERGE*************/

OPROBABLE CAUSE IS MACHINE ROUNDING ERROR*/

2'0THE IMPULSE RESPONSE MAY BE CORRECT"/

3'OCHECK WITH A FREQUENCY RESPONSE*)

RETURN

END

 

 



<< >>