A Subroutine for Finite Wordlength FIR Filter Design
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: IDEFIX - TEST PROGRAM FOR THE SUBROUTINE IDEFIR
C
C AUTHOR: ULRICH HEUTE
C INSTITUT FUER NACHRICHTENTECHNIK,UNIV.ERLANGEN-NUERNBERG
C D-8520 ERLANGEN, W-GERMANY
C INPUT:
C FORMAT...6F8.6 (1 CARD).
C DEL1....PASSBAND TOLERANCE;
C DEL2....STOPBAND TOLERANCE;
C OMGP1...FIRST PASSBAND EDGE FREQUENCY;
C OMGS1...FIRST STOPBAND EDGE FREQUENCY;
C OMGP2...SECOND PASSBAND EDGE FREQUENCY;
C OMGS2...SECOND STOPBAND EDGE FREQUENCY.
C
C NORMALIZED FREQUENCY OMG = F / SAMPLING RATE.
C IN LOWPASS AND HIGHPASS CASES, THE LAST TWO PARAMETERS ARE ZERO
C IN BANDFILTER CASES, THE LAST TWO PARAMETERS HAVE THE VALUES
C OMGP(S)2 = 0.5 - OMGP(S)1
C IN HILBERT CASES, DEL2 = DEL1 IS VALID;OMGP1 IS THE HILBERT BAND
C EDGE FREQUENCY; THE LAST T H R E E PARAMETERS ARE ZERO.
C
C OUTPUT:
C ALL DATA CONCERNING THE DESIGN ARE DELIVERED TO THE STANDARD
C OUTPUT UNIT.
C-----------------------------------------------------------------------
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
COMMON /CHEK/ NL, NG, KA, KB
C
C THE COMMON-DEFINITIONS HAVE TO BE ADAPTED TO THE REQUIREMENTS OF
C THE SUBROUTINE 'DESIGN' CALLED BY THE SUBROUTINE 'IDEFIR'.
C
DIMENSION H0(256), RH0(128)
C
C SET UP MACHINE CONSTANTS.
C
C IRCHN...CHANNEL ASSIGNED TO INPUT DEVICE;
C IWCHN...CHANNEL ASSIGNED TO OUTPUT DEVICE;
C
IRCHN = I1MACH(1)
IWCHN = I1MACH(2)
C
C INPUT OF THE TOLERANCE SCHEME.
C
READ (IRCHN,9999) DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2
C
C A MAXIMUM OF 20 DESIGN ITERATIONS IS ALLOWED.
C
ITMAX = 20
C
DP = DEL1
OP = OMGP1
OS = OMGS1
C
WRITE (IWCHN,9998)
C
WRITE (IWCHN,9997)
C
CALL IDEFIR(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP, NF00,
* NF, IWCMIN, IWC, IWACT, DEL11, DEL21, EMAX1, EMAX2, H0, RH0,
* ITERA, IERR)
C
C OUTPUT OF FILTER TYPE AND DESIGN PARAMETERS.
C
IF (ITYP.EQ.0) WRITE (IWCHN,9995) DEL1, DEL2, OP, OS
IF (ITYP.EQ.1) WRITE (IWCHN,9994) DEL1, DEL2, OP, OS
IF (ITYP.EQ.2) WRITE (IWCHN,9993) DEL1, DEL2, OP, OS, OMGP2, OMGS2
IF (ITYP.EQ.3) WRITE (IWCHN,9992) DEL1, DEL2, OP, OS, OMGP2, OMGS2
IF (ITYP.EQ.4) WRITE (IWCHN,9991) DP, OP
C
C TEST FOR POSSIBLE ERROR RETURNS.
C
GO TO (40, 10, 20, 30), IERR
C
C OUTPUT FOR THE CASE OF INCONSISTENT SPECIFICATIONS.
C
10 WRITE (IWCHN,9981)
STOP
C
C OUTPUT, IF THE MAXIMUM ALLOWED FILTER ORDER IS EXCEEDED.
C
20 WRITE (IWCHN,9980)
STOP
C
C OUTPUT,IF NO SOLUTION WAS FOUND.
C
30 WRITE (IWCHN,9982) ITERA
40 CONTINUE
C
C OUTPUT OF REALIZATION PARAMETERS.
C
IF (ITYP.NE.0) WRITE (IWCHN,9990) DEL1, DEL2, OMGP1, OMGS1
WRITE (IWCHN,9989) IWCMIN, NF00
WRITE (IWCHN,9996) IWC, IWACT, NF
C
C OUTPUT OF THE REDUCED TOLERANCE SCHEME.
C
WRITE (IWCHN,9988) DEL11, DEL21
C
C OUTPUT OF EXACT AND ROUNDED COEFFICIENTS.
C
N1 = (NF+1)/2
WRITE (IWCHN,9985)
WRITE (IWCHN,9987) (H0(I),I=1,N1)
WRITE (IWCHN,9984)
WRITE (IWCHN,9987) (RH0(I),I=1,N1)
C
C OUTPUT OF ACTUAL LOWPASS TOLERANCES.
C
WRITE (IWCHN,9986) EMAX1, EMAX2
WRITE (IWCHN,9983) ITERA
C
STOP
C
9999 FORMAT (6F8.6)
C
9998 FORMAT (1H1)
9997 FORMAT (54H ITERATIVE DESIGN OF LINEAR PHASE FIR FILTERS IN DIREC,
* 6HT FORM/36H WITH MINIMUM COEFFICIENT WORDLENGTH//)
9996 FORMAT (41H COEFFICIENT WORDLENGTH APPLIED WC = , I2/6H ACTUA,
* 48HL WORDLENGTH WITHOUT LEADING ZEROS WC ACT = ,
* I2/31H ACTUAL FILTER LENGTH NF = , I3/)
9995 FORMAT (24H LOWPASS WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
* F8.6/14X, 10HOMEGA P = , F8.6, 11H,OMEGA S = , F8.6/)
9994 FORMAT (25H HIGHPASS WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
* F8.6/15X, 10HOMEGA P = , F8.6, 11H,OMEGA S = , F8.6/)
9993 FORMAT (25H BANDPASS WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
* F8.6/15X, 10HOMEGAP1 = , F8.6, 11H,OMEGAS1 = , F8.6/15X,
* 10HOMEGAP2 = , F8.6, 11H,OMEGAS2 = , F8.6/)
9992 FORMAT (25H BANDSTOP WITH DELTA 1 = , F8.6, 11H,DELTA 2 = ,
* F8.6/15X, 10HOMEGAP1 = , F8.6, 11H,OMEGAS1 = , F8.6/15X,
* 10HOMEGAP2 = , F8.6, 11H,OMEGAS2 = , F8.6/)
9991 FORMAT (34H HILBERT-TRANSFORMER WITH DELTA = , F8.6, 9H,OMEGAP =,
* 1H , F8.6/)
9990 FORMAT (35H EQUIVALENT LOWPASS WITH DELTA 1 = , F8.6, 9H,DELTA 2 ,
* 2H= , F8.6/25X, 10HOMEGA P = , F8.6, 11H,OMEGA S = , F8.6/)
9989 FORMAT (31H EST. MIN. WORDLENGTH WC MIN = , I2, 9H BITS (SI,
* 30HGN+BITS BEHIND THE DUAL POINT)/ 22H ESTIMATED LOWPASS LEN,
* 9HGTH NF = , I3/)
9988 FORMAT (54H REDUCED LOWPASS TOLERANCE SCHEME DELTA 1 PRIME = ,
* F8.6/ 54H DELTA 2 PRIME =
* , F8.6/)
9987 FORMAT (4(1X, F13.6))
9986 FORMAT (1H /50H ACTUAL TOLERANCES OF THE EQUIVALENT LOWPASS WITH ,
* 21HROUNDED COEFFICIENTS /11H DELTA 1 = , F8.6, 11H,DELTA 2 = ,
* F8.6/)
9985 FORMAT (54H EXACT FILTER COEFFICIENTS H0 (I), I = 0...((NF-1)/2)
* /)
9984 FORMAT (1H /37H ROUNDED FILTER COEFFICIENTS H0 (I) /)
9983 FORMAT (31H NUMBER OF DESIGN ITERATIONS = , I3)
9982 FORMAT (34H NO SOLUTION FOUND WITHIN ALLOWED , I2, 11H DESIGN ITE,
* 8HRATIONS-/47H IN THE LAST TRIAL, THE FOLLOWING PARAMETERS WE,
* 9HRE USED /)
9981 FORMAT (1H ///48H INCONSISTENT SPECIFICATIONS,CHECK PARAMETER ORD,
* 6HERING )
9980 FORMAT (1H ///48H FILTER REQUIREMENTS EXCEED MAX. ALLOWED FILTER ,
* 6HORDER )
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: I D E F I R
C
C DESIGN OF LINEAR-PHASE FIR FILTERS IN DIRECT FORM WITH MINIMUM
C COEFFICIENT WORDLENGTH. POSSIBLE FILTER TYPES INCLUDE LOWPASS
C HIGHPASS,SYMMETRICAL BANDPASS AND BANDSTOP FILTERS (WITH CENTER
C FREQUENCY AT OMEGA = 0.25, WHERE OMEGA = 0.5 IS HALF THE SAMPLING
C FREQUENCY), AND HILBERT TRANSFORMERS.
C THE DESIGN IS CARRIED OUT IN THE FOLLOWING WAY - -
C PROBLEM TRANSFORMATION TO AN EQUIVALENT LOWPASS PROBLEM;
C FILTER LENGTH,WORDLENGTH, AND ROUNDING ERROR ESTIMATION;
C LOWPASS DESIGN WITH APPROPRIATELY REDUCED TOLERANCES;
C COEFFICIENT ROUNDING, TOLERANCE CONTROL; ITERATION IF NECESSARY;
C COEFFICIENT TRANSFORMATION BACK TO THE DESIRED FILTER TYPE.
C
C I N P U T
C I N P U T P A R A M E T E R S = TOLERANCE SCHEME WITH
C DEL1...PASSBAND TOLERANCE DELTA 1;
C DEL2...STOPBAND TOLERANCE DELTA 2;
C OMGP1..FIRST PASSBAND EDGE FREQUENCY;
C OMGS1..FIRST STOPBAND EDGE FREQUENCY;
C OMGP2..SECOND PASSBAND EDGE FREQUENCY;
C OMGS2..SECOND STOPBAND EDGE FREQUENCY.
C NORMALIZED FREQUENCY OMG = F / SAMPLING RATE.
C IN LOWPASS AND HIGHPASS CASES,THE LAST TWO PARAMETERS ARE ZERO.
C IN BANDFILTER CASES,THE LAST TWO PARAMETERS HAVE THE VALUES
C OMGP(S)2 = 0.5 - OMGP(S)1
C IN HILBERT CASES, DEL2 = DEL1 IS VALID;OMGP1 IS THE HILBERT BAND
C EDGE FREQUENCY; THE LAST T H R E E PARAMETERS ARE ZERO.
C
C O U T P U T
C O U T P U T P A R A M E T E R S
C DEL1,DEL2,OMGP1,OMGS1...TOLERANCE SCHEME OF THE EQUIVALENT LOWPASS
C PROBLEM; THE ORIGINAL TOLERANCE SCHEME IS REPLACED .
C ITYP...FILTER TYPE INFORMATION- ITYP =0..LOWPASS;
C 1..HIGHPASS;
C 2..BANDPASS;
C 3..BANDSTOP;
C 4..HILBERT-TRANSFORMER;
C NF00...FILTER LENGTH ACCORDING TO TOLERANCES AND ERROR ESTIMATION;
C NF.....FILTER LENGTH ACTUALLY APPLIED;
C IWCMIN...ESTIMATED MIN.COEFF.WORDLENGTH (ACCORDING TO MAX.QUANTI-
C ZATION STEPSIZE + SIGN BIT);
C IWC......WORDLENGTH USED FOR THE DESIGN (IN THE SAME SENSE);
C IWACT....ACTUAL WORDLENGTH (SIGN + NECESSARY BITS BEHIND THE DUAL
C POINT (WITHOUT LEADING ZEROS) FOR THE LARGEST ABS.VALUE);
C DEL11,DEL21...REDUCED TOLERANCES FOR THE EQUIVALENT LOWPASS DESIGN
C EMAX1,2..TOLERANCES (I.E. MAX. ERRORS) OF THE EQUIVALENT LOWPASS
C WITH ROUNDED COEFFICIENTS;
C H0.......COEFFICIENTS BEFORE ROUNDING;
C RH0......ROUNDED COEFFICIENTS;
C ITERA....NUMBER OF DESIGN ITERATIONS;
C IERR.....ERROR PARAMETER -
C IERR= 1...NO ERROR;
C 2...INCONSISTENT SPECIFICATIONS;
C 3...MAX. ALLOWED FILTER ORDER EXCEEDED;
C 4...MAX. ALLOWED NUMBER OF ITERATIONS EXCEEDED
C-----------------------------------------------------------------------
C
SUBROUTINE IDEFIR(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP,
* NF00, NF, IWCMIN, IWC, IWACT, DEL11, DEL21, EMAX1, EMAX2, H0,
* RH0, ITERA, IERR)
C
C COMMON-DEFINITIONS HAVE TO BE INSERTED HERE ACCORDING TO THE SUB-
C ROUTINE 'DESIGN' CALLED BY THIS SUBROUTINE.
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
COMMON /CHEK/ NL, NG, KA, KB
C
DIMENSION H0(256), RH0(128)
IERR = 1
SMALL = 2.*R1MACH(3)
C
C THE FOLLOWING PARAMETERS CONCERN THE MAX. ALLOWED NUMBER OF DESIGN
C ITERATIONS (ITMAX), THE MINIMUM USAGE OF THE TOLERANCES (PROZ),THE
C REDUCTION (CR) AND INCREASE (CI) FACTORS FOR THE TOLERANCE CHANGE.
C
ITMAX = 20
PROZ = 0.9
CR = 0.75
CI = 0.25
C
C TRANSFORMATION TO LOWPASS PROBLEM.
C
CALL TYPTRA(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP, IERR)
C
C IERR = 2 IN THE CASE OF INCONSISTENT SPECIFICATIONS.
C
IF (IERR.EQ.2) RETURN
C
DELOMG = OMGS1 - OMGP1
C
C CHECKING THE FILTER SYMMETRY.
C
OM = ABS(OMGP1+OMGS1-.5)
ISYM = 0
IF ((DEL1.EQ.DEL2) .AND. (OM.LT.SMALL)) ISYM = 1
ITERA = 0
ISOLUT = 0
C
C DETERMINATION OF MINIMUM COEFFICIENT WORDLENGTH.
C
CALL WCMIN(DEL1, DEL2, DELOMG, IWCMIN)
IWC = IWCMIN
C
C O U T E R D E S I G N I T E R A T I O N L O O P
C
10 CONTINUE
DEL11 = DEL1
DEL21 = DEL2
C
C I N N E R D E S I G N I T E R A T I O N L O O P
C
20 CONTINUE
C
C STOP AFTER MAXIMUM NUMBER OF ITERATIONS.
C
IF (ITERA.LT.ITMAX) GO TO 30
IERR = 4
GO TO 80
30 CONTINUE
C
C DETERMINATION OF FILTER LENGTH.
C
CALL LENGTH(DEL11, DEL21, DELOMG, NF)
NN = (NF/2)*2
IF (NN.EQ.NF) NF = NF + 1
C
C IERR = 3 , IF THE MAX. ALLOWED FILTER ORDER IS EXCEEDED.
C
IF (NF.LE.128) GO TO 40
IF ((ITYP.NE.2) .AND. (ITYP.NE.3) .AND. (NF.LE.256)) GO TO 40
IERR = 3
RETURN
40 CONTINUE
C
C EXPECTED EFFICIENT ROUNDING ERROR.
C
IF (ISYM.EQ.1) NF = ((NF+2)/4)*2
CALL EFFERR(DEL11, NF, IWC, S1)
CALL EFFERR(DEL21, NF, IWC, S2)
IF (ISYM.EQ.1) NF = 2*NF - 1
C
C COMPARISON OF ESTIMATED TOTAL ERRORS WITH TOLERANCES.
C
F1 = DEL11 + S1
DD1 = F1 - DEL1
F2 = DEL21 + S2
DD2 = F2 - DEL2
IF ((DD1.GT.SMALL).OR.(DD2.GT.SMALL)) GO TO 70
ITERA = ITERA + 1
IF (ITERA.EQ.1) NF0 = NF
IF (ITERA.EQ.1) NF00 = NF0
C
C EQUIVALENT LOWPASS DESIGN FOR REDUCED TOLERANCES.
C
N1 = (NF+1)/2
N = N1
CALL DESIGN(DEL11, DEL21, OMGP1, OMGS1, H0, NF)
C
C COEFFICIENT ROUNDING TO (IWC-1) BITS BEHIND THE DUAL POINT.
C
Q = 2.**(1-IWC)
DO 50 I=1,N1
HR = H0(I)/Q
RH0(I) = Q*AINT(HR)
50 CONTINUE
C
C DETERMINATION OF REAL MAX. ERRORS.
C
CALL MAXERR(RH0, N1, NF, 0., OMGP1, EMAX1)
CALL MAXERR(RH0, N1, NF, OMGS1, .5, EMAX2)
C
C COMPARISON WITH TOLERANCES.
C
DD1 = EMAX1 - DEL1
DD2 = EMAX2 - DEL2
IF ((DD1.GT.SMALL).OR.(DD2.GT.SMALL)) GO TO 60
ISOLUT = ISOLUT + 1
C
C CHECKING FOR SUFFICIENT USAGE OF TOLERANCES.
C
IF ((DD1/DEL1).GE.(PROZ-1.)) DD1 = 0.
IF ((DD2/DEL2).GE.(PROZ-1.)) DD2 = 0.
IF ((DD1.EQ.0.) .AND. (DD2.EQ.0.)) GO TO 80
IF (NF.LE.NF0) GO TO 80
IF (ISOLUT.EQ.1) NF0 = NF
C
C TOLERANCE CHANGE.
C
60 CONTINUE
IF (DD1.LE.0.) DEL11 = DEL11 - CI*DD1
IF (DD2.LE.0.) DEL21 = DEL21 - CI*DD2
70 CONTINUE
IF (DD1.GT.0.) DEL11 = DEL11 - CR*DD1
IF (DD2.GT.0.) DEL21 = DEL21 - CR*DD2
IF (ISYM.EQ.1) DEL11 = AMIN1(DEL11,DEL21)
IF (ISYM.EQ.1) DEL21 = DEL11
C
C CHECK FOR REALIZABILITY.
C
IF ((DEL11.GT.0.) .AND. (DEL21.GT.0.)) GO TO 20
C
C E N D O F I N N E R D E S I G N I T E R A T I O N L O O P .
C
C
C INCREASE OF WORDLENGTH IN CASE OF UNREALIZABILITY.
C
IWC = IWC + 1
GO TO 10
C
C E N D O F O U T E R D E S I G N I T E R A T I O N L O O P .
C
C
C TRANSFORMATION TO DESIRED FILTER TYPE.
C
80 CONTINUE
IF (ITYP.EQ.0) GO TO 90
CALL BACTRA(ITYP, H0, RH0, NF)
N1 = (NF+1)/2
90 CONTINUE
C
C DETERMINATION OF ACTUAL COEFFICIENT WORDLENGTH.
C
CALL WCACT(RH0, N1, IWC, IWACT)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TYPTRA
C TRANSFORMATION TO EQUIVALENT LOWPASS;SPECIFICATIONS TEST.
C-----------------------------------------------------------------------
C
SUBROUTINE TYPTRA(DEL1, DEL2, OMGP1, OMGS1, OMGP2, OMGS2, ITYP,
* IERR)
C
C TRANSFORMATION TO EQUIVALENT LOWPASS PROBLEM.
C IDENTIFICATION OF FILTER TYPES- ITYP = 0..LOWPASS)
C ITYP = 1..HIGHPASS)
C ITYP = 2..BANDPASS)
C ITYP = 3..BANDSTOP)
C ITYP = 4..HILBERT-TRANSFORMER
C
C L O W P A S S
C
IF ((OMGP2.NE.0.) .OR. (OMGS2.NE.0.)) GO TO 20
IF (OMGP1.GT.OMGS1) GO TO 10
ITYP = 0
GO TO 60
C
C H I G H P A S S
C
10 CONTINUE
IF (OMGS1.EQ.0.) GO TO 40
ITYP = 1
OMGP1 = .5 - OMGP1
OMGS1 = .5 - OMGS1
GO TO 60
20 CONTINUE
C
C CHECKING FOR INCONSISTENT SPECIFICATIONS.
C SET UP A MACHINE CONSTANT.
C
EPSMAC = 2.*R1MACH(3)
DOP = ABS(OMGP1+OMGP2-.5)
DOS = ABS(OMGS1+OMGS2-.5)
DO = AMAX1(DOP,DOS)
IF (DO.GE.EPSMAC) GO TO 50
C
C B A N D P A S S
C
IF (OMGP1.LT.OMGS1) GO TO 30
ITYP = 2
OMGP1 = .5 - 2.*OMGP1
OMGS1 = .5 - 2.*OMGS1
GO TO 60
C
C B A N D S T O P
C
30 CONTINUE
ITYP = 3
OMGP1 = 2.*OMGP1
OMGS1 = 2.*OMGS1
GO TO 60
C
C H I L B E R T - T R A N S F O R M E R
C
40 CONTINUE
C
C CHECKING FOR INCONSISTENT SPECIFICATIONS.
C
IF ((DEL1.NE.DEL2) .AND. (DEL2.NE.0.)) GO TO 50
ITYP = 4
OMGS1 = .25 + OMGP1
OMGP1 = .25 - OMGP1
DEL1 = .5*DEL1
DEL2 = DEL1
GO TO 60
C
C ERROR - INCONSISTENT SPECIFICATIONS.
C
50 IERR = 2
C
60 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: BACTRA
C TRANSFORMATION OF THE DESIGNED LOWPASS TO THE DESIRED FILTER TYPE.
C-----------------------------------------------------------------------
C
SUBROUTINE BACTRA(ITYP, H0, RH0, NF)
C
DIMENSION H0(256), RH0(128)
C
N1 = (NF+1)/2
IF (ITYP-3) 10, 30, 60
10 CONTINUE
C
C H I G H P A S S AND B A N D P A S S
C
N = N1 - 1
DO 20 I=1,N,2
J = N1 - I
H0(J) = -H0(J)
RH0(J) = -RH0(J)
20 CONTINUE
IF (ITYP.EQ.1) GO TO 80
30 CONTINUE
C
C B A N D P A S S AND B A N D S T O P
C
N11 = N1
N1 = NF
NF = 2*NF - 1
NN = 2*N11
I = N1
J = N11
IF (NN-N1) 40, 40, 50
40 CONTINUE
C
C LOWPASS FILTER LENGTH EVEN.
C
H0(I) = 0.
RH0(I) = 0.
I = I - 1
50 CONTINUE
C
C LOWPASS FILTER LENGTH EVEN OR ODD.
C
H0(I) = H0(J)
H0(I-1) = 0.
RH0(I) = RH0(J)
RH0(I-1) = 0.
I = I - 2
J = J - 1
IF (J-1) 80, 80, 50
60 CONTINUE
C
C H I L B E R T - T R A N S F O R M E R
C
I = N1
H0(I) = 0.
RH0(I) = 0.
V = -2.
I = I - 1
70 CONTINUE
V = -V
H0(I) = V*H0(I)
RH0(I) = V*RH0(I)
I = I - 2
IF (I) 80, 80, 70
80 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: WCMIN
C MINIMUM COEFFICIENT WORDLENGTH ESTIMATION BY AN EMPIRICAL FORMULA.
C-----------------------------------------------------------------------
C
SUBROUTINE WCMIN(D1, D2, DOMG, IW)
C
ALD = ALOG10(D1)
ALS = ALOG10(D2)
ALO = ALOG10(DOMG*360.)
W = 7.22116 - 0.589268*ALD
W = W + (0.112354*ALD-2.328539)*ALO
W = W - 3.559051*ALS
IW = IFIX(W)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: LENGTH
C FILTER LENGTH ESTIMATION BY AN EMPIRICAL FORMULA.
C-----------------------------------------------------------------------
C
SUBROUTINE LENGTH(D1, D2, DOMG, NF)
C
DL1 = ALOG10(D1)
DL2 = ALOG10(D2)
DQ1 = DL1*DL1
D = .005309*DQ1 + .07114*DL1 - .4761
D = D*DL2
D = D - .00266*DQ1 - .5941*DL1 - .4278
D = 2. + D/DOMG
NF = IFIX(D)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: EFFERR
C ITERATIVE CALCULATION OF THE APPROXIMATE STATISTICAL ERROR BOUND S
C FOR A GIVEN TOLERANCE D,FILTER LENGTH NF,COEFF.WORDLENGTH IW.
C-----------------------------------------------------------------------
C
SUBROUTINE EFFERR(D, NF, IW, S)
C
QH = 2.**(-IW)
CALL SSTAT(NF, 0., SQ)
ALPHA = D/(SQ*QH)
10 CONTINUE
CALL SSTAT(NF, ALPHA, SQ)
ALPHA1 = D/(SQ*QH)
DALPHA = ALPHA1 - ALPHA
IF (DALPHA-.05) 30, 30, 20
20 CONTINUE
ALPHA = ALPHA1
GO TO 10
30 CONTINUE
S = QH*SQ
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SSTAT
C EVALUATION OF THE APPROXIMATE STATISTICAL ERROR FORMULA.
C-----------------------------------------------------------------------
C
SUBROUTINE SSTAT(NF, ALPHA, SQ)
C
PI = 4.*ATAN(1.0)
PIQ = PI*PI
C
C WEIGHTING FUNCTION M(ALPHA).
C
AM = 1./(1.+SQRT(ALPHA))
C
C ERROR BOUND.
C
ANZ = FLOAT(NF-1)
SQ = (ANZ/PI+.5)*AM
X = (4.-24.*AM/PIQ)*ANZ
X = X + 4. - 3.*AM
X = SQRT(X*AM/3.)
SQ = SQ + X
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: MAXERR
C SEARCH FOR THE MAX.ABS.ERROR OF A LINEAR-PHASE FIR FILTER.
C-----------------------------------------------------------------------
C
SUBROUTINE MAXERR(RH0, N1, NF, OMGA, OMGE, EMAX)
C
DIMENSION RH0(128)
EMAX = 0.
OMG = OMGA
IA = 1
ISTG = 1
DOM1 = 1./(4.*FLOAT(NF))
DOM2 = DOM1/2.
DOM3 = DOM2/2.
DFR = .05
DFRH = .025
CALL ERROR(RH0, N1, NF, OMG, E)
BF1 = ABS(E)
IF (BF1.GT.EMAX) EMAX = BF1
DELOMG = DOM3
OMG = OMG + DELOMG
10 CONTINUE
BF2 = BF1
20 CONTINUE
CALL ERROR(RH0, N1, NF, OMG, E)
BF1 = ABS(E)
IF (BF1.GT.EMAX) EMAX = BF1
IF (IA) 30, 110, 30
30 CONTINUE
DELFR = 2.*(BF1-BF2)/(BF1+BF2)
IF (ISTG) 90, 90, 40
40 CONTINUE
IF (DELFR) 80, 80, 50
50 CONTINUE
ISTG = 1
DELOMG = DOM2
IF (DELFR.LT.DFRH) DELOMG = DOM3
IF (DELFR.GT.DFR) DELOMG = DOM1
OMG = OMG + DELOMG
60 CONTINUE
IF (OMG-OMGE) 10, 70, 70
70 CONTINUE
OMG = OMGE
IA = 0
GO TO 10
80 CONTINUE
ISTG = -1
IF (DELOMG.EQ.DOM1) DELOMG = -3.*DOM3
IF (DELOMG.EQ.DOM2) DELOMG = -DOM3
IF (DELOMG.EQ.DOM3) DELOMG = DOM1
OMG = OMG + DELOMG
IF (DELOMG) 20, 60, 60
90 CONTINUE
IF (DELFR) 100, 100, 50
100 CONTINUE
DELOMG = DOM1
OMG = OMG + DELOMG
GO TO 60
110 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ERROR
C ERROR E OF A LINEAR-PHASE FIR LOWPASS TRANSFER FUNCTION WITH
C RESPECT TO THE IDEAL LOWPASS.
C-----------------------------------------------------------------------
C
SUBROUTINE ERROR(RH0, N1, NF, OMG, E)
C
DIMENSION RH0(128)
C
PI = 4.*ATAN(1.0)
PI2 = PI*2.
NN = 2*N1
IF (NF-NN) 10, 30, 30
10 CONTINUE
C
C LENGTH NF ODD.
C
H = RH0(N1)
N = N1 - 1
DO 20 I=1,N
ANNU = FLOAT(N1-I)
H = H + 2.*RH0(I)*COS(ANNU*OMG*PI2)
20 CONTINUE
GO TO 50
30 CONTINUE
C
C LENGTH NF EVEN.
C
H = 0.
DO 40 I=1,N1
ANNU = FLOAT(N1-I) + .5
H = H + RH0(I)*COS(ANNU*OMG*PI2)
40 CONTINUE
H = H*2.
50 CONTINUE
C
C ERROR VALUE E.
C
E = H
E1 = H - 1.
EA = ABS(E)
EA1 = ABS(E1)
IF (EA1.LT.EA) E = E1
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: WCACT
C DETERMINATION OF ACTUAL WORDLENGTH.
C-----------------------------------------------------------------------
C
SUBROUTINE WCACT(RH0, N1, IWC, IWACT)
C
DIMENSION RH0(128)
HMAX = ABS(RH0(1))
DO 10 I=2,N1
H = ABS(RH0(I))
HMAX = AMAX1(HMAX,H)
10 CONTINUE
L = -2
20 CONTINUE
L = L + 1
IF (HMAX-1.) 30, 40, 40
30 CONTINUE
HMAX = HMAX*2.
GO TO 20
40 CONTINUE
IWACT = IWC - L
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESIGN
C AUXILIARY ROUTINE CALLING THE ODD-LENGTH LOWPASS DESIGN ROUTINE
C 'TOMSB2'.
C-----------------------------------------------------------------------
C
SUBROUTINE DESIGN(DEL11, DEL21, OMGP1, OMGS1, H0, NF)
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
COMMON /CHEK/ NL, NG, KA, KB
DIMENSION H0(256)
C
C AUXILIARY VARIABLES FOR 'TOMSB2'.
C
FP = OMGP1
FS = OMGS1
CALL TOMSB2(DEL11, DEL21, FP, FS, H0, NF)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TOMSB2
C SUBROUTINE ACCORDING TO "PROGRAM TOM" WITH THE
C EXCHANGE ALGORITHM FOR DESIGNING LOW PASS FINITE LENGTH IMPULSE
C RESPONSE DIGITAL FILTERS WITH LINEAR PHASE.
C
C J.H. MCCLELLAN AND T.W. PARKS---RICE UNIVERSITY---OCTOBER 26, 1971
C
C THE FILTER LENGTH (NF), THE PASSBAND CUTOFF FREQUENCY (FP),
C THE STOPBAND FREQUENCY (FS), AND THE WEIGHTING FACTOR (AA) ARE
C INPUTS TO THE ALGORITHM. THEN THE STOPBAND DEVIATION (DS) IS
C MINIMIZED AND THUS THE PASSBAND DEVIATION (DP=AA*DS) IS ALSO
C MINIMIZED. THE PROGRAM SEKS TO DETERMINE THE BEST LOCATION OF
C THE PEAKS OF THE ERROR CURVE IN ORDER TO MINIMIZE THE DEVIATIONS.
C
C REFERENCE--PARKS,T.W. AND J.H. MCCLELLAN,*CHEBYCHEV APPROXIMATION
C OF NON-RECURSIVE DIGITAL FILTERS WITH LINEAR PHASE*, IEEE TRANS.
C ON CIRCUIT THEORY, VOL. CT-19, MARCH, 1972.
C
C SPECIAL PARAMETERS FOR CHANGE OF FS, IF STOPBAND RIPPLE IS
C SPECIFIED
C D1D = STOPBAND RIPPLE D2
C STEP = INITIAL STEPSIZE FOR CHANGING FS& STEP DECREASES WITH
C EXP.& IF SOLUTION FAR FROM OPTIMUM TRY 0.02
C XAC = ACCURACY OF D2
C 1000. MEANS ACCURACY TO 0.1 PERCENT
C 100. MEANS ACCURACY TO 1 PERCENT ETC.
C IF XAC = D1D NO ACTION IS TAKEN AND FS REMAINS UNCHANGED
C
C DIMENSION OF X, AD, Y, IFR, FOPT = N*/2+3 WHERE N* IS
C MAXIMUM LENGTH OF N
C DIMENSION OF AF = 16(N*/2+3)+1
C-----------------------------------------------------------------------
C
SUBROUTINE TOMSB2(DEL1, DEL2, FP, FS, H0, NF)
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
COMMON /CHEK/ NL, NG, KA, KB
DIMENSION H0(256)
C
IWCHN = I1MACH(2)
DO 5 I=1,131
X(I) = 0.
5 CONTINUE
C
C INITIALIZATION OF THE GRID OF 16N POINTS. THE CONSTANT LGRID
C CAN BE CHANGED IF A DIFFERENT GRID DENSITY IS DESIRED.
C
LGRID = 16
XAC = DEL2
AA = DEL1/DEL2
D1D = DEL2
STEP = 0.02
ICH = 1
KITER = 0
D1C = D1D/XAC
10 CONTINUE
PI2 = 8.0*ATAN(1.0)
KCOUNT = 0
TW = FS - FP
NGRID = LGRID*(N-1)
NGRID = NGRID + LGRID/2
KGRID = 2
DELF = 0.5/FLOAT(NGRID)
AF(1) = 0.0
DO 20 J=1,NGRID
AF(J+1) = FLOAT(J)*DELF
20 CONTINUE
KFP = FP/DELF - 1.
KFS = FS/DELF
30 CONTINUE
IF (FP.LT.AF(KFP+1)) GO TO 40
KFP = KFP + 1
GO TO 30
40 CONTINUE
IF (FS.LE.AF(KFS)) GO TO 50
KFS = KFS + 1
GO TO 40
50 CONTINUE
AF(KFP) = FP
AF(KFS) = FS
NGRID = NGRID + 1
DO 60 J=1,NGRID
AF(J) = COS(PI2*AF(J))
60 CONTINUE
C
C INITIAL GUESS FOR THE N+1 OPTIMAL FREQUENCIES
C
KP = (KFP-1)/2
KP = 2*KP + 1
KS = KFS/2
KS = 2*KS + 1
XPP = KP - 1
XSS = NGRID - KS
XAA = XPP + XSS
PT = XAA/FLOAT(N-1)
IF (PT.GT.3.) GO TO 70
WRITE (IWCHN,9999)
9999 FORMAT (20H ERROR IN INPUT DATA)
GO TO 560
70 KR = XPP/PT
IF (KR.EQ.0) KR = KR + 1
PT = XPP/FLOAT(KR)
IFR(1) = 1
DO 80 J=1,KR
XTW = FLOAT(J)*PT
NTW = IFIX(XTW) + 1
NTW = NTW/2
IFR(J+1) = 2*NTW + 1
80 CONTINUE
KR = KR + 1
IFR(KR) = KP
IFR(KR+1) = KS
NTW = N - KR
IF (NTW.GE.0) GO TO 90
WRITE (IWCHN,9999)
GO TO 560
90 IF (NTW.EQ.0) GO TO 110
PT = XSS/FLOAT(NTW)
DO 100 J=1,NTW
XTW = FLOAT(J)*PT + FLOAT(KS) - 1.0
NB = XTW
NB = NB/2
K111 = KR + 1 + J
IFR(K111) = 2*NB - 1
100 CONTINUE
110 CONTINUE
C
C MAIN CALCULATION ROUTINE
C CALCULATE RHO ON THE OPTIMAL SET
C THE FIRST 4 PASSES ARE MADE ON A GRID OF 8N POINTS THEN THE
C FINER GRID OF 16N POINTS IS USED.
C
120 CONTINUE
N111 = N + 1
DO 130 J=1,N111
J111 = IFR(J)
X(J) = AF(J111)
130 CONTINUE
M = (N-1)/15 + 1
M = (N-1)/20 + 1
N111 = N + 1
DO 140 J=1,N111
AD(J) = D(J,N+1,M)
140 CONTINUE
R1 = 0.0
DO 150 J=1,KR
R1 = R1 + AD(J)
150 CONTINUE
R2 = 0.0
K = 1
DO 160 J=1,KR
R2 = R2 + AD(J)*AA*FLOAT(K)
K = -K
160 CONTINUE
K111 = KR + 1
N111 = N + 1
DO 170 J=K111,N111
R2 = R2 + FLOAT(K)*AD(J)
K = -K
170 CONTINUE
RHO = R1/R2
NU = 1
IF (RHO.GT.0.0) NU = -1
RHO = ABS(RHO)
K = NU
XRAY = X(KR+1)
DO 180 J=1,KR
AD(J) = AD(J)*(XRAY-X(J))
Y(J) = 1.0 + AA*FLOAT(K)*RHO
K = -K
180 CONTINUE
K = -K
K111 = KR + 1
DO 190 J=K111,N
X(J) = X(J+1)
AD(J) = AD(J+1)*(XRAY-X(J))
Y(J) = FLOAT(K)*RHO
K = -K
190 CONTINUE
Y(N+1) = FLOAT(K)*RHO
NL = 0
NG = 0
TEST = GEE(IFR(KR+1))
RO = ABS(TEST-RHO)/RHO
IF (RO.LT.0.1 .OR. KCOUNT.LT.4) GO TO 200
WRITE (IWCHN,9998)
9998 FORMAT (20H INTERPOLATION ERROR)
GO TO 560
200 CONTINUE
C
C SEARCH FOR THE NEW EXTREMAL FREQUENCIES
C
NTOT = 0
NPT = IFR(1)
NG = IFR(2)
KA = 1
KB = 1
KSIGN = NU
XT = Y(1)
NFIR = IFR(1)
NOLD = IFR(N+1)
IF (KR.EQ.1) GO TO 330
IF (NPT-1) 210, 260, 210
210 ZT = GEE(NPT-KGRID)
IF (FLOAT(KSIGN)*(ZT-XT)) 260, 260, 220
220 KPT = NPT - KGRID
230 KPT = KPT - KGRID
XT = ZT
IF (KPT-1) 250, 240, 240
240 ZT = GEE(KPT)
IF (FLOAT(KSIGN)*(ZT-XT)) 250, 250, 230
250 IFR(1) = KPT + KGRID
NTOT = NTOT + 1
KSIGN = -KSIGN
GO TO 310
260 ZT = GEE(NPT+KGRID)
IF (FLOAT(KSIGN)*(ZT-XT)) 300, 300, 270
270 KPT = NPT
280 KPT = KPT + KGRID
XT = ZT
ZT = GEE(KPT+KGRID)
IF (FLOAT(KSIGN)*(ZT-XT)) 290, 290, 280
290 IFR(1) = KPT
NTOT = NTOT + 1
300 KSIGN = -KSIGN
310 XFIR = ABS(1.0-XT)/AA
KA = 2
320 IF (KA.GT.(N+1)) GO TO 420
NL = NPT
NG = NGRID + 1
IF (KA.LE.N) NG = IFR(KA+1)
NPT = IFR(KA)
XT = Y(KA)*FLOAT(KB)
IF (KA.EQ.KR) GO TO 330
IF (KA.EQ.(KR+1)) GO TO 340
GO TO 350
330 IF ((KFP-NPT).LT.KGRID) GO TO 410
NPT = KFP
NTOT = 1
GO TO 410
340 KB = -1
IF ((NPT-KFS).LT.KGRID) GO TO 410
NPT = KFS
NTOT = 1
GO TO 410
350 KPT = NPT
ZT = GEE(KPT-KGRID)
IF (FLOAT(KSIGN)*(ZT-XT)) 380, 380, 360
360 XT = ZT
KPT = KPT - KGRID
ZT = GEE(KPT-KGRID)
IF (FLOAT(KSIGN)*(ZT-XT)) 370, 370, 360
370 IFR(KA) = KPT
KSIGN = -KSIGN
NTOT = NTOT + 1
KA = KA + 1
GO TO 320
380 IF ((KPT+KGRID).GT.NGRID) GO TO 410
ZT = GEE(KPT+KGRID)
IF (FLOAT(KSIGN)*(ZT-XT)) 410, 410, 390
390 XT = ZT
KPT = KPT + KGRID
IF ((KGRID+KPT).GT.NGRID) GO TO 400
ZT = GEE(KPT+KGRID)
IF (FLOAT(KSIGN)*(ZT-XT)) 400, 400, 390
400 IFR(KA) = KPT
KSIGN = -KSIGN
NTOT = NTOT + 1
KA = KA + 1
GO TO 320
410 IFR(KA) = NPT
KSIGN = -KSIGN
KA = KA + 1
GO TO 320
420 CONTINUE
XOLD = ABS(XT)*AA
NL = 0
NG = 0
IF (IFR(1).EQ.1) GO TO 450
IF (IFR(N+1).EQ.NGRID) GO TO 430
WRITE (IWCHN,9997)
9997 FORMAT (38H ERROR-NEITHER ENDPOINT IS AN EXTREMUM)
GO TO 560
430 IF (NFIR.EQ.1) GO TO 470
TEST = GEE(1)
TEST = -FLOAT(NU)*(TEST-1.0) - XOLD
IF (TEST.LE.0.0) GO TO 470
NTOT = 1
DO 440 J=1,N
N111 = N + 2 - J
M111 = N + 1 - J
IFR(N111) = IFR(M111)
440 CONTINUE
IFR(1) = 1
KR = KR + 1
GO TO 470
450 IF (IFR(N+1).EQ.NGRID) GO TO 470
IF (NOLD.EQ.NGRID) GO TO 470
TEST = GEE(NGRID)
TEST = FLOAT(KSIGN)*TEST - XFIR
IF (TEST.LE.0.0) GO TO 470
NTOT = 1
DO 460 J=1,N
IFR(J) = IFR(J+1)
460 CONTINUE
IFR(N+1) = NGRID
KR = KR - 1
470 CONTINUE
IF (KR.GT.0) GO TO 480
WRITE (IWCHN,9996)
9996 FORMAT (6H ERROR)
GO TO 560
480 CONTINUE
KCOUNT = KCOUNT + 1
IF (KCOUNT.GT.50) GO TO 500
IF (KCOUNT.EQ.4 .AND. KGRID.EQ.2) GO TO 490
IF (NTOT.NE.0) GO TO 120
IF (KGRID.EQ.2) GO TO 490
GO TO 500
490 KGRID = 1
GO TO 120
500 CONTINUE
IF (ABS(RHO-D1D).LT.D1C) GO TO 520
KITER = KITER + 1
IF (KITER.GT.1000) GO TO 520
IF (RHO.GT.D1D) GO TO 510
IF (ICH.EQ.1) STEP = STEP/2.
FS = FS - STEP
ICH = -1
GO TO 10
510 IF (ICH.EQ.(-1)) STEP = STEP/2.
FS = FS + STEP
ICH = 1
GO TO 10
520 CONTINUE
KRP = N - KR
R1 = AA*RHO
J = N
530 CONTINUE
IF (X(J)) 540, 550, 540
540 J = J + 1
GO TO 530
550 X(J) = -1.
CALL PREF2(H0, NF)
GO TO 570
560 CONTINUE
STOP
570 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: D
C THE SUBROUTINE D(,,,) CALCULATES THE LAGRANGIAN INTERPOLATION
C COEFFICIENTS FOR USE IN THE BARYCENTRIC INTERPOLATION FORMULA
C-----------------------------------------------------------------------
C
FUNCTION D(K, N, M)
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
D = 1.0
Q = X(K)
DO 30 L=1,M
DO 20 J=L,N,M
IF (J-K) 10, 20, 10
10 D = 2.0*D*(X(J)-Q)
20 CONTINUE
30 CONTINUE
D = 1.0/D
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: GEE
C THE SUBROUTINE GEE EVALUATES THE FREQUENCY RESPONSE USING THE
C INTERPOLATION COEFFICIENTS WHICH ARE CALCULATED IN THE ROUTINE D.
C-----------------------------------------------------------------------
C
FUNCTION GEE(K)
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
COMMON /CHEK/ NL, NG, KA, KB
IF (K-NL) 10, 40, 10
10 IF (K-NG) 20, 40, 20
20 P = 0.0
XF = AF(K)
D = 0.0
DO 30 J=1,N
C = AD(J)/(XF-X(J))
D = D + C
P = P + C*Y(J)
30 CONTINUE
GEE = P/D
RETURN
40 K111 = KA + KB
GEE = Y(K111)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: PREF2
C THE SUBROUTINE PREF SAMPLES THE FREQUENCY RESPONSE AT 2**M POINTS
C SO THAT THE FFT (R.C. SINGLETON) CAN BE CALLED TO INVERSE TRANSFORM
C THIS DATA TO OBTAIN THE IMPULSE RESPONSE.
C-----------------------------------------------------------------------
C
SUBROUTINE PREF2(A, NF)
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
COMPLEX Z
DIMENSION A(256), B(256)
DIMENSION Z(256)
EQUIVALENCE (B(1),AF(257)), (Z(1),AF(520))
C
C SET UP A MACHINE CONSTANT.
C
FSH = R1MACH(1)
C
PI2 = 8.0*ATAN(1.0)
NX = 2
MX = 1
10 IF (NX.GE.NF) GO TO 20
NX = 2*NX
MX = MX + 1
GO TO 10
20 NN = NX
MM = MX
NX = NX/2
XNN = NN
DF = 1.0/XNN
IF (IFR(1).EQ.1) GO TO 30
P = GOO(1.0)
GO TO 40
30 P = Y(1)
40 A(1) = P
IF (IFR(N+1).EQ.NGRID) GO TO 50
P = GOO(-1.0)
GO TO 60
50 P = Y(N)
60 A(NX+1) = P
L = 1
N111 = NX - 1
DO 110 J=1,N111
AT = DF*FLOAT(J)
AT = COS(PI2*AT)
70 AS = X(L)
IF (AT.GT.AS) GO TO 90
IF ((AS-AT).LT.FSH) GO TO 80
L = L + 1
GO TO 70
80 A(J+1) = Y(L)
GO TO 100
90 IF ((AT-AS).LT.FSH) GO TO 80
A(J+1) = GOO(AT)
100 L111 = NN + 1 - J
A(L111) = A(J+1)
B(J+1) = 0.0
IF (L.GT.1) L = L - 1
M111 = NN + 1 - J
B(M111) = 0.
110 CONTINUE
B(1) = 0.0
NX1 = NX + 1
B(NX1) = 0.0
DO 120 J=1,NN
Z(J) = CMPLX(A(J),B(J))
120 CONTINUE
CALL BITREV(Z, MM)
CALL FFTFT(Z, MM)
DO 130 J=1,NX
B(J) = REAL(Z(J))
130 CONTINUE
DO 140 J=1,N
JJ = N + 1 - J
A(JJ) = B(J)
140 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: GOO
C THE SUBROUTINE GOO IS THE SAME AS GEE EXCEPT THAT IT IS CALLED BY
C A FREQUENCY VALUE RATHER THAN BY BY A GRID INDEX VALUE.
C-----------------------------------------------------------------------
C
FUNCTION GOO(F)
C
COMMON /CALC/ X(131), AF(2096), AD(131), Y(131), SCALE
COMMON /CHEB/ N, NGRID, KGRID, KFP, KFS, IFR(131)
XF = F
P = 0.0
D = 0.0
DO 20 J=1,N
IF (XF.NE.X(J)) GO TO 10
GOO = Y(J)
GO TO 30
10 C = AD(J)/(XF-X(J))
D = D + C
P = P + C*Y(J)
20 CONTINUE
GOO = P/D
30 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTFT
C FREQUENCY-TO-TIME FFT ROUTINE; INPLACE ALGORITHM; COMPLEX DATA IN
C ARRAY A; INPUT DATA ASSUMED TO BE IN BITREVERSED ORDER.
C-----------------------------------------------------------------------
C
SUBROUTINE FFTFT(A, M)
C
COMPLEX A, U, W, T
DIMENSION A(1)
PI = 4.0*ATAN(1.)
N = 2**M
DO 30 L=1,M
LE = 2**L
LE1 = LE/2
U = (1.0,0.0)
W = CMPLX(COS(PI/FLOAT(LE1)),SIN(PI/FLOAT(LE1)))
DO 20 J=1,LE1
DO 10 I=J,N,LE
IP = I + LE1
T = A(IP)*U
A(IP) = A(I) - T
A(I) = A(I) + T
10 CONTINUE
U = U*W
20 CONTINUE
30 CONTINUE
DO 40 K=1,N
A(K) = A(K)/FLOAT(N)
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: BITREV
C REORDERING OF 2**M DATA IN COMPLEX ARRAY A IN BITREVERSED ORDER.
C-----------------------------------------------------------------------
C
SUBROUTINE BITREV(A, M)
C
DIMENSION A(1)
COMPLEX A, T
C
N = 2**M
NV2 = N/2
NM1 = N - 1
J = 1
DO 40 I=1,NM1
IF (I.GE.J) GO TO 10
T = A(J)
A(J) = A(I)
A(I) = T
10 K = NV2
20 IF (K.GE.J) GO TO 30
J = J - K
K = K/2
GO TO 20
30 J = J + K
40 CONTINUE
RETURN
END
Return to Main Software Page